Hi Friedemann,
I've had a look at the code after hearing of the issue via other channels. Posting my reply here for maximum usefulness for the general community.
As far as I can see it is as you've suspected: the tracers in many modules of WRF-Chem (including GHG module) are converted to mole fractions using the
wet air density. This means that we get wrong values at the emission point an time due to local variability of water vapor. It seems that this correction is usually not large (I've found less then 0.5% discrepancy in mass after 6 hours of my tracer in my domain), but it should be corrected, since the effect will be larger for longer simulations with larger domains.
On the other hand the advection routines seem to work as they should for dry-air tracer (tracers in my run seem to be conserved right), but I was unable to confirm that explicitly from the code.
If I am correct then the errors are affecting only the initial concentration calculated at the emission point, but not the subsequent transport.
Would be very happy is someone from WRF-Chem crew can comment on that. More details below.
Details:
Emissions
Emissions in GHG module are calculated by add_ghg_fluxes subroutine (chem/module_ghg_fluxes.F). This routine uses rho_phy as air density (but it's not denoted in the code whether this is total density, i.e. including water vapour, or dry density). It is, however, very similar to how other modules in WRF-Chem treat the emissions, and to the mechanism mentioned in WRF-Chem FAQ (https://ruc.noaa.gov/wrf/wrf-chem/FAQ.htm). There it says explicitly that it's dry air mole fraction (I've marked critical moments in
bold):
Q: How are the gas phase emissions units converted from mole per km^2 per hour to the change in parts per million used by WRF-Chem?
A: It is assumed that you have examined the anthropogenic emissions web page that discusses how one can build the emissions using the emiss_v03.F program.
The units of all gas species in the binary files created with emiss_v03.F have units of mole/(km*km)/hr. The units of all aerosol species are microgram/(m*m)/s and they are treated differently than the gas phase emissions - hence the different emissions units. (These are also units of elevated (point) sources for each vertical grid).
In module_emissions_anthropogrnics.F WRF-Chem code, the subroutine add_anthropogenics has the variable "conv_rho" that is used to convert the emissions units (gas-phase chemistry) to ppm (note that the model works with them as a mixing ratio). Specifically, conv_rho is used in module_emissions_anthropogenics.F to convert the gas-phase emissions from mole/(km*km)/hr into delta-ppmv, where delta-ppmv is the incremental addition to the gas phase species (units of ppmv).
The calculation of conv_rho needs:
rho_phy: the air density in kg/(m3),
dtstep: the meteorology big time step in seconds,
dz8w: the vertical grid spacing (delta-z) for the lowest layer (or the layer for the emissions) in meters
Let E be emission in mole(/km-km)/hr from binary files (this is the output data from the emiss_v03.F program)
Then
E*(1.E-3)*(1.E-3)/(dz8w(I,K,J)*3600.) = E2 =emission in mole/(m^3)/sec
And then let D be the dry air density in mole/(m3), or
D=rho_phy(i,k,j)/(.02897 kg/mole)
And the mixing ratio in units of parts per million (ppmv) is equal to 1.E6 * mixing ratio where the mixing ratio units is (mole/mole).
So that leads to the conv_rho variable which converts the emissions units to parts per million being:
conv_rho = (1.E6)*dtstep*E2/D = E*dtstep*8.047E-6/(rho_phy(I,k,J)*dz8w(i,k,j))
So, the chemical emissions that are read into the WRF-Chem model and are multiplied by "conv_rho" to change the units.
However, in both subroutine add_anthropogenics (chem/module_emissions_anthropogenics.F) and add_ghg_fluxes (chem/module_ghg_fluxes.F),
rho_phy is passed from the parent subroutine emissions_driver (chem/emissions_driver.F), which in turns gets it as
rho from chem_driver. As you've already pointed out, chem driver calculates this rho in chem_prep (module_chem_utilities.F) using the formula:
rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV)),
where:
1/alt = rho: this is from base state variables; alt = inverse density of dry air* [m3 / kg]
moist(i,k,j,P_QV) is the water vapor mixing ratio [kg / kg]
*Now, I'm only 90% sure that 1/alt = rho_dry, i.e. density of dry-air (dry density), but the formula above only makes sense when it is. Assuming this is the case, then in the emissions in chem modules we actually use total density, including gaseous water vapor, i.e. wet air density. This contradicts the FAQ quoted above, and is not what I would expect to happen.
To me this means that we are in fact using the wrong density when converting fluxes to mole fractions in the atmosphere. It seems that this correction is usually not large (I've found less then 0.5% discrepancy in mass after 6 hours of continuous emissions of a point tracer in my domain), but it should be corrected, since the effect will be larger for longer simulations with larger domains.
Advection
Now the imporant question is whether the model handles the tracers as dry-air mole fractions when advecting. My experiments show that it actually does (e.g. there is no variation in the constant-value-LBC tracer over the domain). Which then means that the error only concerns the emission point (in time and space), after which the mole fractions should be conserved ok. I cannot confirm that with 100% confidence, unfortunately. I'm hoping that WRF-Chem developers can say something in this regard.
Potential solution:
I see this as a minor issue, but still I'd fix it in my working code... On the global level I'll leave the decision to the developers. chem_prep subroutine may need a bugfix, but this may have a domino effect on other chemistry modules. Another option is to leave chem_prep and make sure that chem modules tackle this on their own, but since different code owners do things differently (e.g. see below), this might be tricky.
Side note:
FYI: I found is that in subroutine add_anthropogenics (chem/module_emissions_anthropogenics.F) some compounds are apparently emitted inconsistently, using wet dry density and some using dry-air density. Specifically: https://github.com/wrf-model/WRF/blob/release-v4.0.1/chem/module_emissions_anthropogenics.F, lines 87-111, note conv_rho and conv_rho_aer usage).
Edits: minor language improvements