[BUG/ISSUE] Hg simulation in version 13.4 missing land and ocean emissions

[viral211]

What institution are you from?

Harvard

Description of the problem

I am testing the Hg simulation in version 13.4 and I think this version misses the ocean and land Hg emissions. In previous versions, these were calculated in routine EMISSMERCURY in mercury_mod. This routine is no longer there and I don't see it moved to HEMCO either. @yantosca, do you know where this calculation is done now? Thanks.

tagging @arifein @CharikleiaGRN @jennyfisher

[CharikleiaGRN]

Ciao! I would also like to mention that in the new restart file in the http://geoschemdata.wustl.edu/ for GC-Hg 13.4 the variables SNOW_HG_OCEAN, SNOW_HG_OCEAN, SNOW_HG_LAND, SNOW_HG_LAND_STORED are missing and there are some differences in the concentrations of species comparing to the restart file I had before for Viral's Hg simulation.

[CharikleiaGRN]

One more comment :) I checked the total deposition of the earlier version of @viral211 's new chemistry mechanism vs v13.4.

Total Deposition Before: ~6900 T Total Deposition v13.4: ~3000 T

Could this difference be related to low deposition from missing ocean and land emissions?

[yantosca]

Hi all, thanks for writing. I have been out sick since IGC10 and am just back now. Thanks for your patience.

I believe the land and ocean codes were commented out (maybe in the code I got from @msl3v) since we were only looking into the gas-phase chemistry. Also, I wasn't sure how the species in the land/ocean mercury modules mapped into the new Hg species. Can someone clarify?

[arifein]

Hi @yantosca, hope you're feeling better now! As far as I know, all of these re-emission processes will only emit elemental mercury (Hg0). In the EMISSMERCURY subroutine in mercury_mod, we have:

    !--------------------------------------------------------------
    ! Here we are NOT using the Global Terrstrial Mercury Model...
    !--------------------------------------------------------------
    CALL LAND_MERCURY_FLUX ( EHg0_ln, LHGSNOW, State_Grid, State_Met )
    IF ( prtDebug )  CALL DEBUG_MSG( '### EMISSMERCURY: a LAND_FLUX' )

    CALL VEGEMIS( Input_Opt, State_Met, LVEGEMIS,  EHg0_dist, EHg0_vg,   RC )
    IF ( prtDebug ) CALL DEBUG_MSG( '### EMISSMERCURY: a VEGEMIS' )

    CALL SOILEMIS( EHg0_dist, EHg0_so, State_Grid, State_Met )
    IF ( prtDebug ) CALL DEBUG_MSG( '### EMISSMERCURY: a SOILEMIS' )

    !==========================================================================
    !ENDIF
    !==========================================================================

    CALL SNOWPACK_MERCURY_FLUX( EHg0_snow,  LHGSNOW,  State_Chm, &
                                State_Grid, State_Met )
    IF ( prtDebug ) CALL DEBUG_MSG( '### EMISSMERCURY: a SNOW_FLUX' )

    CALL BIOMASSHG( Input_Opt, EHg0_bb, RC )
    IF ( prtDebug ) CALL DEBUG_MSG( '### EMISSMERCURY: a BIOMASS' )

VEGEMIS (vegetation emissions) has been set to 0 for the last versions of GEOS-Chem that I've seen, so I'm not sure if it's still needed? And BIOMASSHG would be covered already by the HEMCO emissions that you added. So the crucial subroutines are LAND_MERCURY_FLUX, SOILEMIS, SNOWPACK_MERCURY_FLUX, and OFFLINEOCEAN_READMO (the offline ocean Hg0 exchange routine). Both OFFLINEOCEAN_READMO and SOILEMIS only include Hg0. The re-emission routines ( LAND_MERCURY_FLUX and SNOWPACK_MERCURY_FLUX) also only emit to Hg0, but they rely on counting Hg2+ species that have deposited. This is covered by routines in depo_mercury_mod.F90, and as far as I can see it should work with the new Hg2+ species as long as they have the attribute Is_Hg2=T and the new HgP species Is_HgP=T in the species database. I see however that these are not yet in the new species_database.yml file. Let me know if I can clarify something further, it seems from other issues #1284 #1286 that other edits to the species database file will be needed at the same time.

[viral211]

Thanks for pointing this out @arifein. I have uploaded a new species_database file that include the Is_Hg2 and Is_HgP attributes in #1286.

[viral211]

@CharikleiaGRN The lower deposition is indeed due to lower emissions. Also, it would be best to spin up the model for a minimum of two years to minimize the influence of initial conditions. The SNOWHG* fields should be part of the restart files generated after the spin up.

[CharikleiaGRN]

Sure, always 2 years spin-up. I was just confused with the 13.4. I thought it was only me getting low conc. and dep. :S Thank you all!

[yantosca]

Hi all... question for you. In routine SOILEMIS we have:

       IF ( IS_SNOWFREE_LAND ) THEN

          ! attenuate solar radiation based on function of leaf area index
          ! Jacob and Wofsy 1990 equations 8 & 9
          TAUZ = State_Met%MODISLAI(I,J) * 0.5e+0_fp

          ! For very low and below-horizon solar zenith angles, use
          ! same attenuation as for SZA=85 degrees
          SUNCOSVALUE = MAX( State_Met%SUNCOS(I,J), 0.09e+0_fp )

          ! fraction of light reaching the surface is
          ! attenuated based on LAI
          LIGHTFRAC = EXP( -TAUZ / SUNCOSVALUE )

          ! Dry soil Hg concentration, ng Hg /g soil
          DRYSOIL_HG = DRYSOIL_PREIND_HG * EHg0_dist(I,J)

          ! Soil emissions, ng /m2 /h
          SOIL_EMIS = EXP( 0.0011 * State_Met%SWGDN(I,J) * LIGHTFRAC ) * &
                      DRYSOIL_HG * SOIL_EMIS_FAC

          ! Grid box surface area [m2]
          AREA_M2   = State_Grid%Area_M2(I,J)

          ! convert soilnat from ng /m2 /h -> kg /gridbox /s
          EHg0_so(I,J) = SOIL_EMIS * AREA_M2 * 1e-12_fp / &
                         ( 60e+0_fp * 60e+0_fp )

          ! Multiply by fractional land area
          EHg0_so(I,J) = EHg0_so(I,J) * FRAC_SNOWFREE_LAND

where EHg0_dist is the output from routine VEGEMIS. If you are saying that VEGEMIS has been turned off then this would also imply that the soil emissions are zero (by this equation):

SOIL_EMIS = EXP( 0.0011 * State_Met%SWGDN(I,J) * LIGHTFRAC ) * &
                      DRYSOIL_HG * SOIL_EMIS_FAC

[yantosca]

Never mind, I see that EHg0_dist is read from disk. Sorry.

[yantosca]

@viral211 @arifein: The species database update from commit https://github.com/geoschem/geos-chem/commit/3c676ec3ed78f55ad7a9135875056846e46c5052 has Is_Hg0 but not Is_Hg2 and Is_HgP tags. Could you let me know which species should be tagged with Is_Hg2 and Is_HgP? Sorry for the bother.

[viral211]

@yantosca: Here is a species database file that includes tags for the Hg2 and HgP species. species_database.txt

[yantosca]

Thanks @viral211. Now I've run into another issue. Because we had set State_Chm%N_Hg_CATS = 1 to remove the tagged species, I get this error:

ADVECTED SPECIES MENU
------------------------------------------------
  #  Species Name
  1  Hg0                            
  2  HgBr                           
  3  HgBrNO2                        
  4  HgBrHO2                        
  5  HgBrBrO                        
  6  HgBrClO                        
  7  HgBrOH                         
  8  HgBr2                          
  9  HgCl                           
 10  HgClNO2                        
 11  HgClHO2                        
 12  HgClClO                        
 13  HgClBrO                        
 14  HgClBr                         
 15  HgClOH                         
 16  HgOH                           
 17  HgOHNO2                        
 18  HgOHHO2                        
 19  HgOHClO                        
 20  HgOHBrO                        
 21  HgOHOH                         
 22  HgCl2                          
 23  Hg2ClP                         
 24  Hg2ORGP                        
 25  Hg2STRP                        
At line 2547 of file /local/ryantosca/GC/rundirs/epa-kpp/bugfix_Hg/CodeDir/src/GEOS-Chem/Headers/state_chm_mod.F90
Fortran runtime error: Index '1' of dimension 1 of array 'state_chm%hg2_id_list' above upper bound of 0

Error termination. Backtrace:
#0  0xd34731 in init_hg_simulation_fields
        at /local/ryantosca/GC/rundirs/epa-kpp/bugfix_Hg/CodeDir/src/GEOS-Chem/Headers/state_chm_mod.F90:2547

I think this is because the loops in the code over N_Hg_Cats assumed that each of Hg0, Hg2, HgP always had the same number of species. But it loosk like we now have 1 Hg0 species, 1 Hg2 species, and 3 HgP species. How did you work around this??

[viral211]

@yantosca, Hg0_id_list, Hg2_id_list and HgP_id_list were used only for the tagged Hg simulation. I had commented them out everywhere they appeared.

[yantosca]

Thanks @viral211. I am going to get rid of the Hg categories variables everywhere. Then it should work. My code still had those defined.

[yantosca]

OK, now I get:

---> DATE: 2019/01/01  UTC: 00:00  X-HRS:      0.000000
 HEMCO already called for this timestep. Returning.
===============================================================================
TPCORE_FVDAS (based on GMI) Tracer Transport Module successfully initialized
===============================================================================
 ### min, max ocean_conc:    0.00000000       1.07819281E-10
 ### min, max EHg0_ln      0.0000000000000000        0.0000000000000000     
 ### min, max EHg0_oc     -3.9545772816530695E-006   6.2899588501098109E-004
 ### min, max EHg0_snow    0.0000000000000000        0.0000000000000000     
 ### min, max EHg0_so      0.0000000000000000        1.1551907062530933E-004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% USING O3 COLUMNS FROM THE MET FIELDS! %%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
---> DATE: 2019/01/01  UTC: 00:10  X-HRS:      0.166667
---> DATE: 2019/01/01  UTC: 00:20  X-HRS:      0.333333
 ### min, max ocean_conc:    0.00000000       1.07819281E-10
 ### min, max EHg0_ln      0.0000000000000000        8.2085032408869016E-004
 ### min, max EHg0_oc     -5.9728892192051022E-006   6.2899580178582602E-004
 ### min, max EHg0_snow    0.0000000000000000        1.5759451798084504E-007
 ### min, max EHg0_so      0.0000000000000000        1.1551907062530933E-004
---> DATE: 2019/01/01  UTC: 00:30  X-HRS:      0.500000
---> DATE: 2019/01/01  UTC: 00:40  X-HRS:      0.666667
 ### min, max ocean_conc:    0.00000000       1.07819281E-10
 ### min, max EHg0_ln      0.0000000000000000        1.3244678373512488E-004
 ### min, max EHg0_oc     -4.8335696131025117E-006   6.2899579483727155E-004
 ### min, max EHg0_snow    0.0000000000000000        1.7375748068806004E-007
 ### min, max EHg0_so      0.0000000000000000        1.1551907062530933E-004
---> DATE: 2019/01/01  UTC: 00:50  X-HRS:      0.833333
     - Creating file for SpeciesConc; reference = 20190101 000000
        with filename = OutputDir/GEOSChem.SpeciesConc.20190101_0000z.nc4
     - Creating file for Restart; reference = 20190101 010000
        with filename = ./GEOSChem.Restart.20190101_0100z.nc4
---> DATE: 2019/01/01  UTC: 01:00  X-HRS:      1.000000
     - CLEANUP: deallocating arrays now...

so it looks like the Hg0 arrays that would be populated from routine EMISSMERCURY are getting updated. I'll need you all to take a look at this implementation once I push the updates.

[yantosca]

I've pushed my updates (commits 4c8f6cfbc94f6b1a58f84492dd8d19ca7e8931ee, 2c4e9400ad4ecb35b14886b62185785b10e602be, 2792f8549492fbde31ba7bad36e326aca8e8fd77, and e50fc2f00c20749486be99c75c0c4aa6c123a434) to the bugfix/Hg-sim branch.

@arifein @CharikleiaGRN @viral211, please test further and let me know of any other problems, especially in how the Hg0 emissions are updated in the new EMISSMERCURY routine.

[arifein]

Hi @yantosca, thanks for sending those new commits. I'm planning to run the simulation for a year (or more) to evaluate the output. One issue I found so far is that in the subroutine EMISSMERCURY, there is a missing call to the subroutine SRCHg0 to actually add the Hg0 emissions to the species and then save it to HISTORY diagnostics. I saw that you already updated the SRCHg0 subroutine to remove the anthropogenic emissions, but it wasn't called anywhere yet. I send here the edited snippet in context:

    IF ( Input_Opt%LNLPBL ) HG_EMIS = 0.0e+0_fp

    ! Zero arrays for Hg deposition
    CALL RESET_HG_DEP_ARRAYS()
    IF ( prtDebug ) CALL DEBUG_MSG( '### EMISSMERCURY: ' // &
         'a RESET_HG_DEP_ARRAYS' )

+    ! Add Hg(0) source into State_Chm%Species [kg]
+    CALL SRCHg0( Input_Opt,  State_Chm, State_Diag, State_Grid, State_Met, RC )
+   IF ( prtDebug ) CALL DEBUG_MSG( '### EMISSMERCURY: a SRCHg0' )

    ! Convert species units back to original unit
    CALL Convert_Spc_Units( Input_Opt, State_Chm, State_Grid, State_Met, &
                            OrigUnit,  RC )

[arifein]

Hello, another question from my side about the emissions for HEMCO_Diagn.rc. The units for Hg2ClP are listed as kg/m1/a by default–> should this be changed to kg/m2/s for consistency?

As another note, in previous versions, we had added the particulate Hg emissions to the gas phase HgII tracer. I saw that currently the emissions are added directly to HgP (Hg2ClP). I'm not sure if there is any significant difference in this behaviour, but might be something to keep in mind if someone is looking into Hg2/HgP.

###############################################################################
#####   Hg2 emissions                                                     #####
###############################################################################
EmisHgCl2_Anthro   HgCl2    0    -1   -1   2  kg/m2/s  HgCl2_emission_flux_from_anthropogenic

###############################################################################
#####   Hg2 emissions                                                     #####
###############################################################################
EmisHg2ClP_Anthro  Hg2ClP   0    -1   -1   2  kg/m1/a  Hg2ClP_emission_flux_from_anthropogenic

[jennyfisher]

Hi @arifein -

I can't comment on the units (but I would guess you're right) so hopefully @yantosca or @viral211 will answer that one.

On the other code, that's from the HEMCO_Diagn file, so it just affects how things get written out to the diagnostics, not how they get read into the model. For that you'd need to look at the HEMCO_Config.rc. In the preliminary version I had from @viral211, all Hg2 and HgP emissions (regardless of inventory) were added to the HgCl2 species. e.g.:

0 STREETS_HG2 $ROOT/MERCURY/v2020-07/Streets/Streets2019_Hg.nc Hg2 2000-2015/1/1/0 C xy kg/m2/s HGCL2   - 1 1
0 STREETS_HGP $ROOT/MERCURY/v2021-08/Streets/Streets2019_Hg.nc HgP 2000-2015/1/1/0 C xy kg/m2/s HGCL2   - 1 1

Can you confirm whether this behaviour is changed in your HEMCO_Config.rc? If so it would be worth getting Viral's opinion on that (as I thought gas-particle partitioning was happening online post-emission). It might also be worth checking the output of those diagnostics -- do they both have values, or is the particulate one (EmisHg2ClP_Anthro) actually just filled with zeros?

-Jenny

[viral211]

@arifein, the units for the EmisHg2ClP_Anthro in HEMCO_Diagn.rc should indeed be kg/m2/s. Regarding the allocation of HgP emissions, I had assigned all HgP emissions to HgCl2(g) and allowed the Hg(II) gas/particle partitioning to happen online. But assigning HgP emissions to HgClP if ok too. It wouldn't make much of a difference since once emitted the Hg(II) gas to particle ratio would quickly reach the equilibrium values calculated online.

[arifein]

Hi @viral211 @jennyfisher, yes has changed in the new HEMCO_Config.rc. Good to know that it shouldn't have a big effect on results!

(((STREETS
0 STREETS_HG0 $ROOT/MERCURY/v2021-09/Streets/Streets2019_Hg.nc Hg0 2015/1/1/0 C xy kg/m2/s Hg0    - 1 1
0 STREETS_HG2 $ROOT/MERCURY/v2021-09/Streets/Streets2019_Hg.nc Hg2 2015/1/1/0 C xy *       HgCl2  - 1 1
0 STREETS_HGP $ROOT/MERCURY/v2021-09/Streets/Streets2019_Hg.nc HgP 2015/1/1/0 C xy *       Hg2ClP - 1 1
)))STREETS

[yantosca]

! Add Hg(0) source into State_Chm%Species [kg]
CALL SRCHg0( Input_Opt,  State_Chm, State_Diag, State_Grid, State_Met, RC )
IF ( prtDebug ) CALL DEBUG_MSG( '### EMISSMERCURY: a SRCHg0' )

D'oh! Thanks @arifein for pointing this out. This is now fixed in 10003b80cec646f5c6c47107c15ee2174d5dc95b.

[arifein]

Hi all, the first year of simulation is finished and the spinup so far looks like it's going in the right direction. Since the restart file seems to contain quite low Hg concentrations (especially over the ocean), the required spinup will be at least 2-3 years, especially for the Southern Hemisphere.

Another question that came up when reviewing the Hg0 emission diagnostics from the ocean. The variable for net Hg0 emissions (EmigHg0ocean) from the ocean in the MercuryEmis diagnostics collection does not seem to match the difference between the one-way up-and-down fluxes (FluxHg0fromOceanToAir – FluxHg0fromAirToOcean), which are diagnostics from the MercuryOcean collection. For the first year, I calculate global fluxes of:

FluxHg0fromOceanToAir = 2559.7 Mg/yr FluxHg0fromAirToOcean = 669.6 Mg/yr Net Hg0 emission flux = 1890.1 Mg/yr

However, the EmigHg0ocean output = 3780.2 Mg/yr, which is exactly double the difference of the one-way fluxes. Looking at a relative difference map, these two methods of calculating the net Hg0 flux differ by a factor of 2 in every ocean grid box, so it seems like a there is a unit error or missing factor of 2 somewhere. If I do the same for Viral's 12.8 code, I get the same value globally for (FluxHg0fromOceanToAir – FluxHg0fromAirToOcean) = EmigHg0ocean, so it seems like this difference was introduced somewhere in the implementation. I've been trying to find the issue in mercury_mod.F90, but I haven't found the culprit yet.

I think that the EmigHg0ocean variable is showing the correct flux and the one-way fluxes are too low by a factor of 2. If this is the case, this would just be an issue of output diagnostics rather than incorrect ocean emissions being introduced into the simulation.

[jennyfisher]

Hi Ari,

This is great, thanks so much!

I'm a bit concerned about the ocean fluxes! I'm not convinced that if the one-way fluxes are wrong this is just a problem for the diagnostics. FUP and FDOWN are getting saved in the diagnostic array directly in mercury_mod.F90, and so if there is a problem with these I would expect it could affect the simulation... This might be why the model is stabilising for the SH mid-latitudes at much lower concentrations than we've ever seen before (since the ocean source is so important here)...

I had a quick look through the code and I also don't see anything that looks off in mercury_mod.F90, but I can't figure out where else this would come in.

@yantosca - any idea what is going on here? Can you help dig into where each of these diagnostics gets their data?

Cheers, Jenny

[arifein]

Hi @jennyfisher, I was also worried that SH mid-latitudes were stabilizing at much lower concentrations, but it's not actually fully equilibrated after the 1st year. I ran a second year of the simulation: below, you can see that the Hg0 at 30 S is still increasing in the second year (but the seasonal cycle makes it look like it's equilibrated). I will continue running additional years to try to get a spun-up restart file, but it looks like the model is headed in the right direction. This is why I thought that the net emission flux is correct, but the FUP and FDOWN diagnostics are wrong. I will try to dig into this in the coming days with print statements from the model.

[arifein]

I just did some tests with writing out outputs (as printouts and netcdf) for EmigHg0ocean, FluxHg0fromOceanToAir, and FluxHg0fromAirToOcean. It seems like when I printout from the code or use instantaneous diagnostics, EmigHg0ocean = FluxHg0fromAirToOcean - FluxHg0fromOceanToAir, which is what we'd expect. Weirdly, if I use time-averaged netcdf diagnostics, EmigHg0ocean = 2 * (FluxHg0fromAirToOcean - FluxHg0fromOceanToAir). @yantosca have you seen this behaviour before?

I also found a couple of small typos in Headers/state_diag_mod.F90, but they don't have any effect on this behaviour: Line 8930:

CASE( 34 )
-                diagId = 'FluxHg0froimAirToOcean'
+                diagId = 'FluxHg0fromAirToOcean'

Line 11940:

     ELSE IF ( TRIM( Name_AllCaps ) == 'FLUXHG0FROMAIRTOOCEAN' ) THEN
        IF ( isDesc    ) Desc  = &
-            'Volatization flux of Hg0 from the ocean to the atmosphere'
+            'Volatization flux of Hg0 from the atmosphere to the ocean'
        IF ( isUnits   ) Units = 'kg s-1'
        IF ( isRank    ) Rank  =  2

     ELSE IF ( TRIM( Name_AllCaps ) == 'FLUXHG0FROMOCEANTOAIR' ) THEN
        IF ( isDesc    ) Desc  = &
-            'Deposition flux of Hg0 from the atmosphere to the ocean'
+            'Deposition flux of Hg0 from the ocean to the atmosphere'
        IF ( isUnits   ) Units = 'kg s-1'
        IF ( isRank    ) Rank  =  2

[yantosca]

@arifein, I'll take a look at this as soon as I can. I wonder if we need to zero out the diagnostic arrays on every call to mercury_mod.F90 so we don't carry around junk.

[jennyfisher]

Thanks Ari & Bob!

[arifein]

Another update from my side: I finished 3 years of spinup and use the 4th year for analysis. I send here a benchmark comparison between Shah et al. (2021) chemistry in GCv12.8 (Reference Model Version) and GCv14–with Bob's bug fixes (New Model Version). The distribution of TGM looks reasonable compared to observations, although v14 generally shows lower Hg concentrations in the Northern Midlatitudes. See the pdf for the full comparison.

I think the ocean emissions are working properly in the simulation (although the one way diagnostics are wrong), since the Southern Hemisphere Hg is similar in magnitude to v12.8. I wonder if the remaining differences are due to updates in the photolysis rates and oxidant concentrations read into HEMCO_Config.rc? For example, there are updates to the read-in files for:

(((JBrO
* JBrO           $ROOT/GCClassic_Output/13.0.0/$YYYY/GEOSChem.JValues.$YYYY$MM$DD_0000z.nc4          Jval_BrO          2010-2019/1-12/1/0 C xyz 1        * - 1 1
)))JBrO
(((LRED_JNO2
* JNO2           $ROOT/GCClassic_Output/13.0.0/$YYYY/GEOSChem.JValues.$YYYY$MM$DD_0000z.nc4          Jval_NO2          2010-2019/1-12/1/0 C xyz 1        * - 1 1
)))LRED_JNO2

vs. Viral's HEMCO_Config.rc:

* JBrO $ROOT/MERCURY/v2014-09/JVALUES/jBrO.daytime.geos5.2x25.nc JBrO 2002/1-12/1/0 C xyz s-1 * - 1 1
* JNO2 $ROOT/MERCURY/v2014-09/JVALUES/jvalues.noon.geos5.2x25.nc JNO2 2005/1-12/1/0 C xyz s-1 * - 1 1

Another question that I had, does anyone know why all oxidants are read from Oxidants_for_Hg_sim_2016.nc4, but then Br/BrO are also read from the GCv13 output? Is Br_GC or GLOBAL_Br ultimately used in the model?

# --- 2016 oxidant fields [molec/cm3] ---
# Consider supplying other years as necessary
* GLOBAL_Br    $ROOT/MERCURY/v2022-04/Oxidants_for_Hg_sim_2016.nc4  Br           2016/1-12/1/0 C xyz 1  * - 1 1
* GLOBAL_BrO   $ROOT/MERCURY/v2022-04/Oxidants_for_Hg_sim_2016.nc4  BrO          2016/1-12/1/0 C xyz 1  * - 1 1
* GLOBAL_CH4   $ROOT/MERCURY/v2022-04/Oxidants_for_Hg_sim_2016.nc4  CH4          2016/1-12/1/0 C xyz 1  * - 1 1
* GLOBAL_Cl    $ROOT/MERCURY/v2022-04/Oxidants_for_Hg_sim_2016.nc4  Cl           2016/1-12/1/0 C xyz 1  * - 1 1
* GLOBAL_ClO   $ROOT/MERCURY/v2022-04/Oxidants_for_Hg_sim_2016.nc4  ClO          2016/1-12/1/0 C xyz 1  * - 1 1
* GLOBAL_CO    $ROOT/MERCURY/v2022-04/Oxidants_for_Hg_sim_2016.nc4  CO           2016/1-12/1/0 C xyz 1  * - 1 1
* GLOBAL_HO2   $ROOT/MERCURY/v2022-04/Oxidants_for_Hg_sim_2016.nc4  HO2          2016/1-12/1/0 C xyz 1  * - 1 1
* GLOBAL_NO    $ROOT/MERCURY/v2022-04/Oxidants_for_Hg_sim_2016.nc4  NO           2016/1-12/1/0 C xyz 1  * - 1 1
* GLOBAL_NO2   $ROOT/MERCURY/v2022-04/Oxidants_for_Hg_sim_2016.nc4  NO2          2016/1-12/1/0 C xyz 1  * - 1 1
* GLOBAL_O3    $ROOT/MERCURY/v2022-04/Oxidants_for_Hg_sim_2016.nc4  O3           2016/1-12/1/0 C xyz 1  * - 1 1
* GLOBAL_OH    $ROOT/MERCURY/v2022-04/Oxidants_for_Hg_sim_2016.nc4  OH           2016/1-12/1/0 C xyz 1  * - 1 1

...

#==============================================================================
# --- BrOx fields --
#==============================================================================
(((BrOx_GC
* Br_GC          $ROOT/GCClassic_Output/13.0.0/$YYYY/GEOSChem.SpeciesConc.$YYYY$MM$DD_0000z.nc4      SpeciesConc_Br    2010-2019/1-12/1/0 C xyz 1        * - 1 1
* BrO_GC         $ROOT/GCClassic_Output/13.0.0/$YYYY/GEOSChem.SpeciesConc.$YYYY$MM$DD_0000z.nc4      SpeciesConc_BrO   2010-2019/1-12/1/0 C xyz 1        * - 1 1
)))BrOx_GC

benchmark_0102_0203_2.pdf

[jennyfisher]

This is great work Ari.

@viral211 - can you comment on Ari's questions above? I would think you are best placed to do so.

Looking at the PDF you linked to it seems like there are still maybe some problems with the ocean evasion - even the net looks a lot lower than in the reference version... Thoughts on that?

[arifein]

Hi Jenny, thanks for looking into that. The ocean evasion/net evasion issues are still the ones that I commented on before, which I think are related to diagnostics only. The net fluxes in the benchmark were plotted as FluxHg0fromAirToOcean - FluxHg0fromOceanToAir, which seem to be 2× too low everywhere. You can see here when I make the difference plot for EmisHg0ocean, the output looks better. So the one-way diagnostics fluxes still need to be fixed (seems like instantaneous outputs are accurate while time-averaged are wrong), but I don't think it affects model behaviour.

[viral211]

I wonder if the remaining differences are due to updates in the photolysis rates and oxidant concentrations read into > HEMCO_Config.rc?

The JBrO and JNO2 fields are no longer used. Neither are the Br_GC and BrO_GC fields. These lines were carried over from the previous version of HEMCO_Config.rc and are not needed anymore. @arifein's benchmark file shows lower Hg2 emissions in the new simulation. I am not sure why that is, but it is one reason for the lower Hg0 concentrations in the mid latitudes in the new simulation. It doesn't seem to be the only reason though.

[arifein]

@viral211 thanks for clarifying that. I didn't see any other differences so far in the HEMCO_Config.rc files.

One thing that I am still thinking about regarding the Hg2 emissions is the units question in HEMCO. One change since your 12.8 code version is the issue you brought up in #1286, that:

Change mol. wts. of all Hg(I) and Hg(II) species to that of Hg(0) (200.59). This is because the Hg(I) and Hg(II) masses are tracked on the basis of kgHg... In earlier model versions, we could use specify two separate values of MWs, one for unit conversion (emMW) and one for dry dep. and other processes (MW), but this is not supported anymore.

If you look at the difference for global EmisHg2HgPanthro between Ref (539.65 Mg) and New (399.02 Mg), the ratio (1.35) is the same as the ratio between the mol. wts. of HgCl2 (271.52) and Hg (200.59). 399 Mg is what I would expect from AMAP anthropogenic Hg(II) emissions from GC v12 runs, before Hg(II) species were added.

I guess I can check if this is the only difference between the two versions by artificially scaling Hg(II) emissions by 1.35 to see whether the benchmark is the same?

[jennyfisher]

@arifein that sounds like a useful test!

@yantosca any updates with trying to figure out what is going wrong with the ocean diagnostics?

Cheers, Jenny

[jennyfisher]

P.S. We should add Hannah Horowitz to this ticket (newest GCHg co-chair) -- does anyone know her Gitub handle?

[msulprizio]

P.S. We should add Hannah Horowitz to this ticket (newest GCHg co-chair) -- does anyone know her Gitub handle?

Tagging @hmhorow

[yantosca]

Hi all, have been swamped working on docs for the 14.0.0 release. Not sure when I can get back to this. Also I learned that we need to prioritize the carboncycle mechanism for another project.

[arifein]

Hi all, I wanted to give an update on the comparison between Hg in the working GEOS-Chem code version (v14-bugfix) and v12 with the chemistry from Shah et al. (2021). The synopsis is that the model doesn't look too bad compared to the observations, but there continue to be unexplained differences from Viral's paper. I attach here the comparison between the two model versions at the surface for Hg0, and the vertical profiles of Hg0 and HgII. The full benchmark is attached below.

Two features I want to highlight: 1) In v14 Hg0 is lower in tropics/mid-latitudes and higher in polar regions at the surface 2) Although Hg0 is lower in v14 throughout the vertical profile at the equator, HgII is similar and even higher for v14 at some points in the profile.

My current thought is that something in the redox chemistry is driving these differences. I tried a simulation where the HgII emissions were scaled up by 35%, and it did not have a big effect on these features. I calculated the overall oxidation rate coefficients for Hg0 in both simulations: for v14 it seems about 11% faster than v12. I'm still trying to look into what might be causing this difference, but haven't seen any obvious causes yet. I'm wondering whether to accept these differences and add in the new dry deposition code, or to keep trying to puzzle this out.

benchmark_0102_0205_2015-2014.pdf

[viral211]

Thank you @arifein. This is very instructive. The two versions indeed look very similar in the tropospheric burden and spatial patterns of Hg0 and HgII. The difference in the stratosphere is kind of large, but I guess it is because of insufficient spin-up time in the new run. I agree that the remaining differences in the troposphere seem to be driven by the higher (net) Hg0 oxidation in version 14. It's not clear to me either why there would be a difference, but there were a couple of bug fixes that the GCST had made while incorporating my version 12 code in version 14 (#1196). I don't think that it's worth spending more time on trying to figure out what is causing this difference. However, if there is a significant imbalance in the total Hg deposition and Hg emissions in the new version, we could turn down the net oxidation by increasing the Hg2OrgP reduction rate constant. We generally use that parameter to balance the Hg budget.

[yantosca]

We can close this issue as PR #1356 has now been merged into the 14.0.1 development branch.

[yantosca]

GEOS-Chem integration tests all passed, except for CO2:

Summary of test results:
------------------------------------------------------------------------------
Execution tests passed: 47
Execution tests failed: 1
Execution tests not yet completed: 0

The error in the CO2 simulation was tracked down to an incorrect entry in HEMCO_Config.rc for CEDS_CO2_SHP, and has now been fixed in #1356.