Issue with units of O3 loss in tagged simulation #2094
Description
Hello,
@xpye98 and I have previously reported issues with the tagged O3 run in GEOS-Chem version 13 and 14. The fix in #1109 corrected the issues involved with unit conversion from archived Prod/Loss of O3 (units of molec/cm3/s). P(O3) and L(O3) are now both converted from molec/cm3/s to units of kilograms (kg) in tagged_o3_mod.F90 and concentration of O3 is calculated in lines 397 and 461 as shown below:
! Add P(O3) [kg] to the total O3 species
Spc(id_O3)%Conc(I,J,L) = Spc(id_O3)%Conc(I,J,L) + PO3_kg
! Apply chemical loss [kg]
Spc(N)%Conc(I,J,L) = Spc(N)%Conc(I,J,L) - LO3_kg
In the second equation, where each term is in units of kilograms, as long as LO3_kg is greater than PO3_kg at a given timestep, the SpcConc will have a value of zero, which is what I see in my 14.2.0 simulation for all tagged species (zero everywhere). @xpye98 has corrected this issue in his v13.2.1 run by modifying the tagged_o3_mod.F90 file as such:
! Apply chemical loss [kg], LO3 is [1], Spc is [kg]
Spc(I,J,L,N) = Spc(I,J,L,N) - LO3 * Spc(I,J,L,N)
In this case, LO3 has been changed to units of [1] and multiplied by SpcConc [kg], which fixes the issue. Since the SpcConc has an initial value of 0, this means that P(O3) will always be greater than LO3*SpcConc, and SpcConc will slowly reach a steady state concentration after some spin up. Outputting LO3 in units of [1] rather than [kg] requires dividing L(O3) by Ox mass in the fullchem_mod.F90 file. Is there anyway we can incorporate this fix in the next release of version 14?
Activity
yantosca commentedon Jan 2, 2024
Thanks for the reminder @ncolombi! I'll work on getting this into the no-diff-to-benchmark branch, so that it can go into version 14.
ncolombi commentedon Jan 2, 2024
Thanks Bob!
xpye98 commentedon Jan 2, 2024
Thanks for tagging me. I believe previous GC (before 12.0.0 or so) used this kind of different units for Ox_prod and Ox_loss, which could produce the correct result. #973 also reported similar bug and fixed it by outputting different unit of Ox_loss.
yantosca commentedon Jan 2, 2024
Question for you @xpye98 @ncolombi: In
GeosCore/fullchem_mod.F90
the loss is computed from KPP as:for each of the loss families that are defined in
KPP/fullchem/fullchem.kpp
:https://github.com/geoschem/geos-chem/blob/3cd802eb4a9b7a226254fc3938778c8db8c07cf2/KPP/fullchem/fullchem.kpp#L13C1-L20C14
that is
LOx
,LCO
,LCH4
. These are in molec/cm3, so we divide by DT to get molec/cm3/s.Because
LOx
is a family species, is it appropriate to divide it by the concentration of O3? That would probably be approximate but am wondering if we'd have to also take into account the other species in theLOx
family into the computation.We could also add a
LO3
family to KPP, which would just be the loss of the species O3, and then divide by the concentration of O3 to get the loss rate in the units that you want. Again, I'm not sure if that would be scientifically appropriate for the purposes of the tagged O3 simulation.Also I'm thinking to add this to a different diagnostic, as the
Prod_LOx
might be used as input to other simulations, and thus if we change the units we would break those simulations.ncolombi commentedon Jan 2, 2024
Hi Bob, I believe you would need to divide LOx by the mass of the Ox family (not just the mass of O3). @xpye98 does that in his flexchem_mod.F90 (from version 13) in lines 1111-1289:
@xpye98, let me know if I misinterpreted anything here.
yantosca commentedon Jan 2, 2024
@xpye98 e98 @ncolombi: Thanks. I'll see how I can implement this. A few things have changed since the 13.1 code, especially the update to KPP 3.0.0. Also, I think we can also use
C(ind_O3)
instead ofSpc(id_O3)%Conc(I,J,L)
, etc. to simplify the code somewhat. Stay tuned.xpye98 commentedon Jan 2, 2024
@yantosca @ncolombi
Yes, in my practice (where I got the correct tag results), I divide it by the mass of Ox family, not just O3, but note that the definitions of Ox family are different for different GC versions.
I agree with you that changing the unit would affect other analysis, maybe outputting Ox_loss with both units (one for consistency with Ox_pord, and the other for tagged O3) would be good.
Outputting just PO3 and LO3 (not Ox family) to KPP won't work. Actually we've tried this before and got very strange results. The primary formation pathway of O3 is NO2 -> O3 + NO and the primary loss is O3 photolysis. Considering only O3 will get very high PO3 and LO3, and somehow will break the tagO3 simulation. Scientifically, it is also meaningful to consider Ox family to avoid fast reactions between O3 and NOx. But this is just my opinion and might need further test.
xpye98 commentedon Jan 2, 2024
@yantosca Also, I'm more curious of why at the very beginning we used different units for prod and loss in previous GC versions. Is this the only way we can tag correctly?
Loss_Ox
diagnostic from molec/cm3/s to unitless (to fix an issue in tagO3 simulations) #209617 remaining items