Provide tropospheric surface area density calculations to musica chemistry

[K20shores]

CAM-SIMA, through the abstract aerosol interface, now supports some aerosol properties. Those aerosol properties can be used to drive some gas-phase aerosol reactions in musica. This issue details a plan to make those connections at runtime for tropospheric surface area density calculations.

• CAM-SIMA can provide the surface area density and effective radius for a particular list of aerosol species types
  • It provides them for a list of species types, aggregated over modes, with the surface of each internally mixed mode partitioned by the mass fraction of the listed types (modal_aerosol_state_mod.F90:848-1010).
• MUSICA needs the number concentration and effective particle radius for a particular list of aerosol species types
  • MUSICA needs these values to be set in our rate parameters using the reaction name for our surface reaction type
  • the N we back out is not the mode's particle number concentration. It is the number of particles of radius d_m/2 that would carry the mode's reactive surface, so it already includes the lognormal surface enhancement exp(2 ln^2 sigma_g) and the mass fraction of the mapped species types. Call it a "surface-equivalent number" so nobody confuses it with a physical number concentration
• Therefore we need
  • a mapping from a list of CAM-SIMA aerosol species types to a MUSICA surface reaction
    • Or, we map the CAM-SIMA aerosol species to species in MUSICA. From there, we can map the musica species to the particular reactions, which will then allow us to calculate and set the aerosol information we need. Either one works, as long as the aerosol side stays a list of species types and not a single species

Acceptance criteria

• Configured CAM-SIMA aerosol species can be mapped to configured MUSICA surface reactions
• Properly mapped species have their aerosol properties for surface area density and effective radius either used in the musica ccpp wrapper to set the number concentration and effective radius for each musica surface reaction at each chemistry timestep
• Tests are added which have valid gas-phase surface reactions that match CAM-SIMA aerosols for MAM
  • CARMA is not yet implemented. BAM is not the target of this work. Focus on MAM
  • mechanism-parameter check. Diff the mechanism configuration against the constants in CAM's mo_usrrxt: reaction probability per reaction equals CAM's gamma_* constant; diffusion coeff is 1e-5 m2/s for every reacting gas (CAM's dg = 0.1 cm2/s), except effectively infinite for glyoxal, whose CAM rate has no diffusion term; molecular weight matches the value behind CAM's hard-coded mean-speed coefficient for that gas
  • summed-rate check. Given a MAM4 state, the sum over modes of MICM's rates should be close enough to CAM's hetrxtrate to 1e-3 relative (mo_usrrxt.F90:3215-3231). Exact agreement is not possible: CAM rounds its mean-speed coefficients (for example 1.40e3 for N2O5) to three digits, which limits agreement to about 1e-4 on its own
    • Consider saving the current rates calculated in CAM to a file (describe exactly how they were obtained, which commit, etc.) and add a regression test against these saved values.
• Tests are added which catch missing/invalid/incorrect mappings of CAM-SIMA aerosols and musica mechanisms

Ideas

• Create one MICM SURFACE reaction per mode per gas. Create one configuration file for each mam configuration
  • naming convention: <base_reaction_name>.mode<l> for mode l of nbins (for example N2O5_sulfate.mode1 ... .mode4), all sharing the same gas-phase reactant/products, reaction probability, and species-type list — only (N_l, r_l) differs per mode
  • we can use the mechanism configuration file list format to reduce the amount of duplication. Place each surface reaction in a file for each mam configuration and share the same core parts of TS1
  • no special handling needed on the MICM side to sum these back together, multiple reactions consuming the same gas-phase species already add in the solver, same as CAM's sum over modes does
• Set the gas-phase diffusion coefficient in the mechanism/rate-constant config to 1e-5 m2 s-1
  (matching CAM's constant D_g = 0.1 cm2/s, used for every species) for parity with CAM, unless a
  species-specific, physically-based value is deliberately chosen instead
• Calculate the number concentration from the surface area density $SAD=4\pi r_{eff}^2 N$, where $N$ is the number concentration
  • Do this per mode, not on the aggregate sad/reff outputs. CAM's mo_usrrxt/hetrxtrate never uses the aggregate sad/reff. Those are diagnostic-only in CAM (SAD_TROP, REFF_TROP history fields). The real rate sums per mode, using the optional per-mode sfc(:,:,l) (cm2/cm3) and dm_aer(:,:,l) (cm) outputs of surf_area_dens. Per mode: r_l = dm_l / 2 and N_l = sfc_l / (pi * dm_l^2) reproduces CAM's per-mode term exactly.
  • note surf_area_dens skips the primary_carbon mode entirely, so black-c only shows up through the accumulation mode
  • make sure to get the units right. CAM-SIMA provides surface area density in [cm2 cm-3] and effective radius in [cm]. Musica needs the number concentration in [particles m-3], and the effective radius in [m]
• Create a mapping for CAM-SIMA aerosols and musica reactions
  • consider implementing a similar mapping to what we do for tuvx
    • we create a mapping file stored alongside our mechanism and tuvx data and then read that file to create a mapping from tuvx to musica
  • we could do something similar, either in the namelist for cam or just another json file since we already have APIs to read these, that map CAM-SIMA MAM aerosol species names to the corresponding surface reaction name in a particular mechanism they should be mapped to
  • default the mapping to the two groupings CAM already uses today, for parity: sad_chem_spec_types (default sulfate, s-organic, p-organic, black-c, ammonium) for the general tropospheric heterogeneous reactions, and sad_seasalt_spec_types (default seasalt) for the sea-salt-specific ones (dehalogenation, ClNO2 yield). both are CAM aerosol_nl namelist variables today and are runtime configurable, we should keep that
  • each mapping entry should map to a list of CAM-SIMA aerosol species types, not a single species name. aero_state%surf_area_dens's types_list argument (aerosol_state_mod.F90:340) is itself an array. A MAM mode is internally mixed, and modal_aerosol_state_mod.F90 computes surface area as the mass-fraction-weighted share attributable to whichever species types are listed. Some surface reactions may need more than one species type combined into one aggregate surface (for example sulfate and nitrate together), so the mapping schema should be reaction_name -> [species_type, ...] from the start, for example {"N2O5_sulfate": ["sulfate"]}, rather than {"N2O5_sulfate": "sulfate"}, which cannot express the combined case later without a schema change
• Add a new file inside the musica tree, for example atmospheric_physics/schemes/musica/aerosol/musica_ccpp_surface_area_density.F90.
  • Use the aerosol interface to fetch the MAM aero_state / aero_props, similar to what ndrop does.
  • Alternatively, select and cache the MAM aerosol-model index once, at init, the way nucleate_ice_ccpp_init does: then loop aerosol_instances_get_num_models() once and keep the matching index as module state.
  • At each chemistry timestep, use the cached index for a direct aerosol_instances_get_props/aerosol_instances_get_state lookup (no search), the way nucleate_ice_ccpp_run does at lines 377-379
  • call surf_area_dens once per species-type list, not once per label. Each call loops every column, level, mode, and species, and allocates work arrays. CAM makes only two calls per timestep total (one per list). TS1 has about 40 uptake reactions times four modes, so 160 labels. Calling once per label costs about 80x more than grouping the labels by their shared species-type list and filling every label in a group from one call
  • make this a separate CCPP scheme, not an internal call from musica_ccpp_run. musica_ccpp is now split into musica_ccpp_photolysis and musica_ccpp_chemistry (atmospheric_physics#443). Place the new scheme between them, taking micm_rate_parameters as inout. The slot is already earmarked as Slot for a SAD -> MICM rate-parameter provider: in test/test_suites/suite_modal_aero_wet_musica.xml, right above mam_surfarea_diagnostics
  • validate the mode count at init. For each base reaction, check that the number of .mode<l> labels equals aero_props%nbins(), and abort if not. A silent mode-count mismatch would otherwise only show up as a wrong answer at run time
  • do not special-case the primary_carbon mode out of the mechanism config. surf_area_dens returns zero surface for that mode, so its .mode<l> reaction just receives N = 0. Generating it and feeding zero is correct, and simpler than excluding it
  • mapping the sea-salt species type to the sea-salt surface is necessary, but it is not sufficient for the sea-salt-specific reactions. The ClNO2 yield depends on total sea-salt mass, with an HCl-limited branch. The dehalogenation rates carry land masks and a dehalogenation factor. Defer this extra physics, but say so explicitly, do not imply the mapping alone finishes these reactions
• A MAM sandbox for this work already exists: jimmielin/CAM-SIMA branch hplin/modal_aero_sad_sandbox (clone, then run bin/git-fleximod update). Its test suite modal_aero_wet_musica (atmospheric_physics/test/test_suites/) runs the full MAM wet-process block plus mam_surfarea_diagnostics, which already calls aero_state%surf_area_dens the way CAM does, writes SAD_TROP, REFF_TROP, and per-mode sad_a01..04/dmaer_a01..04 to the history tape, and runs the identity check from the acceptance criteria below. It has been verified to run to completion on Derecho with wet-process bit-for-bit checks intact
• Below is a list of relevant PRs/files that could be helpful for this work or provide more context

Relevant PRs

• cam6_4_195: Aerosol surface area densities refactor CAM#1628
• Various fixes for MUSICA; split TUV-x and MICM into separate CCPP schemes #443 (splits musica_ccpp into musica_ccpp_photolysis and musica_ccpp_chemistry, the slot for the new SAD scheme)

Relevant code locations

In cam-sima

• cime_config/namelist_definition_cam.xml
  • this namelist has the definition of species that will be included for different aerosols
  • this particular section seems to show that we are defining these BAM species
    • sulf, dust1, dust2, dust3, dust4, bcar1, bcar2, ocar1, ocar2, sslt1, sslt2, sslt3, sslt4
    • when a configuration like this is present, we will need to connect them to species in musica and get the properties we need for our surface reactions
    • this means we will need a way to connect names configured in musica configurations with names configured in the nameilst. We might do this similar to how we connect rates in tuvx with rates in micm
• cime_config/namelist_definition_cam.xml#L430-L441
  • mode_defs is the MAM equivalent of the BAM species list above, but it defaults to empty ('') in cam-sima right now. there's no shipped MAM mode/species example in this repo yet, we'll need to write our own for testing (something like the mam4 mode_defs/rad_climate example from classic CAM would work as a starting point)
• test/unit/fortran/src/aerosol/modal_test_helpers.F90
  • a small hand-built 2-mode mock (accum: sulfate + black-c, coarse: dust) that our own modal aerosol unit tests already use. good starting point for testing this too
• src/aerosol/radiative_aerosol_definitions.F90
  • This shows that there are only a couple of allowed species types in cam-sima.
• src/aerosol/aerosol_state_mod.F90
  • specifically the surface area density and effective radius attributes
• src/aerosol/modal_aerosol_state_mod.F90
  • mam implementaiton
• src/aerosol/aerosol_instances_mod.F90
  • allows us to get the aerosol state
• src/aerosol/aerosol_spec_utils.F90
  • potentially useful for finding a specific species maybe?

In atmospheric physics

• schemes/musica/musica_ccpp.F90
• schemes/musica/micm/musica_ccpp_micm.F90

In musica

• micm state
  • rate parameters
  • rate parameter ordering
  • You can use these to set a particular rate paramter
    • getting an index
    • setting a rate for an index
    • If you happen to make an array that contains every rate parameter correctly ordered, you can pass it directly to musica, state%rate_parameters(:) = my_rate_values(:)
      • because we can solve multiple grid cells simultaneously, this would need to contain all the different rate parameters across all grid cells, corectly ordered

[jimmielin]

Thanks @K20shores. I'm also posting the acceptance criteria previously outlined in our discussion in the original issue NCAR/musica#1030 (comment)

Proposed acceptance criteria

What has to be equal

MICM computes each surface reaction from two numbers per grid cell that the host supplies: a radius r and a number concentration N.

CAM computes the same first-order rate from each aerosol mode's surface area density A_m (the part attributable to the mapped species types; i.e., sfc(:,:,m) from aero_state%surf_area_dens) and the mode's wet number-median diameter d_m (dm_aer(:,:,m)):

CAM,  mode m:     
     k_m = A_m           / ( (d_m / 2) / D_g + 4 / (v gamma) )

MICM, one reaction: 
     k   = 4 pi r^2 N    / (  r        / D_g + 4 / (v gamma) )

D_g is the gas diffusivity, v the mean molecular speed of the gas, and
gamma the reaction probability. The two expressions are the same function.
With one MICM surface reaction per mode, the per-mode identity

     r = d_m / 2      
     N = A_m / (4 pi r^2)

makes MICM reproduce CAM's k_m exactly, and the sum over modes reproduces CAM's hetrxtrate.

n.b. N is not the particle number concentration of the mode.
It is the number of particles of radius d_m/2 that would carry the mode's reactive surface, so it includes the lognormal surface enhancement exp(2 ln^2 sigma_g) and the mass fraction of the mapped species types. It is better called a "surface-equivalent number".

1. Identity test (built into mam_surfarea_diagnostics for your convenience)

For every surface reaction in the mechanism, read back the r and N that
were written into MICM's rate parameters and require, at every grid cell:

• 4 pi r^2 N equals A_m for exactly one mode m, to tolerance,
  after unit conversion
  (MICM uses m and m-3; the aerosol interface returns cm2/cm3 and cm);
• r equals d_m / 2 for that same mode, to tolerance;
• both are zero above the tropopause. -- a tropopause mask is needed here as CAM does not compute trop SADs above the tropopause. tropopause_vertical_layer_index is available as a ccpp input

(here tolerance is 1e-10 relative)

A unit test can be developed with synthetic per-mode inputs.

An integration test is in the mam_surfarea_diagnostics, which writes per surface reaction, the recovered 4 pi r^2 N and the index of the MAM mode it matches, with 0 for "both zero" and -1 for "no match").

Reactions deliberately mapped to other species types (seasalt) will show -1 against CAM's chemistry-type surface.

2. the mechanism parameters

The remaining inputs to k live in the mechanism configuration. Diff them against the constants in CAM's mo_usrrxt:

• reaction probability per reaction equals CAM's gamma_* constant;
• diffusion coeff is 1e-5 m2/s for every reacting gas (CAM's single
  dg = 0.1 cm2/s), and effectively infinite for glyoxal, whose CAM rate
  has no diffusion term;
• molecular weight matches the value behind CAM's hard-coded mean-speed
  coefficient for that gas.

3. the rate itself

The sum over modes of MICM's rates equals CAM's hetrxtrate to 1e-3 relative. CAM's mean-speed coefficients (for example 1.40e3 for N2O5) are rounded to three digits, so exact agreement is impossible by about 1e-4.