EMEP-WRF global run: 4% discrepancy between reduced N emission and deposition budgets

[yaoead96]

Hello,

I have an enquiry regarding the global mass budgets of reduced N (NH3+NH4). The 'Mass Budgets Summary' file outputted by the model indicates that the global NH3 emission in 2015 is 53.1 TgN, while the total deposition of reduced N (NH3+NH4) is 55.5 TgN, which results in a ~4.5% relative difference (i.e. (Deposition-Emission)/Emission ). My own calculation of the annual global reduced N budgets also has the same discrepancy. The details of the model simulation is listed below. I am wondering is there any potential explanation for unbalanced budget? By the way, may I know how does the global EMEP model deal with the BIC as well?
MassBudgetSummary_2015.txt

Model version: EMEP4.34 + WRF4.1.1
Domain: Global run excluding the North Pole belt with resolution of 1° x 1°.
Year: 2015
Emission inventory: ECLIPSE v6a

I have also enclosed the 'Mass Budgets Summary' file for your information. Please feel free to contact me should you need any further information.

Thank you very much for your help.

Yao (PhD student)
UKCEH & University of Edinburgh
UKCEH Supervisor: Dr Massimo Vieno

[yaoead96]

Hello,

This is for your update.

The model runs in my previous post are using EQSAM gas/particle portioning scheme. I have previously run the same global simulations but with MARS scheme. The model outputs with MARS show a positive difference between total RDN emission and deposition as well as a smaller relative difference (0.91%) in 2015. Detailed numbers can be found in attached 2 screenshots. It seems that it is the aerosol scheme that drives the budget difference.

[mifads]

Hi Yao,
Sorry for the slow reply. There seems to be a problem, yes, and we are looking at this. Note though that the mass budget should consist of more than emis and deposition. One needs to consider initial and final masses in the domain, and any fluxes across the side or top boundaries (not such an issue for global NHx simulations, but important for e.g. European ones).

[mvieno]

Hi @mifads thanks for the update :)

This budget was done for a global run, and the run log of the run shows that EMEP did recognised it as a global.

Does the fluxin and fluxout in the MassBalance not including any lateral or top boundary flux?

[gitpeterwind]

The budgets include fluxes from top and lateral boundaries if they exist.
The loss of N can be a sign that the main timesplitting timestep is too large. You can adjust it using config_emep.nml settings, for example:
dt_advec=600,

[mvieno]

@gitpeterwind ok, @yaoead96 can you try that. I still not clear how to explain the difference between the two aerosols equilibrium schemes.

[yaoead96]

Hi, this is for your update. Following my previous posts, I have also calculated the same global quantities from the model run without Forest Fire Emissions (FFE). The results are tabulated below. It seems that the exclusion of forest fire emissions has no impact on the discrepancies between global RDN emissions and deposition.

[mifads]

Not sure I can help with the code, but I still don't understand why you only have emission and deposition in your tables. The perfect mass budget should have something like

[ inputs ] - [outputs ] = 0.0
[initial_mass + emission + fluxin (sides+top)] - [  deposition + fluxout (sides+top) + final_mass ] = 0.0

We are not getting zero though, and for S we get a fracmass (=[outputs]/[inputs]) of ca. 1.02 for global runs, e.g.

Mass balance   1    Sulphur
++++++++++++++++++++++++++++++++++++++++++++++++
               sumini      summas     fluxout      fluxin    fracmass
 Sulphur   3.6409E+08  3.1595E+08  2.8777E+08  3.9019E+08  1.0229E+00

There is a problem, but less than 4-5% I hope. What do you get for this Sulphur family output (in RunLog or the standard output).

[mvieno]

@mifads see below the EMEP mass balance log for our GLOBAL run. The flux in for NH3 is zero and flux out is very small. Similarly for NH4, the flux in and out are relatively small compared to ddep and wdep. @yaoead96 does the 4% only apply to RDN (NHx)?

#n Spec usedMW emis ddep wdep init sum_mass fluxout fluxin frac_mass
90 NH3 17 6.17E+10 2.77E+10 1.04E+10 0.00E+00 1.30E+08 2.37E+05 0.00E+00 6.20E-01
93 NH4_f 18 0.00E+00 3.50E+11 1.11E+13 1.36E+13 2.23E+12 1.06E+08 2.50E+08 1.01E+00

++++++++++++++++++++++++++++++++++++++++++++++++
Mass balance 1 Sulphur
++++++++++++++++++++++++++++++++++++++++++++++++
sumini summas fluxout fluxin fracmass
Sulphur 4.5065E+13 8.6681E+12 1.5555E+11 7.1499E+10 1.0005E+00
ifam totddep totwdep totem
1 4.413E+13 2.092E+14 2.169E+14
++++++++++++++++++++++++++++++++++++++++++++++++
Mass balance 2 Nitrogen
++++++++++++++++++++++++++++++++++++++++++++++++
sumini summas fluxout fluxin fracmass
Nitrogen 5.5175E+13 1.6386E+13 6.0889E+11 8.8449E+11 1.0163E+00
ifam totddep totwdep totem
2 3.764E+12 3.634E+13 1.215E+11
++++++++++++++++++++++++++++++++++++++++++++++++
Mass balance 3 Carbon
++++++++++++++++++++++++++++++++++++++++++++++++
sumini summas fluxout fluxin fracmass
Carbon 3.3025E+17 3.2867E+17 3.0812E+13 3.5453E+13 9.9597E-01
ifam totddep totwdep totem
3 7.980E+12 2.539E+14 5.656E+11
++++++++++++++++++++++++++++++++++++++++++++++++

[gitpeterwind]

It seems you have large amount of NHx at the start ("init"=1.36E+13) and end ("sum_mass"=2.23E+12+1.30E+08) of the run, much more than what is emitted. For NH4+NH3:
ini+emis = 1.36E+13 + 6.17E+10 = 13661700000000
sum_mass+dep = 2.23E+12+1.30E+08 + 2.77E+10+1.04E+10+3.50E+11+1.11E+13 =13718230000000
That is ok.

[mvieno]

@gitpeterwind, interesting that with that amount of NH4 in "init" it does not have a large effect. We start the model from the 1st of Jan with whatever default NH4, NO3 and SO4 value are in the EMEP model. I done tests switching off anthropogenic NH3 emission and the model results are as expected. There is some NH3 emitted from forest fires.

[mifads]

Using your budget terms, I get a fracmass of 1.004 for reduced nitrogen, so 0.4% discrepancy. In python3:

import numpy as np
#Spec         usedMW      emis     ddep     wdep     init sum_mass  fluxout   fluxin frac_mass
nh3 = np.array([ 17., 6.17E+10,2.77E+10,1.04E+10,0.00E+00,1.30E+08,2.37E+05,0.00E+00,6.20E-01 ])
nh4 = np.array([ 18., 0.00E+00,3.50E+11,1.11E+13,1.36E+13,2.23E+12,1.06E+08,2.50E+08,1.01E+00 ])

RDNterms = 14.0/17*nh3 + 14.0/18*nh4  # budget for N in kg
mw, emis, ddep, wdep, init, sum_mass, fluxout, fluxin, frac_mass = RDNterms
fracmass = (ddep + wdep + fluxout + sum_mass ) / ( init + emis + fluxin )
print(fracmass)

=> 1.00402512481

So, much better, though a little worse than the fracmass found for the Sulphur family in your results above. This still leaves the mystery of why the total Nitrogen family has 1.016.

[mifads]

@gitpeterwind

It seems you have large amount of NHx at the start ("init"=1.36E+13) and end ("sum_mass"=2.23E+12+1.30E+08) of the run, much more than what is emitted. For NH4+NH3:
ini+emis = 1.36E+13 + 6.17E+10 = 13661700000000
sum_mass+dep = 2.23E+12+1.30E+08 + 2.77E+10+1.04E+10+3.50E+11+1.11E+13 =13718230000000

According to the header of MassBudget, the terms are in kg, so to balance the N in reduced N one needs to multiply by 14/17 for NH3 and 14/18 for NH4.

[gitpeterwind]

There is clearly something wrong in the initial values. You don't use any nesting (IC) options for example?
If I try a global run (1 day) I get:
84 NH3 17.0 1.2189E+08 2.4438E+07 2.2990E+07 0.0000E+00 5.1878E+07 4.5391E-03 0.0000E+00 8.1471E-01
87 NH4_f 18.0 0.0000E+00 3.3079E+06 4.5922E+07 1.5252E+08 1.2779E+08 7.7454E+03 1.3829E+04 1.1606E+00
i.e. 1.5252E+08 for the intial value. That is more reasonable and a factor 1e5 smaller than your value.(I used latest version though)

[mifads]

And I just checked some rv4.30 results:

 # Mass Budget. Units are kg with MW used here. ADJUST! if needed
 #n    Spec        usedMW        emis        ddep        wdep        init    sum_mass     fluxout      fluxin   frac_mass
 87 NH3              17.0  6.3141E+10  1.5472E+10  2.4050E+10  0.0000E+00  9.0604E+07  1.4785E+02  0.0000E+00  6.2737E-01
 90 NH4_f            18.0  0.0000E+00  3.2576E+09  2.2072E+10  1.4948E+08  1.2319E+08  1.1206E+08  1.8132E+08  7.7282E+01

Same type of init as in Peter's latest rv4.42.

[mvieno]

@mifads , @gitpeterwind . Many thanks! So as we use the 4.36 or 4.34 it must be the way we run the model. @gitpeterwind I have not setup any specific IC. Hers is my namelist.

config_emep.txt

[mvieno]

These are the only calls I get about BIC in the logs file:

BCs_Mybcmap: CH4 settings (iyr,nml,ch4,trend): 2015 -999 -1 1870.0 1.051
rdCMX: START /Common_4.34/ZCM_EmChem19//CMX_BoundaryConditions.txt
rdCMX:# CMX_BoundaryConditions.txt for EmChem19a
rdCMX:# Provides mapping betweeen default boundary and initial conditions
rdCMX:# (BICs) and EMEP species. A numerical factor may be applied, and a wanted
rdCMX:# column is also provided so that BICs may we switched off. In many cases
rdCMX:# one may prefer to provide results from another model as BICs; in this
rdCMX:# case the mapping here is not used.
rdCMX:# Provided as part of the default GenChem system.
rdCMX:# For list of possible BICs, see defBICs in BoundaruConditions_mod.f90.
rdCMX:# For SeaSalt and Dust, see extra_mechanisms
rdCMX:#
rdCMX:# globBC emep fac wanted
rdCMX:#
rdCMX:SET: 'O3 ','O3 ' , 1.0, T => 18 19 1.000 T
rdCMX:SET: 'HNO3 ','HNO3 ' , 1.0, T => 14 27 1.000 T
rdCMX:SET: 'SO2 ','SO2 ' , 1.0, T => 1 70 1.000 T
rdCMX:SET: 'SO4 ','SO4 ' , 1.0, T => 2 104 1.000 T
rdCMX:SET: 'PAN ','PAN ' , 1.0, T => 5 66 1.000 T
rdCMX:SET: 'CO ','CO ' , 1.0, T => 6 29 1.000 T
rdCMX:SET: 'C2H6 ','C2H6' , 1.0, T => 10 31 1.000 T
rdCMX:SET: 'C4H10 ','NC4H10 ' , 1.0, T => 11 32 1.000 T
rdCMX:SET: 'NO ','NO ' , 1.0, T => 3 20 1.000 T
rdCMX:SET: 'NO2 ','NO2 ' , 1.0, T => 4 21 1.000 T
rdCMX:SET: 'HCHO ','HCHO ' , 1.0, T => 12 41 1.000 T
rdCMX:SET: 'CH3CHO ','CH3CHO ' , 1.0, T => 13 42 1.000 T
rdCMX:SET: 'NO3_f ','NO3_f ' , 1.0, T => 15 106 1.000 T
rdCMX:SET: 'NO3_c ','NO3_c ' , 1.0, T => 16 107 1.000 T
rdCMX:SET: 'NH4_f ','NH4_f ' , 1.0, T => 17 108 1.000 T
rdCMX:SET: 'H2O2 ','H2O2 ' , 1.0, T => 19 25 1.000 T
rdCMX:# CMX_BoundaryConditions.txt for generic emep run, Sea-salt
rdCMX:# globBC emep fac wanted
rdCMX:#
rdCMX:SET: 'SeaSalt_f ','SeaSalt_f', 1.0, T => 7 115 1.000 T
rdCMX:SET: 'SeaSalt_c ','SeaSalt_c', 1.0, T => 8 116 1.000 T
rdCMX:# CMX_BoundaryConditions.txt for generic emep, dust
rdCMX:# globBC emep fac wanted
rdCMX:#
rdCMX:SET: 'DUST_f', 'Dust_wb_f', 1.0, T #Dust => 21 75 1.000 T
rdCMX:SET: 'DUST_c', 'Dust_wb_c', 1.0, T #Dust => 20 76 1.000 T
rdCMX: DONE nposs, nused 20 20
BC:trends O3,CO,VOC,SOx,NOx,NH3: 2015 1.000 1.000 1.000 0.2033 0.3928 0.9230
BC: O3 Mace Head correction for year 2015
WARNING SET LATFUNC TO CONSTANT 1
reading IBC for O3 from ./Logan_P.nc
Mace Head correction for O3, trend and Mace Head value 8.764 1.000 41.000
coarse DUST BIC read from climatological file
Changing hyam from hPa to Pa in ./Dust.nc
fine DUST BIC read from climatological file

[gitpeterwind]

Actually in the first post of this issue, the file for MassBudgetSummary_2015.txt show correct initial values. So there must be something special with your last run?

[mvieno]

@gitpeterwind, ah sorry I may have grab the wrong file. I need to check with Yao. The original post mass balance did use EQSAM, where my post used MARS (I need to double check this).

[gitpeterwind]

... but anyway there is also a discrepancy between emitted and deposited N in the first post values :(

[mvieno]

@gitpeterwind the large INIT was probably the north pole (with WRF landcover 0), which was not excluded in my run.

[gitpeterwind]

Ah, yes, wrf would set the mapping factor corresponding to an area almost infinite! (it is corrected for in the latest version). That may give such results.