Issue with units of O3 loss in tagged simulation

[ncolombi]

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?

[yantosca]

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]

Thanks Bob!

[xpye98]

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]

Question for you @xpye98 @ncolombi: In GeosCore/fullchem_mod.F90 the loss is computed from KPP as:

       !====================================================================
       ! HISTORY (aka netCDF diagnostics)
       !
       ! Prod and loss of families or species [molec/cm3/s]
       !
       ! NOTE: KppId is the KPP ID # for each of the prod and loss
       ! diagnostic species.  This is the value used to index the
       ! KPP "C" array (in module gckpp_Global.F90).
       !====================================================================

       ! Chemical loss of species or families [molec/cm3/s]
       IF ( State_Diag%Archive_Loss ) THEN
          DO S = 1, State_Diag%Map_Loss%nSlots
             KppId = State_Diag%Map_Loss%slot2Id(S)
             State_Diag%Loss(I,J,L,S) = C(KppID) / DT
          ENDDO
       ENDIF

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 the LOx 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]

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:

   !====================================================================
   ! HISTORY (aka netCDF diagnostics)
   !
   ! Prod and loss of families or species [molec/cm3/s]
   !
   ! NOTE: KppId is the KPP ID # for each of the prod and loss
   ! diagnostic species.  This is the value used to index the
   ! KPP "VAR" array (in module gckpp_Global.F90).
   !====================================================================

   ! Chemical loss of species or families [molec/cm3/s]
   ! Change unit to [1/m3/s] for TagO3 simulation
   ! xpye, 2022/02/12
   IF ( State_Diag%Archive_Loss ) THEN
      id_O3 = Ind_('O3')
      id_NO2 = Ind_('NO2')
      id_NO3 = Ind_('NO3')
      id_PAN = Ind_('PAN')
      id_PPN = Ind_('PPN')
      id_MPAN = Ind_('MPAN')
      id_HNO4 = Ind_('HNO4')
      id_N2O5 = Ind_('N2O5')
      id_HNO3 = Ind_('HNO3')
      id_BrO = Ind_('BrO')
      id_HOBr = Ind_('HOBr')
      id_BrNO2 = Ind_('BrNO2')
      id_BrNO3 = Ind_('BrNO3')
      id_MPN = Ind_('MPN')
      id_ETHLN = Ind_('ETHLN')
      id_MVKN = Ind_('MVKN')
      id_MCRHN = Ind_('MCRHN')
      id_MCRHNB = Ind_('MCRHNB')
      id_PROPNN = Ind_('PROPNN')
      id_R4N2 = Ind_('R4N2')
      id_PRN1 = Ind_('PRN1')
      id_PRPN = Ind_('PRPN')
      id_R4N1 = Ind_('R4N1')
      id_HONIT = Ind_('HONIT')
      id_MONITS = Ind_('MONITS')
      id_MONITU = Ind_('MONITU')
      id_OLND = Ind_('OLND')
      id_OLNN = Ind_('OLNN')
      id_IHN1 = Ind_('IHN1')
      id_IHN2 = Ind_('IHN2')
      id_IHN3 = Ind_('IHN3')
      id_IHN4 = Ind_('IHN4')
      id_INPB = Ind_('INPB')
      id_INPD = Ind_('INPD')
      id_ICN = Ind_('ICN')
      id_IDN = Ind_('IDN')
      id_ITCN = Ind_('ITCN')
      id_ITHN = Ind_('ITHN')
      id_ISOPNOO1 = Ind_('ISOPNOO1')
      id_ISOPNOO2 = Ind_('ISOPNOO2')
      id_INO2B = Ind_('INO2B')
      id_INO2D = Ind_('INO2D')
      id_INA = Ind_('INA')
      id_IDHNBOO = Ind_('IDHNBOO')
      id_IDHNDOO1 = Ind_('IDHNDOO1')
      id_IDHNDOO2 = Ind_('IDHNDOO2')
      id_IHPNBOO = Ind_('IHPNBOO')
      id_IHPNDOO = Ind_('IHPNDOO')
      id_ICNOO = Ind_('ICNOO')
      id_IDNOO = Ind_('IDNOO')
      id_MACRNO2 = Ind_('MACRNO2')
      id_ClO = Ind_('ClO')
      id_HOCl = Ind_('HOCl')
      id_ClNO2 = Ind_('ClNO2')
      id_ClNO3 = Ind_('ClNO3')
      id_Cl2O2 = Ind_('Cl2O2')
      id_OClO = Ind_('OClO')
      id_O = Ind_('O')
      id_O1D = Ind_('O1D')
      id_IO = Ind_('IO')
      id_HOI = Ind_('HOI')
      id_IONO = Ind_('IONO')
      id_IONO2 = Ind_('IONO2')
      id_OIO = Ind_('OIO')
      id_I2O2 = Ind_('I2O2')
      id_I2O3 = Ind_('I2O3')
      id_I2O4 = Ind_('I2O4')

     ! Error check
      IF ( id_O3 <= 0 ) THEN
        PRINT *, "Error"
     ENDIF

     ! Get ozone molecular weight from species database
      O3_MW_g  = State_Chm%SpcData(id_O3)%Info%MW_g  ! g/mol

      DO S = 1, State_Diag%Map_Loss%nSlots
         KppId = State_Diag%Map_Loss%slot2Id(S)
         State_Diag%Loss(I,J,L,S) = VAR(KppID) / DT
         ! Calculate Ox mass, according to Ox family in KPP 
         Ox_Mass        = State_Chm%Species(I, J, L, id_O3     ) &  
                        + State_Chm%Species(I, J, L, id_NO2    ) &
                      + State_Chm%Species(I, J, L, id_NO3    )*2.0_fp &
                      + State_Chm%Species(I, J, L, id_PAN    ) &
                      + State_Chm%Species(I, J, L, id_PPN    ) &
                      + State_Chm%Species(I, J, L, id_MPAN    ) &
                      + State_Chm%Species(I, J, L, id_HNO4   ) &
                      + State_Chm%Species(I, J, L, id_N2O5   )*3.0_fp &
                      + State_Chm%Species(I, J, L, id_HNO3   ) &
                      + State_Chm%Species(I, J, L, id_BrO    ) &
                      + State_Chm%Species(I, J, L, id_HOBr   ) &
                      + State_Chm%Species(I, J, L, id_BrNO2  ) &
                      + State_Chm%Species(I, J, L, id_BrNO3  )*2.0_fp &
                      + State_Chm%Species(I, J, L, id_MPN    ) &
                      + State_Chm%Species(I, J, L, id_ETHLN  ) &
                      + State_Chm%Species(I, J, L, id_MVKN  ) &
                      + State_Chm%Species(I, J, L, id_MCRHN  ) &
                      + State_Chm%Species(I, J, L, id_MCRHNB  ) &
                      + State_Chm%Species(I, J, L, id_PROPNN  ) &
                      + State_Chm%Species(I, J, L, id_R4N2   ) &
                      + State_Chm%Species(I, J, L, id_PRN1   ) &
                      + State_Chm%Species(I, J, L, id_PRPN   ) &
                      + State_Chm%Species(I, J, L, id_R4N1   ) &
                      + State_Chm%Species(I, J, L, id_HONIT   ) &
                      + State_Chm%Species(I, J, L, id_MONITS   ) &
                      + State_Chm%Species(I, J, L, id_MONITU   ) &
                      + State_Chm%Species(I, J, L, id_OLND   ) &
                      + State_Chm%Species(I, J, L, id_OLNN   ) &
                      + State_Chm%Species(I, J, L, id_IHN1   ) &
                      + State_Chm%Species(I, J, L, id_IHN2   ) &
                      + State_Chm%Species(I, J, L, id_IHN3   ) &
                      + State_Chm%Species(I, J, L, id_IHN4   ) &
                      + State_Chm%Species(I, J, L, id_INPB   ) &
                      + State_Chm%Species(I, J, L, id_INPD   ) &
                      + State_Chm%Species(I, J, L, id_ICN  ) &
                      + State_Chm%Species(I, J, L, id_IDN  )*2.0_fp &
                      + State_Chm%Species(I, J, L, id_ITCN  ) &
                      + State_Chm%Species(I, J, L, id_ITHN  ) &
                      + State_Chm%Species(I, J, L, id_ISOPNOO1  ) &
                      + State_Chm%Species(I, J, L, id_ISOPNOO2  ) &
                      + State_Chm%Species(I, J, L, id_INO2B  ) &
                      + State_Chm%Species(I, J, L, id_INO2D  ) &
                      + State_Chm%Species(I, J, L, id_INA  ) &
                      + State_Chm%Species(I, J, L, id_IDHNBOO  ) &
                      + State_Chm%Species(I, J, L, id_IDHNDOO1  ) &
                      + State_Chm%Species(I, J, L, id_IDHNDOO2 ) &  
                      + State_Chm%Species(I, J, L, id_IHPNBOO ) &  
                      + State_Chm%Species(I, J, L, id_IHPNDOO ) &  
                      + State_Chm%Species(I, J, L, id_ICNOO   ) &  
                      + State_Chm%Species(I, J, L, id_IDNOO   ) *2.0_fp&  
                      + State_Chm%Species(I, J, L, id_MACRNO2 ) &  
                      + State_Chm%Species(I, J, L, id_ClO     ) &  
                      + State_Chm%Species(I, J, L, id_HOCl     ) &  
                      + State_Chm%Species(I, J, L, id_ClNO2     ) &  
                      + State_Chm%Species(I, J, L, id_ClNO3  )*2.0_fp&
                      + State_Chm%Species(I, J, L, id_Cl2O2  )*2.0_fp&
                      + State_Chm%Species(I, J, L, id_OClO   )*2.0_fp&
                      + State_Chm%Species(I, J, L, id_O      ) &
                      + State_Chm%Species(I, J, L, id_O1D    ) &
                      + State_Chm%Species(I, J, L, id_IO     ) &
                      + State_Chm%Species(I, J, L, id_HOI    ) &
                      + State_Chm%Species(I, J, L, id_IONO   ) &
                      + State_Chm%Species(I, J, L, id_IONO2  )*2.0_fp &
                      + State_Chm%Species(I, J, L, id_OIO    )*2.0_fp &
                      + State_Chm%Species(I, J, L, id_I2O2   )*2.0_fp &
                      + State_Chm%Species(I, J, L, id_I2O3   )*3.0_fp &
                      + State_Chm%Species(I, J, L, id_I2O4   )*4.0_fp 
         ! Grid box volume [cm3]
         BOXVL          = State_Met%AIRVOL(I, J, L) * 1e+6_fp
         ! Convert Ox_Mass from [molec/cm3] to [molec]
         Ox_Mass        = Ox_Mass*BOXVL
         ! Divide L(Ox) [molec/cm3/s] by Ox mass [molec] 
         ! and then convert to [1/m3/s] for HEMCO
         IF ( Ox_Mass > 0.0_fp ) THEN
                 State_Diag%Loss(I, J, L, S)  = &
           ( State_Diag%Loss(I, J, L, S) / Ox_Mass ) *  1.0e+6_fp
         ENDIF
      ENDDO
   ENDIF

   ! Chemical production of species or families [molec/cm3/s]
   ! Change unit to [kg/m3/s], for TagO3 simulation
   ! xpye, 2022/02/12
   IF ( State_Diag%Archive_Prod ) THEN
      DO S = 1, State_Diag%Map_Prod%nSlots
         KppID = State_Diag%Map_Prod%slot2Id(S)
         State_Diag%Prod(I,J,L,S) = VAR(KppID) / DT

         ! Convert P(Ox) from [molec/cm3/s] to [kg/cm3/s]
         State_Diag%Prod(I, J, L, S) = State_Diag%Prod(I, J, L, S) &
                                       / AVO*1.e-3_fp * O3_MW_g
         ! Convert to [kg/m3/s], for HEMCO 
         State_Diag%Prod(I, J, L, S)     = State_Diag%Prod(I, J, L, S)* 1e+6_fp
      ENDDO
   ENDIF

@xpye98, let me know if I misinterpreted anything here.

[yantosca]

@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 of Spc(id_O3)%Conc(I,J,L), etc. to simplify the code somewhat. Stay tuned.

[xpye98]

@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]

@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?

[yantosca]

@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?

I think it may just be due to historical code. In older versions you had to use bpch diagnostic output to archive the tagO3 but we removed that in favor of the netCDF History diagnostics. Also, the tagO3 simulation has probably been a bit stale (until you and others have recently started looking into it). Thanks for bringing this all up, we'll get it fixed!

[msulprizio]

@ncolombi @yantosca Is this issue resolved or should we reopen it? (The draft PR for a fix closed the issue automatically)

[ncolombi]

Hi Melissa, Apologies for the late response, but it is correct that this issue is not yet resolved-thanks for reopening and let me know if you need any clarifications about the issue! Thanks, Nadia

…

On Fri, Feb 9, 2024 at 2:54 PM Melissa Sulprizio ***@***.***> wrote: Reopened #2094 <#2094>. — Reply to this email directly, view it on GitHub <#2094 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ANIHFE2YKRA244RO4RV6JTTYS2LHRAVCNFSM6AAAAABBKSBE6GVHI2DSMVQWIX3LMV45UABCJFZXG5LFIV3GK3TUJZXXI2LGNFRWC5DJN5XDWMJRG43DENZYG4ZDCOI> . You are receiving this because you were mentioned.Message ID: ***@***.***>