Scheduled Downtime
On Friday 21 April 2023 @ 5pm MT, this website will be down for maintenance and expected to return online the morning of 24 April 2023 at the latest

rho_phy and unit of tracers into WRF atmosphere


New member
I'm using WRF to compute transport of passive greenhouse gas tracers (chem_opt=17, emiss_opt=17). Those are usually expressed as dry air mole fractions, and I assumed WRF does the same. A student's question prompted me to look into this and now I'm not sure whether emissions are actually emitted as wet air mole fractions.
The lines in question are in module_ghg_fluxes.F, lines 77f (in WRF 4.1.1):

         conv_rho=8.0461e-6/rho_phy(i,k,j)*dtstep/dz8w(i,k,j)  ! 8.0461e-6=molar_mass(air)/3600, [g/mol/s]                                                                                
         chem(i,k,j,p_co2_ant)= chem(i,k,j,p_co2_ant) + conv_rho* emis_ant(i,k,j,p_e_co2)
The question boils down to whether rho_phy is the density of dry air or the density of wet air. I traced the definition of rho_phy via emissions_driver.F to subroutine chem_prep in module_chem_utilities.F:

rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))

I didn't track down 'alt' (the input to chem_prep is grid%alt), but this looks like density of wet air to me. Comments in e.g. module_mosaic_cloudchem.F and module_sorgam_vbs_aqchem.F suggest that, too. On the other hand, a comment in module_chem_plumerise_scalar.F says it's dry air density ('dn0' in subroutine plumerise). And if it's wet air, wouldn't the equation be inconsistent because 8.0461e-6 is computed with the molar mass of dry air...?
In emissions_driver, conversion factors for some species (called 'conv'; analogous to 'conv_rho' above) contain the factor 1/rho_phy just as conv_rho, while others contain the factor alt.

Can someone explain which density rho_phy actually is?

Best regards,
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.



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 ( 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)

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)),
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.


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:, lines 87-111, note conv_rho and conv_rho_aer usage).

Edits: minor language improvements