Update TS1 config file

[K20shores]

Closes #1041

I converted the music box configuration listed in the issue with the musica-cli. To do so successfully, I had to add three species to the species list.

Then, I worked with claude to add in the 'unsupported reactions' section of the configuration file from chemistry cafe. I added those if they could have real values, or with types that let us set their rate to zero so we can figure out how to handle them later. Below lists what each reaction is, how it's handled, and what we might need to do

Correctly implemented, I think:

• usr_N2O5_aer, usr_NO3_aer, usr_NO2_aer, usr_HO2_aer, usr_GLYOXAL_aer, usr_ISOPNITA_aer, usr_ISOPNITB_aer, usr_ONITR_aer, usr_HONITR_aer, usr_TERPNIT_aer, usr_NTERPOOH_aer, usr_NC4CHO_aer, usr_NC4CH2OH_aer
  • Tropospheric aerosol uptake, one gas-phase reactant each
  • Implemented as SURFACE with constant reaction probabilities from CAM
• het7, het8, het11
  • N2O5, ClONO2, BrONO2 hydrolysis on NAT polar stratospheric cloud particles
  • Implemented as SURFACE with constant reaction probabilities from CAM
• het12, het13, het14
  • Same three reactions, on ice polar stratospheric cloud particles
  • Implemented as SURFACE with constant reaction probabilities from CAM

Deferred for later, but still included

• het1–3
  • The uptake rate depends on temperature and aerosol composition
  • Implementing CAM's composition/temperature-dependent formula would fix this properly.
    • Or, if these can be computed using the surface area density/effective radius provided to use from CAM/CAM-SIMA, we should do that instead
  • Implemented as SURFACE with reaction probability 0
• het4–6
  • MICM's SURFACE type only allows one gas-phase reactant. These need 2.
  • Fixing this needs both the composition-dependent formula above and a new reaction type that can consume two reactants.
  • Implemented as USER_DEFINED with scaling factor 0
• het9, het10
  • Same two-reactant limitation
  • Implemented as USER_DEFINED with scaling factor 0
• het15–17
  • Same situation as het9-10
  • Implemented as USER_DEFINED with scaling factor 0
• usr_CO_OH: CO + OH → CO2 + HO2. JPL's official rate law is two coupled pieces.
  • I guess we used to be able to represent this before JPL19 as the sum of two reactions
  • The updated formula in JPL 19 prevents that because it's a Troe term, plus a correction term that subtracts a fraction of that same Troe term
  • no existing MICM reaction type can compute that subtraction.
  • It's USER_DEFINED with scaling factor 0 for now, and needs a new MICM reaction type, or somehow we need to represent this as a mix of several existing types if it can be
  • Implemented as USER_DEFINED with scaling factor 0

AI usage: Claude helped me correctly set the names of photolysis rates and to either implement the het/user define rates or add stubs

[codecov-commenter]

Codecov Report

❌ Patch coverage is 81.25000% with 3 lines in your changes missing coverage. Please review.
✅ Project coverage is 72.12%. Comparing base (05997c6) to head (a5896d9).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
.../examples/miem_cams_finn_box_model_real_fixture.py	0.00%	1 Missing ⚠️
python/musica/examples/miem_nox_box_model.py	0.00%	1 Missing ⚠️
...musica/examples/miem_nox_box_model_real_fixture.py	0.00%	1 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #1044   +/-   ##
=======================================
  Coverage   72.12%   72.12%           
=======================================
  Files         147      147           
  Lines       12178    12185    +7     
=======================================
+ Hits         8783     8789    +6     
- Misses       3395     3396    +1     

Flag	Coverage Δ
cpp_fortran	65.06% <ø> (ø)
javascript	94.73% <ø> (ø)
python	72.23% <81.25%> (+0.01%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.
• 📦 JS Bundle Analysis: Save yourself from yourself by tracking and limiting bundle sizes in JS merges.

[jimmielin]

Thanks @K20shores for sending this over.

A few general chemistry notes that stood out to me then specifically for CAM-SIMA.

1. usr_SO3_H2O

musica/configs/v1/ts1/t1s2.json

Lines 3985 to 4013 in 10d504b

	{
	"type": "ARRHENIUS",
	"A": 5.1188196459999996e-23,
	"B": 0.0,
	"C": 6539.99676963515,
	"D": 300.0,
	"E": 0.0,
	"reactants": [
	{
	"name": "H2O",
	"coefficient": 1.0
	},
	{
	"name": "SO3",
	"coefficient": 1.0
	}
	],
	"products": [
	{
	"name": "H2SO4",
	"coefficient": 1.0
	},
	{
	"name": "H2O",
	"coefficient": 1.0
	}
	],
	"gas phase": "gas"
	},

I think the A here is mechanistically wrong. It appears that 5.1188e-23 was arrived here by 8.5e-41 * 6.022e17 from CAM

!-----------------------------------------------------------------
!       ... so3 + 2*h2o --> h2so4 + h2o
!       Note: this reaction proceeds by the 2 intermediate steps below
!           so3 + h2o --> adduct
!           adduct + h2o --> h2so4 + h2o
!               (Lovejoy et al., JCP, pp. 19911-19916, 1996)
!	The first order rate constant used here is recommended by JPL 2011.
!	This rate involves the water vapor number density.
!-----------------------------------------------------------------

       if ( usr_SO3_H2O_ndx > 0 ) then
          call comp_exp( exp_fac, 6540.0_r8*tinv(:), ncol )
          if( h2o_ndx > 0 ) then
             fc(:) = 8.5e-21_r8 * m(:,k) * h2ovmr(:,k) * exp_fac(:)
          else
             fc(:) = 8.5e-21_r8 * invariants(:,k,inv_h2o_ndx) * exp_fac(:)
          end if
          rxt(:,k,usr_SO3_H2O_ndx) = 1.0e-20_r8 * fc(:)
       end if

where 6.022e17 is the conversion factor for cm3 molec-1 s-1 to (mol m-3)-1 s-1 ("m") for a two-body rxn constant but the actual mechanics of this reaction is a fast pre-equilibrium

SO3 + H2O <-> SO3.H2O (adduct fast pre-equilibrium)
SO3.H2O + H2O -> H2SO4 + H2O

the steady state observed rate is k(T) [SO3] [H2O]^2, where k(T) is the fc(:) here. Note how it depends on [H2O]^2, but one of the [H2O] factors got entirely dropped in the MUSICA mechanism if you write it out as just SO3 + H2O -> H2SO4 + H2O (as in the .json, the "2" also got dropped from the reactants and needs to be put back)... we can't just multiply the constants, leaving out one of the H2Os is off by a factor of [H2O]... so about 17 orders of magnitude.

2. usr_CO_OH being stubbed out makes the mechanism basically unrealistic as written because CO + OH is the dominant global OH sink (and the only CO oxidation path in TS1, I think?) -- so all the tropospheric oxidants are meaningless until this is fixed.

Either a new MICM type has to be made or fall back to the pre-JPL19 two-term representation. Zero will really mess up the mechanism...

3. the het stubs are fine but this means no strat het chemistry and no ozone-hole chemistry.

More notes on SIMA forthcoming...

[jimmielin]

For CAM-SIMA implementation...

Right now the config has one SURFACE reaction per gas, so there's only one pair of (radius, number). This is at odds with what we agreed on for #1030 where we can't collapse per-mode SAD into one aggregate, so if you're not modifying MICM internals, per-mode duplicates of the same reactions needs to be made to match CAM.

Thinking more about the het.. reactions I realize that the het1-3 being deferred might be problematic, since N2O5 hydrolysis is a dominant pathway to produce HNO3 globally.

I think het1-17 would be better wired up as USER_DEFINED. For het1-3, have the host supply k which is a function of SAD, so you don't have to have variable gamma for SURFACE reactions. For het4 it's trickier as there's a limiting reagent construction

!-----------------------------------------------------------------------
!     	... ClONO2 + HCl(liq) =  Cl2  + HNO3  Sulfate Aerosol Reaction
!-----------------------------------------------------------------------
                  if( hclvmr > small_div .and. clono2vmr > small_div ) then
                     if ( hclvmr > clono2vmr ) then
                        rxt(i,k,rid_het4) = max( 0._r8,wrk*av_clono2*gprob_cnt_hcl(i,k) )*hcldeni
                     else
                        rxt(i,k,rid_het4) = max( 0._r8,wrk*av_clono2*gprob_cnt_hcl(i,k) )*cntdeni
                     end if
                  end if

For het13, 14 I see it's implemented now but my hunch is that MICM's rates may be too slow for ice PSCs. MICM's SURFACE formula is the full resistor but for CAM stratospheric het13, 14, it is the free-molecular construction with no gas-diffusion term. This differs from the trop rates in TS1 (e.g., usr_N2O5_aer...)

[K20shores]

Thanks @K20shores for sending this over.

A few general chemistry notes that stood out to me then specifically for CAM-SIMA.

1. usr_SO3_H2O

musica/configs/v1/ts1/t1s2.json

Lines 3985 to 4013 in 10d504b

	{
	"type": "ARRHENIUS",
	"A": 5.1188196459999996e-23,
	"B": 0.0,
	"C": 6539.99676963515,
	"D": 300.0,
	"E": 0.0,
	"reactants": [
	{
	"name": "H2O",
	"coefficient": 1.0
	},
	{
	"name": "SO3",
	"coefficient": 1.0
	}
	],
	"products": [
	{
	"name": "H2SO4",
	"coefficient": 1.0
	},
	{
	"name": "H2O",
	"coefficient": 1.0
	}
	],
	"gas phase": "gas"
	},

I think the A here is mechanistically wrong. It appears that 5.1188e-23 was arrived here by 8.5e-41 * 6.022e17 from CAM

!-----------------------------------------------------------------
!       ... so3 + 2*h2o --> h2so4 + h2o
!       Note: this reaction proceeds by the 2 intermediate steps below
!           so3 + h2o --> adduct
!           adduct + h2o --> h2so4 + h2o
!               (Lovejoy et al., JCP, pp. 19911-19916, 1996)
!	The first order rate constant used here is recommended by JPL 2011.
!	This rate involves the water vapor number density.
!-----------------------------------------------------------------

       if ( usr_SO3_H2O_ndx > 0 ) then
          call comp_exp( exp_fac, 6540.0_r8*tinv(:), ncol )
          if( h2o_ndx > 0 ) then
             fc(:) = 8.5e-21_r8 * m(:,k) * h2ovmr(:,k) * exp_fac(:)
          else
             fc(:) = 8.5e-21_r8 * invariants(:,k,inv_h2o_ndx) * exp_fac(:)
          end if
          rxt(:,k,usr_SO3_H2O_ndx) = 1.0e-20_r8 * fc(:)
       end if

where 6.022e17 is the conversion factor for cm3 molec-1 s-1 to (mol m-3)-1 s-1 ("m") for a two-body rxn constant but the actual mechanics of this reaction is a fast pre-equilibrium

SO3 + H2O <-> SO3.H2O (adduct fast pre-equilibrium) SO3.H2O + H2O -> H2SO4 + H2O

the steady state observed rate is k(T) [SO3] [H2O]^2, where k(T) is the fc(:) here. Note how it depends on [H2O]^2, but one of the [H2O] factors got entirely dropped in the MUSICA mechanism if you write it out as just SO3 + H2O -> H2SO4 + H2O (as in the .json, the "2" also got dropped from the reactants and needs to be put back)... we can't just multiply the constants, leaving out one of the H2Os is off by a factor of [H2O]... so about 17 orders of magnitude.

2. usr_CO_OH being stubbed out makes the mechanism basically unrealistic as written because CO + OH is the dominant global OH sink (and the only CO oxidation path in TS1, I think?) -- so all the tropospheric oxidants are meaningless until this is fixed.

Either a new MICM type has to be made or fall back to the pre-JPL19 two-term representation. Zero will really mess up the mechanism...

3. the het stubs are fine but this means no strat het chemistry and no ozone-hole chemistry.

More notes on SIMA forthcoming...

Ah, looks like the coefficient was dropped. Not sure how that happened, but fixing that bug will correct this rate

[jimmielin]

Ah, looks like the coefficient was dropped. Not sure how that happened, but fixing that bug will correct this rate

A also has to be rescaled to 8.5e-41 * m^2. If you're taking m ~ 6.022e17 then A would be around 3e-5... Claude can probably work out the right coeff. based on the CAM code.

[K20shores]

Ah, looks like the coefficient was dropped. Not sure how that happened, but fixing that bug will correct this rate

A also has to be rescaled to 8.5e-41 * m^2. If you're taking m ~ 6.022e17 then A would be around 3e-5... Claude can probably work out the right coeff. based on the CAM code.

I don't think we need Claude for this. This conversion is done using avogadro's number and I'm pretty sure it's just a bug somewhere in mechanism configuration that shouldn't be too hard to find.

[K20shores]

Ah, looks like the coefficient was dropped. Not sure how that happened, but fixing that bug will correct this rate

A also has to be rescaled to 8.5e-41 * m^2. If you're taking m ~ 6.022e17 then A would be around 3e-5... Claude can probably work out the right coeff. based on the CAM code.

I don't think we need Claude for this. This conversion is done using avogadro's number and I'm pretty sure it's just a bug somewhere in mechanism configuration that shouldn't be too hard to find.

or I'm misunderstanding what you're saying entirely

[K20shores]

Thanks @K20shores for sending this over.

A few general chemistry notes that stood out to me then specifically for CAM-SIMA.

1. usr_SO3_H2O

musica/configs/v1/ts1/t1s2.json

Lines 3985 to 4013 in 10d504b

	{
	"type": "ARRHENIUS",
	"A": 5.1188196459999996e-23,
	"B": 0.0,
	"C": 6539.99676963515,
	"D": 300.0,
	"E": 0.0,
	"reactants": [
	{
	"name": "H2O",
	"coefficient": 1.0
	},
	{
	"name": "SO3",
	"coefficient": 1.0
	}
	],
	"products": [
	{
	"name": "H2SO4",
	"coefficient": 1.0
	},
	{
	"name": "H2O",
	"coefficient": 1.0
	}
	],
	"gas phase": "gas"
	},

I think the A here is mechanistically wrong. It appears that 5.1188e-23 was arrived here by 8.5e-41 * 6.022e17 from CAM

!-----------------------------------------------------------------
!       ... so3 + 2*h2o --> h2so4 + h2o
!       Note: this reaction proceeds by the 2 intermediate steps below
!           so3 + h2o --> adduct
!           adduct + h2o --> h2so4 + h2o
!               (Lovejoy et al., JCP, pp. 19911-19916, 1996)
!	The first order rate constant used here is recommended by JPL 2011.
!	This rate involves the water vapor number density.
!-----------------------------------------------------------------

       if ( usr_SO3_H2O_ndx > 0 ) then
          call comp_exp( exp_fac, 6540.0_r8*tinv(:), ncol )
          if( h2o_ndx > 0 ) then
             fc(:) = 8.5e-21_r8 * m(:,k) * h2ovmr(:,k) * exp_fac(:)
          else
             fc(:) = 8.5e-21_r8 * invariants(:,k,inv_h2o_ndx) * exp_fac(:)
          end if
          rxt(:,k,usr_SO3_H2O_ndx) = 1.0e-20_r8 * fc(:)
       end if

where 6.022e17 is the conversion factor for cm3 molec-1 s-1 to (mol m-3)-1 s-1 ("m") for a two-body rxn constant but the actual mechanics of this reaction is a fast pre-equilibrium

SO3 + H2O <-> SO3.H2O (adduct fast pre-equilibrium) SO3.H2O + H2O -> H2SO4 + H2O

the steady state observed rate is k(T) [SO3] [H2O]^2, where k(T) is the fc(:) here. Note how it depends on [H2O]^2, but one of the [H2O] factors got entirely dropped in the MUSICA mechanism if you write it out as just SO3 + H2O -> H2SO4 + H2O (as in the .json, the "2" also got dropped from the reactants and needs to be put back)... we can't just multiply the constants, leaving out one of the H2Os is off by a factor of [H2O]... so about 17 orders of magnitude.

2. usr_CO_OH being stubbed out makes the mechanism basically unrealistic as written because CO + OH is the dominant global OH sink (and the only CO oxidation path in TS1, I think?) -- so all the tropospheric oxidants are meaningless until this is fixed.

Either a new MICM type has to be made or fall back to the pre-JPL19 two-term representation. Zero will really mess up the mechanism...

3. the het stubs are fine but this means no strat het chemistry and no ozone-hole chemistry.

More notes on SIMA forthcoming...

Ah, I found it. Not a bug in musica, but the old chemistry cafe. It wrote out this reaction H2O + SO3 -> H2SO4 + H2O, so it was missing the stoichiometric coefficient of 2. When I correct it to be 2H2O + SO3 -> H2SO4 + H2O and run our conversion tool again, we get this

{
            "type": "ARRHENIUS",
            "A": 3.082625243326537e-05,
            "B": 0.0,
            "C": 6539.99676963515,
            "D": 300.0,
            "E": 0.0,
            "reactants": [
                {
                    "name": "H2O",
                    "coefficient": 2.0
                },
                {
                    "name": "SO3",
                    "coefficient": 1.0
                }
            ],
            "products": [
                {
                    "name": "H2SO4",
                    "coefficient": 1.0
                },
                {
                    "name": "H2O",
                    "coefficient": 1.0
                }
            ],
            "gas phase": "gas"
        }

which seems correct, i think

[jimmielin]

When I correct it to be 2H2O + SO3 -> H2SO4 + H2O and run our conversion tool again,

Ah, neat, since it gets to "A": 3.082625243326537e-05, the same as my back of envelope calculation I think that's correct. I didn't know you generated the coeffs programmatically too, that's really neat. If I did it by hand I would probably asked Claude to confirm :-)

[K20shores]

When I correct it to be 2H2O + SO3 -> H2SO4 + H2O and run our conversion tool again,

Ah, neat, since it gets to "A": 3.082625243326537e-05, the same as my back of envelope calculation I think that's correct. I didn't know you generated the coeffs programmatically too, that's really neat. If I did it by hand I would probably asked Claude to confirm :-)

Our conversion tool reads what Chemistry Cafe outputs with our mechanism configuration library. It converts from the old camp format to the newer format that uses SI units everywhere. It applies those conversions automatically. Louisa creates this information in Chemistry Cafe. This is likely a case where it's correct in Chemistry Cafe for CAM and potentially incorrect for the other output targets of Chemistry Cafe (kpp, box mox, music box, something else) because the main target is CAM, so it relies on Fortran code being correct in CAM rather than relying entirely on data, as we do for the mechanism in musica.