Discrepancy of mass between HEMCO and SpeciesConc output in CH4 simulation

[ohsangik98]

Your name

Sang-Ik, Oh

Your affiliation

Seoul National University, South Korea

What happened? What did you expect to happen?

I appreciate your efforts for GEOS-Chem.

The sums of methane mass have a difference between the one calculated from the HEMCO_diagonstics file, and the one derived from the SpeciesConc output.

I ran a CH4 specialty simulation for a month (2010.01.01~2010.02.01) without any chemical decaying. Then, I calculated the total mass of methane which was emitted during a month in two ways.

1. HEMCO_diagnostics : emission rate [kg/m2/s] area [m2] time [s] = total 39.25 Tg
2. SpeciesConcMND : concentration [molec/cm3] 1e6 air volume [m3] / AVO * 16.04 = total 39.11 Tg

I archived the output of StatMet and SpeciesConc in an 'instantaneous' format. It is thought to be the same between two mass values.

I checked the log file with verbose mode and HEMCO/GC source files... but I couldn't find a potential reason, unfortunately.

What are the steps to reproduce the bug?

Run CH4 simulation with the options

• UseTotalPriorEmis : true
• CH4_LOSS_FREQ : false
• GLOBAL_OH : false
• GLOBAL_CL : false = No chemical sink, no soil absorption (which is not included in the prior emission file of my simulation)

Please attach any relevant configuration and log files.

HEMCO_Config.txt log.txt check_tag_sum_py.txt

I attach configuration, log, and simple Python code.

What GEOS-Chem version were you using?

14.3.1

What environment were you running GEOS-Chem on?

Local cluster

What compiler and version were you using?

ifort 15.0.0

Will you be addressing this bug yourself?

Yes

In what configuration were you running GEOS-Chem?

GCClassic

What simulation were you running?

CH4

As what resolution were you running GEOS-Chem?

2x2.5

What meterology fields did you use?

MERRA-2

Additional information

No response

[ohsangik98]

Additionally, there is a discrepancy in CH4 molar mass in global_ch4_mod.F90 and species_database.yml.

1. global_ch4_mod.F90 : XNUMOL_CH4 = AVO / 16d-3 this variable is used in calculating CH4 loss by oxidation.
2. species_database.yml : CH4_PROP, MW_g : 16.04 this variable seems to be used in converting units in hco_interface_gc_mod.F90, mixing_mod.F90.

This doesn't affect my issue, which does not include any loss process to examine mass conservation. But it may cause a little issue in calculating the CH4 budget.

Thank you.

[msulprizio]

Thanks for reporting this @ohsangik98. We are also looking into this on our end and will update here.

[nicholasbalasus]

Hi @ohsangik98 - I am able to replicate this lack of mass conservation in version 14.4.2. I ran 1 month simulations for January 2019 with all losses turned off (both the soil sink and 3D losses). In the first case, I also turn off emissions to show the minor noise from transport. In the other case, I leave emissions on and the mass balance gets progressively worse (linear plus some noise). Here, mass balance = mass(t=t) - mass(t=0) - sum(emissions(t=0:t)).

I'm not sure what is causing this. Maybe it is related to HEMCO precision (described here and here), but I'm not sure if we would expect a precision issue to accumulate linearly like this (@jimmielin). What do you think @msulprizio @yantosca @lizziel? Happy to take a look elsewhere or try different tests if you have suggestions.

In case this is an accounting issue, here is the code used for running GC and for calculating the mass balance.

GEOS-Chem

``` #!/bin/bash for emissions in yes no; do dir="${SCRATCH}/jacob_lab/${USER}/test-mass-conservation_emissions_${emissions}" rm -rf "${dir}" mkdir -p "${dir}" cd "${dir}" # Compile GCClassic for CH4 / MERRA-2 / 4 x 5 / 47 L git clone https://github.com/geoschem/GCClassic.git cd GCClassic git checkout 14.4.2 git submodule update --init --recursive source run/runScriptSamples/operational_examples/harvard_cannon/gcclassic.gcc10_cannon_rocky.env cd run gc_dir="gc" c="9\n1\n1\n2\n${dir}\n${gc_dir}\nn\n" printf ${c} | ./createRunDir.sh sed -i -e "s|17e-3_fp|17.01e-3_fp|g" \ -e "s|16d-3|16.04e-3_fp|g" "${dir}/GCClassic/src/GEOS-Chem/GeosCore/global_ch4_mod.F90" cd "${dir}/${gc_dir}/build" cmake ../CodeDir -DRUNDIR=.. make -j make install cd "${dir}/${gc_dir}" # Run for January 2019 and use BoundaryConditions to archive CH4 + dry air sed -i "/BoundaryConditions.fields: 'SpeciesBC_?ADV? ',/a\ 'Met_AD'," HISTORY.rc sed -i -e "s|CH4.frequency: 00000100 000000|CH4.frequency: 'End',|g" \ -e "s|CH4.duration: 00000100 000000|CH4.duration: 'End',|g" \ -e "s|#'BoundaryConditions',|'BoundaryConditions',|g" \ -e "s|BoundaryConditions.frequency: 00000000 030000|BoundaryConditions.frequency: 00000000 010000|g" \ -e "s|BoundaryConditions.duration: 00000001 000000|BoundaryConditions.duration: 00000000 010000|g" \ -e "s|'StateMet',|#'StateMet',|g" \ -e "s|'SpeciesConc',|#'SpeciesConc',|g" \ -e "s|'Metrics',|#'Metrics',|g" HISTORY.rc # Save emissions hourly sed -i -e "s|DiagnFreq: Monthly|DiagnFreq: Hourly|g" HEMCO_Config.rc # Turn off GFED since it varies daily sed -i -e "s|GFED : on CH4|GFED : off CH4|g" HEMCO_Config.rc # Optionally turn off emissions if [ $emissions = no ]; then sed -i -e "s|EMISSIONS : true|EMISSIONS : false|g" HEMCO_Config.rc fi # Turn off losses sed -i -e "s|2 OH_pert_factor 1.0 - - - xy 1 1|999 zero_scale_factor 0.0 - - - xy 1 1|g" \ -e "s|SpeciesConc_Cl 2010-2019/1-12/1/0 C xyz 1 \* - 1 1|SpeciesConc_Cl 2010-2019/1-12/1/0 C xyz 1 \* 999 1 1|g" \ -e "s|OH 1985/1-12/1/0 C xyz kg/m3 \* 2 1 1|OH 1985/1-12/1/0 C xyz kg/m3 \* 999 1 1|g" \ -e "s|CH4loss 1985/1-12/1/0 C xyz s-1 \* - 1 1|CH4loss 1985/1-12/1/0 C xyz s-1 \* 999 1 1|g" \ -e "s|CH4uptake 2009/1-12/1/0 C xy kg/m2/s CH4 1 14 2|CH4uptake 2009/1-12/1/0 C xy kg/m2/s CH4 999 14 2|g" HEMCO_Config.rc # Modify the run script then submit it cp runScriptSamples/operational_examples/harvard_cannon/geoschem.run . sed -i -e "s|-c 8|-c 32|g" \ -e "s|15000|96000|g" \ -e "s|sapphire,huce_cascade,seas_compute,shared|huce_ice,sapphire,huce_cascade,seas_compute|g" geoschem.run sbatch geoschem.run done ```

Mass Balance

``` import re import glob import numpy as np import pandas as pd import xarray as xr import matplotlib.pyplot as plt def get_methane_mass_and_emis(dir): files = sorted(glob.glob(dir + "GEOSChem.Boundary*.nc4")) hours = [re.search(r'(\d{8}_\d{4})', f).group(0) for f in files][0:-1] methane_mass = np.zeros((len(hours))) methane_emis = np.zeros((len(hours))) for idx,hour in enumerate(hours): state_file = dir + f"GEOSChem.BoundaryConditions.{hour}z.nc4" emis_file = dir + f"HEMCO_diagnostics.{hour.replace('_','')}.nc" with xr.open_dataset(state_file) as ds: mol_dry_air = ds["Met_AD"] * 1e3/28.9644 mol_methane = ds["SpeciesBC_CH4"] * mol_dry_air methane_mass[idx] = mol_methane.sum().values * 16.04/1e12 with xr.open_dataset(emis_file) as ds: for var in ds.variables: if ("Emis" not in var) | ("Total" in var): continue methane_emis[idx] += (ds[var] * 3600 * ds["AREA"]).sum().values/1e9 mass_balance = np.zeros((len(hours))) for idx in range(len(mass_balance)): mass_balance[idx] = methane_mass[idx] - methane_mass[0] - np.sum(methane_emis[0:idx]) return mass_balance no_emi = "/n/holyscratch01/jacob_lab/nbalasus/test-mass-conservation_emissions_no/gc/OutputDir/" mass_bal_no_emis = get_methane_mass_and_emis(no_emi) yes_emi = "/n/holyscratch01/jacob_lab/nbalasus/test-mass-conservation_emissions_yes/gc/OutputDir/" mass_bal_yes_emis = get_methane_mass_and_emis(yes_emi) fig,axs = plt.subplots(2,1,figsize=(7,6),sharex=True,layout="compressed") axs[0].plot(range(len(mass_bal_no_emis)), mass_bal_no_emis, linewidth=0.75) axs[1].plot(range(len(mass_bal_yes_emis)), mass_bal_yes_emis, linewidth=0.75) for ax in axs: ax.set_ylabel("Mass Balance [Tg]") axs[1].set_xlabel("Hour") axs[0].set_title("Emissions turned off", loc="left") axs[1].set_title("Emissions turned on", loc="left") fig.savefig("test.png", dpi=300, bbox_inches="tight") ```

[yantosca]

Thanks @nicholasbalasus for looking into this. We saw this during the implementation of flexible precision, which allows you to toggle between REAL*4 and REAL*8 at compile time. But in practice the REAL*4 precision leads to truncation that causes mass non-conservation.

• https://wiki.seas.harvard.edu/geos-chem/index.php/Flexible_precision_in_GEOS-Chem#Validation

One of our development items is to allow HEMCO to save output at either REAL*4 or REAL*8. Right now it can only save out REAL*4. This would likely remove the issue. @jimmielin @lizziel please feel free to comment further.

[ohsangik98]

Thanks for all the advice from @yantosca and the GEOS-Chem team.

I just checked mass discrepancies between the three options below :

1. REAL 4 / Non-local PBL : -0.355% (default, case above)
2. REAL 8 / Non-local PBL : -0.263%
3. REAL 8 / Full-mixed PBL : -0.00005%, which shows mass conservation :) ( Values above are relative differences between HEMCO-mass and SpeciesConc-mass )

Thanks again for solving the issue.

[nicholasbalasus]

Agreed! Thanks for finding this @ohsangik98!

[yantosca]

Hi all. It seems that moving the HEMCO-specific code out of the vdiff_mod.F90 may be contributing to this error. Especially when emissions are turned off.

@msulprizio wrote:

Jintai Lin's original fix (in VDIFF, to correct for mass conservation) was fairly straightforward and can be seen here.

One thing that stands out to me is the description of the variable sum_qp0 added by Jintai which says “total mass in the PBL (ignoring the v/v -> m/m conversion) including pre-mixing mass and surface flux (emis+drydep)”. I reviewed all changes to vdiff_mod.F90 since then and found that the HEMCO code (i.e. emissions/drydep) was removed in #311. This is part of an ongoing issue I believe (see #327). I wonder if this is contributing to the reintroduction of the mass conservation issue?

For now, we would recommend using TURBDAY for the CH4, CO2, or carbon simulations until we can provide a fix.

A quick fix might be to activate the HEMCO-specific code on each timestep but bracket out the HEMCO-specific part if emissions are turned off. Then that should allow for the drydep to be applied.

[yantosca]

@ohsangik98 @nicholasbalasus @msulprizio: I've printed out the values of State_Chm%SurfaceFlx (which is computed in the routine Surface_Flux_for_VDIFF as emissions flux - drydep flux. These do not change with time for the CH4 simulation:

444: sflx TRACER            1 :    7.3325212413926914E-008
451: sflx TRACER            1 :    7.3325212413926914E-008
458: sflx TRACER            1 :    7.3325212413926914E-008
465: sflx TRACER            1 :    7.3325212413926914E-008
472: sflx TRACER            1 :    7.3325212413926914E-008
479: sflx TRACER            1 :    7.3325212413926914E-008
487: sflx TRACER            1 :    7.3325212413926914E-008
494: sflx TRACER            1 :    7.3325212413926914E-008
501: sflx TRACER            1 :    7.3325212413926914E-008
508: sflx TRACER            1 :    7.3325212413926914E-008
515: sflx TRACER            1 :    7.3325212413926914E-008
522: sflx TRACER            1 :    7.3325212413926914E-008
530: sflx TRACER            1 :    7.3325212413926914E-008
537: sflx TRACER            1 :    7.3325212413926914E-008
544: sflx TRACER            1 :    7.3325212413926914E-008
551: sflx TRACER            1 :    7.3325212413926914E-008
558: sflx TRACER            1 :    7.3325212413926914E-008
565: sflx TRACER            1 :    7.3325212413926914E-008
579: sflx TRACER            1 :    7.3325212413926914E-008
586: sflx TRACER            1 :    7.3325212413926914E-008
593: sflx TRACER            1 :    7.3325212413926914E-008
600: sflx TRACER            1 :    7.3325212413926914E-008
607: sflx TRACER            1 :    7.3325212413926914E-008
614: sflx TRACER            1 :    7.3325212413926914E-008
622: sflx TRACER            1 :    7.3325212413926914E-008
...

So this would lead me to believe that the non-conservation issue is in vdiff_mod.F90

[yantosca]

I was able to check out and build the version of the code in which Jintai's updates went in (13.0.0-alpha.11). I had to manually comment out the places where we obtained the pointers to HEMCO in global_ch4_mod.F90, Here we see that mass is conserved:

GC.log:972: ### CH4 after:   0.23785176364278848     
GC.log:978: ### CH4 after:   0.23785206996431299     
GC.log:984: ### CH4 after:   0.23785237635125683     
GC.log:990: ### CH4 after:   0.23785268278374674     
GC.log:996: ### CH4 after:   0.23785298925073003     
GC.log:1002: ### CH4 after:   0.23785329574364100     
GC.log:1009: ### CH4 after:   0.23785360244275555     
GC.log:1015: ### CH4 after:   0.23785390892200831     
GC.log:1021: ### CH4 after:   0.23785421537156171     
GC.log:1027: ### CH4 after:   0.23785452180508559     
GC.log:1033: ### CH4 after:   0.23785482822573739     
GC.log:1039: ### CH4 after:   0.23785513463368330     
GC.log:1046: ### CH4 after:   0.23785544268393111     
GC.log:1052: ### CH4 after:   0.23785574982498067     
GC.log:1058: ### CH4 after:   0.23785605663713910     
GC.log:1064: ### CH4 after:   0.23785636328042969     
GC.log:1070: ### CH4 after:   0.23785666982252773     
GC.log:1076: ### CH4 after:   0.23785697629723396     
GC.log:1089: ### CH4 after:   0.23785728537366271     
GC.log:1095: ### CH4 after:   0.23785759530267203     
GC.log:1101: ### CH4 after:   0.23785790312041294     
GC.log:1107: ### CH4 after:   0.23785821073887994     
GC.log:1113: ### CH4 after:   0.23785851823728288
... 

Also, the surface fluxes (emissions - drydep) that are passed into VDIFF are constant, as they are in the current code:

971: sflx TRACER            1 :    7.5654374078286164E-008
977: sflx TRACER            1 :    7.5654374078286164E-008
983: sflx TRACER            1 :    7.5654374078286164E-008
989: sflx TRACER            1 :    7.5654374078286164E-008
995: sflx TRACER            1 :    7.5654374078286164E-008
1001: sflx TRACER            1 :    7.5654374078286164E-008
1008: sflx TRACER            1 :    7.5654374078286164E-008
1014: sflx TRACER            1 :    7.5654374078286164E-008
1020: sflx TRACER            1 :    7.5654374078286164E-008
1026: sflx TRACER            1 :    7.5654374078286164E-008
1032: sflx TRACER            1 :    7.5654374078286164E-008
1038: sflx TRACER            1 :    7.5654374078286164E-008
1045: sflx TRACER            1 :    7.5654374078286164E-008
1051: sflx TRACER            1 :    7.5654374078286164E-008
1057: sflx TRACER            1 :    7.5654374078286164E-008
1063: sflx TRACER            1 :    7.5654374078286164E-008
... 

@msulprizio @nicholasbalasus: I am trying to find what is different in vdiff_mod.F90 between these 2 versions but if you have any insights, please let me know. I may try putting debug printout in both to see what happens.

Also note: emissions are different in the 2 versions but that shouldn't matter as much as we are looking for trends in mass.

[nicholasbalasus]

This looks like it still exists back to 13.0.0 at least.

[msulprizio]

Just a quick clarification that Jintai's original mass conservation fix for non-local PBL mixing went into 12.1.0, not 13.0.0. See the original wiki post archived here.

The removal of the HEMCO code from vdiff_mod.F90 did go into 13.0. Since @nicholasbalasus found that the conservation issue is present in that version that still supports the theory that those updates may a have reintroduced the discrepancies in mass.

[nicholasbalasus]

This unfortunately predates those changes and is present in at least 12.9.0. Thanks to @yantosca for knowing how to set up a 12.9.0 run directory :)

[ohsangik98]

Thanks for everyone's effort.

When I printed the scale factor calculated by Jintai's mass conservation fix part, I found that it has a smaller order than mass unbalance scale, leading me to believe that the Jintai part is not the cause of our issue.

(Sorry for not attaching any plot or file.. I'm not available right now)

[yantosca]

We can now close this issue as PR #2428 has been merged into the 14.5.0 development stream.