Dry deposition: Difference between revisions

From Geos-chem
Jump to navigation Jump to search
Line 779: Line 779:
<blockquote>In module <tt>GeosCore/drydep_mod.F</tt>, I discovered a bug in the subroutine DEPVEL. For the calculation of aerodynamic resistance RA under very stable atmospheric conditions (line 1811), the integration of the stability function phi_h (aka dimensionless vertical temperature gradient) from the roughness length to grid box center doesn't take into account the discontinuity occurring at z/L = 1 where phi_h switches from 1 + 5(z/L) to 5 + z/L  (Holtslag, 1993). The result is too high of a value for the integral of phi_h and subsequently RA.  
<blockquote>In module <tt>GeosCore/drydep_mod.F</tt>, I discovered a bug in the subroutine DEPVEL. For the calculation of aerodynamic resistance RA under very stable atmospheric conditions (line 1811), the integration of the stability function phi_h (aka dimensionless vertical temperature gradient) from the roughness length to grid box center doesn't take into account the discontinuity occurring at z/L = 1 where phi_h switches from 1 + 5(z/L) to 5 + z/L  (Holtslag, 1993). The result is too high of a value for the integral of phi_h and subsequently RA.  


A straightforward solution is to calculate RA under stable conditions (L>0) using an integral form of the stability function, such as equation (12) of Holtslag, 1988. I implemented these changes in drydep_mod.F (from GCv 10-01, public release) and have attached the updated module; locations of code change are flagged with 'blb'. To see the effect of the bug, namely a discontinuity in RA, refer to '''Figure 1''', which plots RA as a function of stability (1/L).  
A straightforward solution is to calculate RA under stable conditions (L>0) using an integral form of the stability function, such as equation (12) of Holtslag, 1988. I implemented these changes in drydep_mod.F (from GCv 10-01, public release). To see the effect of the bug, namely a discontinuity in RA, refer to '''Figure 1''', which plots RA as a function of stability (1/L).  


[[Image:Depvel_figure_1.png]]
[[Image:Depvel_figure_1.png]]

Revision as of 18:27, 17 February 2017

This page describes the current dry deposition scheme used in GEOS-Chem.

Overview

Here is a description of the GEOS-Chem dry deposition scheme from several journal articles:

From Section 4 of Alexander et al [2005]:

Dry deposition velocities for sea-salt aerosols (and sulfate formed in sea-salt aerosols) are computed with the size-dependent scheme of Zhang et al. [2001] integrated over each model size bin and accounting for hygroscopic growth as a function of relative humidity [Gerber, 1985]. Dry deposition velocities for all other species are computed with a standard resistance-in-series scheme based on Wesely [1989] as described by Wang et al. [1998].

From Section 2.4 of Bey et al [2001]:

Dry deposition of oxidants and water soluble species is computed using a resistance-in-series model based on the original formulation of Wesely [1989] with a number of modifications [Wang et al 1998]. The dry deposition velocities are calculated locally using GEOS data for surface values of momentum and sensible heat fluxes, temperature, and solar radiation.

From Section 6 of Wang et al [1998]:

We use a resistance-in-series model [Wesely and Hicks, 1977] to compute dry deposition velocities of O3, NO2, HNO3, PANs and H2O2. The deposition velocity Vi for species i is computed as:

     Vi = 1 / ( Ra + Rb,i + Rc,i )

where Ra is the aerodynamic resistance to transfer to the surface, Rb,i is the boundary resistance, and Rc,i is the canopy surface resistance. Ra and Rb,i are calculated from the GCM meteorological variables [Jacob et al 1993]. Surface resistances Rc,i are based largely on the canopy model of Wesely [1989] with some improvements, including explicit dependence of canopy stomatal resistances on LAI [Gao and Wesely, 1995] and on direct and diffuse PAR within the canopy [Baldocchi et al, 1987]. The same radiative transfer model for direct and diffuse PAR in the canopy is used as in the formulation of isoprene emissions. Surface resistances for deposition to tropical rain forest and tundra are taken from Jacob and Wofsy [1990] and Jacob et al [1992], respectively. The surface resistance for deposition of NO2 is taken to be the same as that of ozone [Erisman & Pul, 1994; Kramm et al, 1995; Eugster and Hesterberg, 1996] and hence lower than specified by Wesely [1989]. Dry deposition of CO and hydrocarbons is negligibly small and not included in the model [Mueller and Brasseur, 1995].

Aerosol dry deposition

Havala Pye reexamined the code, and she found that there are three routines basically using the Zhang et al.(2001) scheme:

  1. AERO_SFCRSII: Aerodynamic resistance for seasalt tracers. Hygroscopic growth is accounted for. Implemented by Rokjin Park ~ 2004? Used for: SALA, SALC, SO4S, NITS
  2. DUST_SFCRSII: Aerodynamic resistance of dust aerosol tracers. No hygroscopic growth. Implemented by bec ~2005? Used for: DST1, DST2, DST3, DST4
  3. ADUST_SFCRSII: Aerodynamic resistance of non-size resolved aersosols. No hygroscopic growth. Based on DUST_SRFCRSII and activated by Pye 2007. Assume particle diameter is 0.5 microns, density is 1.5 g/cm3 (1500 kg/m3). Used for all other aerosols that are not listed above (aerosols are indicated by a T/F flag in INIT_DRYDEP called AIROSOL).

This was briefly mentioned in Pye et al. (2009) paper.

Reference

  1. Pye, H. O. T., H. Liao, S. Wu, L. J. Mickley, D. J. Jacob, D. K. Henze, and J. H. Seinfeld (2009), Effect of changes in climate and emissions on future sulfate-nitrate-ammonium aerosol levels in the United States, Journal of Geophysical Research: Atmospheres, 114(D1), D01205, doi:10.1029/2008JD010701.
  2. Zhang, L. M., S. L. Gong, J. Padro, and L. Barrie (2001), A size-segregated particle dry deposition scheme for an atmospheric aerosol module, Atmos. Environ., 35(3), 549-560, doi:10.1016/s1352-2310(00)00326-5.

Source code and data

The source code of the various routines are located in drydep_mod.f.

  1. Routine DO_DRYDEP is the main driver for the dry deposition. This calls various setup subroutines.
  2. Routine DEPVEL (which is called by DO_DRYDEP) computes dry deposition dry deposition velocities in [m/s].
  3. Dry deposition frequencies [1/s] for all species are computed in DO_DRYDEP, after the call to DEPVEL.

For more information about the input data files required by the drydep scheme, please see the README file: GEOS_1x1/drydep_200203/README

In GEOS-Chem v10-01 and higher versions

NOTE: This update was made to facilitate the implementation of both the HEMCO emissions component and FlexChem chemical solver.

In GEOS-Chem v10-01 and higher versions, the following changes were made to how dry deposition is applied:

  1. Dry deposition is now applied to all tracers in module GeosCore/mixing_mod.F90. Prior to this, dry deposition had been applied in several locations throughout the code, which was confusing and prone to error.
  2. Dry deposition fluxes (i.e. category DRYD-FLX) are now archived to the ND44 diagnostic in module GeosCore/mixing_mod.F90. Furthermore, drydep fluxes are now archived with units of molec/cm2/s for all GEOS-Chem simulations. In prior GEOS-Chem versions, some simulations (in particular Rn-Pb-Be) had archived drydep fluxes in units of kg/s. Therefore, when analyzing output from GEOS-Chem v10-01 (or higher versions), please make sure to modify your plotting or budget programs to expect drydep fluxes in units of molec/cm2/s.
  3. Routines DRYFLX and DRYFLXRnPbBe have now been removed from module GeosCore/drydep_mod.F. The equivalent functionality is now provided by GeosCore/mixing_mod.F90.

--Bob Y. (talk) 21:34, 17 June 2015 (UTC)

In GEOS-Chem v9-02 and earlier versions

NOTE: The information below is now obsolete as of GEOS-Chem v10-01. We shall leave the sections below for reference.

Full-chemistry simulation

In the NOx-Ox-HC-aerosol "full-chemistry" simulation, dry deposition is handled as follows:

  1. The DEPVEL array is passed to routine CALCRATE (in source code file calcrate.f). The dry deposition frequencies of the gas-phase species (NO2, O3, PAN, HNO3, H2O2, PMN, PPN, R4N2, CH2O, N205, etc.) are then entered into the master reaction rate array, RRATE.
  2. The RRATE is passed to either the SMVGEAR or KPP solver (whichever you have selected). The dry deposition frequencies are considered as loss rates in the chemical mechanism.
  3. Dry deposition of the sulfur aerosol tracers (SO2, SO4, SO4s, MSA, NH3, NH4, NIT, NITs) is applied in sulfate_mod.f. Deposited tracer is removed directly from the main tracer array (STT).
  4. Dry deposition of the carbon and/or SOA aerosol tracers (BCPI, BCPO, OCPI, OCPO, etc.) is applied in carbon_mod.f. Deposited tracer is removed directly from the main tracer array (STT).
  5. Dry deposition of the mineral dust tracers (DST1, DST2, DST3, DST4) is applied in dust_mod.F. Deposited tracer is removed directly from the main tracer array (STT).
  6. Dry deposition of the sea salt aerosol tracers (SALA, SALC) is applied in dust_mod.F. Deposited tracer is removed directly from the main tracer array (STT).

Offline simulations

In general, dry deposition for offline simulations is handled in the relevant module (e.g. RnPbBe_mod.F for the Radon-Lead-Beryllium simulation, etc. The dry deposition frequencies stored in DEPVEL are used to remove dry deposited tracer directly from the main tracer array (STT).

--Bob Y. 15:34, 21 October 2011 (EDT)

Input values for dry deposition

Here we list some of the various input values used by the dry deposition routines. Full citations may be found in the references section below.

Land cover parameters

The land cover input variables for the GEOS-Chem dry deposition module are listed in the table below:

Variable Read from Description Values Reference
DRYCOEFF
  • v9-01-03 and higher:
    Olson_1992_Drydep_inputs.nc

  • Prior to v9-01-03:
    drydep.coef
Local dependence of stomal resistance on the intensity I of light impinging the leaf; this is expressed as a mutliplicative factor I/(I+b) to the stomatal resistance where b = 50 W m-2 (equation (7) of Baldocchi et al. [1987]) -0.358, 3.02, 3.85, -0.0978, -3.66, 12, 0.252, -7.8, 0.226, 0.274, 1.14, -2.19, 0.261, -4.62, 0.685, -0.254, 4.37, -0.266, -0.159, -0.206
  • Baldocchi et al, AE, 1987.
IDRYDEP
  • v9-01-03 and higher:
    Olson_1992_Drydep_inputs.nc

  • Prior to v9-01-03:
    drydep.table
Indices for the 11 dry deposition land types 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
  • Wesely, AE 1989;
  • Jacob et al., JGR 1990, 1992.
IOLSON
  • v9-01-03 and higher:
    Olson_1992_Drydep_inputs.nc

  • Prior to v9-01-03:
    drydep.table
Indices for the 72 Olson land types 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74
  • Wesely, AE 1989;
  • Jacob et al., JGR 1990, 1992.
IDEP
  • v9-01-03 and higher:
    Olson_1992_Drydep_inputs.nc

  • Prior to v9-01-03:
    drydep.table
Indices that map each of the 72 Olson land types to the corresponding dry deposition land type 11, 10, 5, 1, 1, 1, 2, 1, 8, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1, 3, 3, 3, 3, 2, 2, 2, 3, 2, 2, 4, 4, 2, 6, 1, 1, 9, 4, 4, 4, 5, 5, 5, 5, 5, 9, 5, 5, 5, 5, 8, 8, 5, 7, 6, 2, 2, 2, 2, 2, 3, 3, 3, 5, 5, 11, 11, 11, 11, 8, 1, 8, 9, 11
  • Wesely, AE 1989;
  • Jacob et al., JGR 1990, 1992.
IWATER
  • v9-01-03 and higher:
    Olson_1992_Drydep_inputs.nc

  • Prior to v9-01-03:
    drydep.table
Indices of the Olson land types denoting water (i.e. ocean/lake/ice) surfaces 1, 66, 67, 68, 69, 74
  • Wesely, AE 1989;
  • Jacob et al., JGR 1990, 1992.
IZO
  • v9-01-03 and higher:
    Olson_1992_Drydep_inputs.nc

  • Prior to v9-01-03:
    drydep.table
Roughness heights (1e-4 m) for each of the 72 Olson land types 10, 25000, 100, 1000, 1000, 1000, 10000, 1000, 10, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1, 1000, 1000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 1000, 10000, 1000, 1000, 2000, 10000, 1000, 1000, 100, 1000, 1000, 1000, 100, 100, 100, 100, 100, 100, 1000, 1000, 1000, 1000, 10, 10, 100, 50, 10000, 2000, 2000, 2000, 2000, 2000, 10000, 10000, 10000, 2000, 50, 100, 100, 100, 100, 10, 1, 1, 500, 10
  • Wesely, AE 1989;
  • Jacob et al., JGR 1990, 1992.

--Bob Yantosca (talk) 17:43, 21 July 2016 (UTC)

Aerodynamic resistances

The aerodynamic resistances (IRI, IRLU, IRAC, IRGSS, IRGSO, IRCLS, IRCLO) and maximum deposition velocity for aerosol (IVSMAX) for each of the 11 dry deposition land types are read from this file:

  • v9-01-03 and higher: Olson_1992_Drydep_inputs.nc
  • Prior to v9-01-03: drydep.table

and have the following values:

DD type Description IRI
[s m-1]
IRLU
[s m-1]
IRAC
[s m-1]
IRGSS
[s m-1]
IRGSO
[s m-1]
IRCLS
[s m-1]
IRCLO
[s m-1]
IVSMAX
[10-2 cm s-1]
Reference
1 Snow/Ice 9999 9999 0 100 3500 9999 1000 100 Wesely, AE, 1989
2 Deciduous forest 200 9000 2000 500 200 2000 1000 100 Wesely, AE, 1989
3 Coniferous forest 400 9000 2000 500 200 2000 1000 100 Wesely, AE, 1989
4 Agricultural land 200 9000 200 150 150 2000 1000 100 Wesely, AE, 1989
5 Shrub/grassland 200 9000 100 350 200 2000 1000 100 Wesely, AE, 1989
6 Amazon forest 200 1000 2000 200 200 9999 9999 100 Jacob & Wofsy, JGR 1990
7 Tundra 200 4000 0 340 340 9999 9999 100 Jacob et al., JGR 1992
8 Desert 9999 9999 0 1000 400 9999 9999 10 Wesely, AE, 1989
9 Wetland 200 9000 300 0 1000 2500 1000 100 Wesely, AE, 1989
10 Urban 9999 9999 100 400 300 9999 9999 100 Wesely, AE, 1989
11 Water 9999 9999 0 0 2000 9999 9999 10 Wesely, AE, 1989

--Bob Yantosca (talk) 17:43, 21 July 2016 (UTC)

OVOCs dry deposition

These values were further changed as part of the new isoprene scheme. The most up-to-date values are on the new isoprene scheme page.

--Jmao, JAF 19 March 2014

This update was tested in the 1-month benchmark simulation v9-01-03c and approved on 19 Dec 2011.

In GEOS-Chem v11-01 and higher versions, these parameters are now defined in the GEOS-Chem species database, and are then copied into the local dry deposition arrays in routine INIT_DRYDEP.

According to Karl et al.(2010), I changed the reactivity factor(f0) for all OVOCs from 0 to 1.

Now dry dep is set up as:

Species H*(moles L-1 atm-1 ) f0 Reference
NO2 0.01 0.1
Ox 0.01 1.0
PAN 3.6 0.1
HNO3 1.0d+14 0.0
H2O2 1.0d+5 1.0
PMN as PAN
PPN as PAN
PYPAN as PAN
ISN2 as HNO3
R4N2 as PAN
CH2O 6.0e+3 change from 0 to 1.0
GLYX 3.6d+5 change from 0 to 1.0
MGLY 3.7d+3 change from 0 to 1.0
GLYC 4.1d+4 change from 0 to 1.0
MPAN,GPAN, GLPAN as PAN
N2O5 as HNO3
HCOOH 1.67d+5 change from 0 to 1.0
ACTA 1.14d+4 change from 0 to 1.0
ISOPND 1.7d+4 0.0
ISOPNB 1.7d+4 0.0
MVKN+MACRN 1.7d+4 0.0
PROPNN 1.0d+3 0.0 NITROOXYACETONE IN SANDER TABLE
RIP as H2O2
IEPOX as H2O2
MAP 8.4d+2 1.0
MVK 4.4d1 change from 0 to 1.0 from R.Sander
MACR 6.5d0 change from 0 to 1.0 from R.Sander
SO2 1.0d+5 0.0

Karl, T., Harley, P., Emmons, L., Thornton, B., Guenther, A., Basu, C., Turnipseed, A., and Jardine, K.: Efficient Atmospheric Cleansing of Oxidized Organic Trace Gases by Vegetation, Science, 330, 816-819, 10.1126/science.1192534, 2010.

Notes: To run with non-local mixing, you don't need to modify globchem.dat with " MVK=DRYDEP+DRYMVK", because drydep here is done out of chemistry step. But you need to modify that for full mixing scheme, which is using chemical solver to calculate drydep.

--Jmao 15:06, 15 December 2011 (EST)

Code updates

Updates to dry deposition when using the Olson 2001 land map

This update was validated in the 1-month benchmark simulation GEOS-Chem v10-01a and approved on 05 Feb 2014.

It was discovered that ozone in the eastern US in GEOS-Chem v9-02 was too high when using the Olson 2001 land map. Patrick Kim implemented the following updates to correct this issue:

1. Force the use of MODIS LAI from 2008 if the last year is beyond 2008. There is a large difference in the 2009 file that still needs to be investigated.
In routine Init_Modis_Lai (modis_lai_mod.F90):
   IF ( USE_OLSON_2001 ) THEN
      I_MODIS     = 1440             ! For Olson 2001, use MODIS LAI
      J_MODIS     = 720              ! on the 0.25 x 0.25 native grid
      MODIS_START = 2005             ! First year of MODIS data  
      MODIS_END   = 2008             ! Force to 2008 (skim, 1/29/14)
2. Change the internal resistance for coniferous forests to match the internal resistance for deciduous forests.
In routine READ_DRYDEP_INPUTS (drydep_mod.F):
     !----------------------------------------
     ! VARIABLE: IRI
     !----------------------------------------

     ! Variable name
     v_name = "IRI"

     ! Read IRI from file
     st1d   = (/ 1         /)
     ct1d   = (/ NDRYDTYPE /)
     CALL NcRd( IRI, fId, TRIM(v_name), st1d, ct1d )

     ! Read the IRI:units attribute
     a_name = "units"
     CALL NcGet_Var_Attributes( fId,TRIM(v_name),TRIM(a_name),a_val )

     ! If using Olson 2001 land map, replace IRI for coniferous forests
     ! with IRI for deciduous forests (skim, mps, 2/3/14)
     IF ( USE_OLSON_2001 ) THEN
        IRI(3) = 200
     ENDIF

     ! Echo info to stdout
     IF ( am_I_Root ) THEN
        WRITE( 6, 130 ) TRIM(v_name), TRIM(a_val)
     ENDIF

--Melissa Sulprizio 11:31, 3 February 2014 (EST)

Updates to facilitate grid-independent operation

This update was tested in the 1-month benchmark simulation v9-01-03j and approved on 17 Apr 2012.

A number of structural updates were made to the dry deposition routines in order to remove bottlenecks that would have prevented grid-independent functionality. These include:

  1. Introduction of a new module to read Olson land map data from native resolution and to regrid to the currently-selected model resolution. The regridding is done in such a way so as to preserve backwards compatibility and to populate the existing arrays (IREG, ILAND, IUSE, FRCLND) in the same manner as before.
  2. Introduction of a new module to read MODIS LAI data from native resolution and regrid to the current resolution in a consistent manner.
  3. Removal of obsolete variables:
    1. From Headers/CMN_VEL_mod.F: IJREG, IJUSE, IJLAND
    2. From GeosCore/drydep_mod.F: NNNTYPE, NNNVEGTYPE
  4. Replacement of obsolete ASCII input files with netCDF files.

--Bob Y. 15:35, 17 December 2012 (EST)

Update dry deposition to use local surface pressure

This update was tested in the 1-month benchmark simulation v9-01-03d and approved on 12 Jan 2012.

Lyatt Jaeglé wrote:

The only difference comes from the pressure value (PRESS variable) used in the dry deposition code to calculate the diffusivity of a molecule in air. Pressure comes into play in the mean free path calculation in DIFFG. For some reason the pressure was set at a constant value of 1500hPa (I do not know why this value was chosen in the first place, any ideas?). When I was going over the dry deposition code for sea salt I had replaced this fixed value with local sea level pressure (SLP). I think that it makes sense to use the actual local pressure instead of this fixed value. The overall effect is small.

--Melissa Payer 11:19, 11 January 2012 (EST)

Aerosol dry deposition velocities over snow and ice surfaces

This update was tested in the 1-month benchmark simulation v9-01-02i and approved on 15 Aug 2011.

Modeled aerosol dry deposition velocities over snow and ice surfaces in the Arctic are much higher than estimated from measured values (e.g., Ibrahim et al. [1983]; Duan et al. [1988]; Nilsson and Rannik [2001]). In GEOS-Chem v9-01-02 we have imposed a dry deposition velocity of 0.03 cm/s for all aerosols over snow and ice surfaces.

This update affects the following aerosol tracers: SO4s, NITs, DST1, DST2, DST3, DST4, SALA, SALC.

--Jenny Fisher 11:08, 20 January 2011 (EST)
--Bob Y. 10:26, 16 August 2011 (EDT)

Removal of CMN_DEP_mod.F

This update was validated with 1-month benchmark simulation v9-02d and approved on 19 Dec 2012.

As described in this wiki post, we have removed the module file Headers/CMN_DEP_mod.F from GEOS-Chem. The dry deposition arrays that were contained in this module are now carried by the the Meteorology State object (aka State_Met).

--Bob Y. 11:40, 4 March 2013 (EST)

Application of dry deposition with mixing_mod.F90

This discussion has been moved to a section above.

--Bob Yantosca (talk) 18:42, 17 February 2016 (UTC)

References

  1. Alexander, B., R.J. Park, D.J. Jacob, Q.B. Li, R.M. Yantosca, J. Savarino, C.C.W. Lee, and M.H. Thiemens, Sulfate formation in sea-salt aerosols: Constraints from oxygen isotopes, J. Geophys. Res., 110, D10307, 2005. PDF
  2. Baldocchi, D.D., B.B. Hicks, and P. Camara, A canopy stomatal resistance model for gaseous deposition to vegetated surfaces, Atmos. Environ. 21, 91-101, 1987.
  3. Brutsaert, W., Evaporation into the Atmosphere, Reidel, 1982.
  4. Businger, J.A., et al., Flux-profile relationships in the atmospheric surface layer, J. Atmos. Sci., 28, 181-189, 1971.
  5. Bey I., D. J. Jacob, R. M. Yantosca, J. A. Logan, B. Field, A. M. Fiore, Q. Li, H. Liu, L. J. Mickley, and M. Schultz, Global modeling of tropospheric chemistry with assimilated meteorology: Model description and evaluation, J. Geophys. Res., 106, 23,073-23,096, 2001. PDF
  6. Corbitt, E. S., How does the current implementation of dry deposition work in GEOS-Chem?, Graduate Student Forum Talk, 20 Oct 2011. PDF
  7. Duan, B., Fairall, C.W., and Thomson, D.W., Eddy correlation measurements of the dry deposition of particles in wintertime, J. Appl. Met., 27, 642-652, 1988.
  8. Dwight, H.B., Tables of integrals and other mathematical data, MacMillan, 1957.
  9. Erisman, J.W, and A.V. Pul, Parameterization of surface resistance for the quantification of atmospheric deposition of acidifying pollutants and ozone, Atmos. Environ, 28, 2595-2607, 1994.
  10. Eugster, W. and R. Hesterberg, Transfer resistances of NO2 determined from eddy correlation flux measurements over a litter meadow at a rural site on the Swiss plateau, Atmos. Environ, 30, 1247-1254, 1996.
  11. Gao, W. and M.L. Wesely, Modeling gaseous dry deposition over regional scales with satellite observations, 1. Model development, Atmos. Environ, '29, 727-737, 1995.
  12. Gerber, H. E. (1985), Relative-humidity parameterization of the Navy aerosol model (NAM), NRL Rep. 8956, Natl. Res. Lab., Washington, D. C., 1985.
  13. Guenther, A., et al, A global model of natural volatile organic compound emissions, J. Geophys. Res., 100, 8873-8892, 1995.
  14. Hicks, B.B., and P.S. Liss, Transfer of SO2 and other reactive gases across the air-sea interface, Tellus, 28, 348-354, 1976.
  15. Ibrahim, M., Barrie, L.A., and Fanaki, F., An experimental and theoretical investigation of the dry deposition of particles to snow, pine trees, and artificial collectors, Atmos. Environ., 17, 781-788, 1983.
  16. Jacob, D.J., and S.C. Wofsy, Budgets of reactive nitrogen, hydrocarbons, and ozone over the Amazon forest during the wet season, J. Geophys. Res., 95, 16737-16754, 1990.
  17. Jacob, D.J., et al, Deposition of ozone to tundra, J. Geophys. Res., 97, 16473-16479, 1992.
  18. Jacob, D.J, et al, Summertime photochemistry of the troposphere at high northern latitudes, J. Geophys. Res., 98, 14,797-14,816, 1993.
  19. Kramm, G. R. Dlugi, G.J. Dollard, T. Foken, N. Moelders, H. Mueller, W. Seiler, and H. Sievering, On the dry deposition of ozone and reactive nitrogen species, Atmos. Environ, 29, 3208-3231, 1995.
  20. Levine, I.N., Physical Chemistry, 3rd ed., McGraw-Hill, New York, 1988.
  21. Mueller, J-F, and G. Brasseur, IMAGES, A three-dimensional chemical transport model of the global troposphere, J. Geophys. Res., 100, 16,445-16,490, 1995.
  22. Munger, J.W., et al, Atmospheric deposition of reactive nitrogen oxides and ozone in a temperate deciduous forest and a sub-arctic woodland, J. Geophys. Res., in press, 1996.
  23. Nilsson, E.D. and Rannik, U. Turbulent aerosol fluxes over the Arctic Ocean: 1. Dry deposition over sea and pack ice, J. Geophys. Res, 106, 32,125-32,137, 2001.
  24. Price, H., L. Jaeglé, A. Rice, P. Quay, P.C. Novelli, R. Gammon, Global Budget of Molecular Hydrogen and its Deuterium Content: Constraints from Ground Station, Cruise, and Aircraft Observations, submitted to J. Geophys. Res., 2007
  25. Pye, H. O. T., H. Liao, S. Wu, L. J. Mickley, D. J. Jacob, D. K. Henze, and J. H. Seinfeld (2009), Effect of changes in climate and emissions on future sulfate-nitrate-ammonium aerosol levels in the United States, J. Geophys. Res., 114(D1), D01205, doi:10.1029/2008JD010701.
  26. Walcek, C.J., R.A. Brost, J.S. Chang, and M.L.Wesely, SO2, sulfate, and HNO3 deposition velocities computed using regional landuse and meteorological data, Atmos. Environ., 20, 949-964, 1986.
  27. Wang, Y., D.J. Jacob, and J.A. Logan, Global simulation of tropospheric O3-NOx-hydrocarbon chemistry, 1. Model formulation, J. Geophys. Res., 103, D9,10,713-10,726, 1998. PDF
  28. Wesely, M.L, Improved parameterizations for surface resistance to gaseous dry deposition in regional-scale numerical models, Environmental Protection Agency Report EPA/600/3-88/025, Research Triangle Park (NC), 1988.
  29. Wesely, M. L., Parameterization of surface resistance to gaseous dry deposition in regional-scale numerical models, Atmos. Environ., 23, 1293-1304, 1989.
  30. Wesely, M.L, and B.B. Hicks, Some factors that affect the deposition rates of sulfur dioxide and similar gases on vegetation, J. Air. Pollut. Control Assoc., 27, 1110-1116, 1977.
  31. Zhang, L. M., S. L. Gong, J. Padro, and L. Barrie (2001), A size-segregated particle dry deposition scheme for an atmospheric aerosol module, Atmos. Environ., 35(3), 549-560, doi:10.1016/s1352-2310(00)00326-5.

--Bob Y. 15:24, 19 February 2010 (EST)

Previous issues that have now been resolved

Dependency between dry deposition and soil NOx emissions

This issue has been resolved in the new soil NOx algorithm that was added to GEOS-Chem v9-02. The new module GeosCore/canopy_nox_mod.F computes the NOx from the canopy independently of the dry deposition module. We shall leave this post here for reference.

The canopy reduction algorithm (originally done in GeosCore/canopy_nox_mod.F) has been moved to the HEMCO soil NOx extension module (HEMCO/Extensions/hcox_soilnox_mod.F90). Module GeosCore/canopy_nox_mod.F was removed from GEOS-Chem v10-01 and higher versions.

In GEOS-Chem there is a code dependency between the dry deposition routines soil NOx emissions routines. This is purely historical baggage that goes back to the days of the old 9-layer Harvard-GISS CTM (from which these routines were taken).

The dry deposition routine DEPVEL (in drydep_mod.f) computes a quantity called CANOPYNOX, which is the bulk surface resistance of the canopy to NOx deposition. This is computed in the following lines of code:

C** Get the bulk surface resistance of the canopy, RSURFC, from the network
C** of resistances in parallel and in series (Fig. 1 of Wesely [1989])
               DTMP1=1.D0/RIXX
               DTMP2=1.D0/RLUXX
               DTMP3=1.D0/(RAC(LDT)+RGSX)
               DTMP4=1.D0/(RDC+RCLX)
               RSURFC(K,LDT) = 1.D0/(DTMP1 + DTMP2 + DTMP3 + DTMP4)
C  Save the within canopy depvel of NOx, used in calculating the 
C  canopy reduction factor for soil emissions.
               ! Remove hardwire for CANOPYNOX (bmy, 1/24/03)
               IF ( K == DRYDNO2 ) THEN
                  CANOPYNOX(IJLOOP,LDT)=DTMP1+DTMP2+DTMP3+DTMP4
               ENDIF

However, the CANOPYNOX variable is not defined within drydep_mod.f, but is found instead within the soil NOx emissions code header file commsoil.h:

     REAL*8 CANOPYNOX(MAXIJ,NTYPE) !track NOx within canopy dry dep.

and is then also used within the soil NOx emissions function SOILCRF, as follows:

     IF ((XLAI(IREF,JREF,K).GT.0.D0).AND.
    &    (CANOPYNOX(IJLOOP,K).GT.0.D0))THEN

        VFNEW=VFNEW*SQRT(WINDSQR/9.D0*7.D0/XLAI(IREF,JREF,K))*
    *        (SOILEXC(2)/SOILEXC(NN))
        SOILCRF=CANOPYNOX(IJLOOP,K)/(CANOPYNOX(IJLOOP,K)
    *        +VFNEW)
     ELSE
    
        SOILCRF=0.D0
     END IF

This is not necessarily a problem as long as dry deposition is always turned on whenever you are doing a full-chemistry simulation that requires soil NOx emissions. In the driver routine main.f, dry deposition is always done before emissions. However, this is bad coding style and a potential source of error. (This also must be corrected for the GEOS-Chem column code, in which all inputs must come thru the argument list.)

In the new soil NOx emissions code being developed by Rynda Hudman (much of which is based on the GEOS-Chem column code being developed by Bob Yantosca and Philippe Le Sager), there is a new module called canopy_nox_mod.f. This module computes the CANOPYNOX quantity independently of routine DEPVEL in drydep_mod.f, which allows for a totally clean separation between dry deposition routines and emissions routines.

At present, the new canopy_nox_mod.f is not available in the standard mainline GEOS-Chem code. It is slated to be incorporated with the soil NOx emissions update.

--Bob Y. 13:51, 18 February 2010 (EST)

Prevent floating-point invalid in drydep_mod.F

NOTE: This fix was included in GEOS-Chem v9-02p.

In routine DEPVEL (in GeosCore/drydep_mod.F we now skip to the next species when the molecular weight for drydep species K (aka XMW(K)) is zero.

Change these lines of code:

!** If LDEP(K)=F, species does not deposit.
!** Deposition is applied only to species with LDEP=T.

      DO K = 1,NUMDEP
         LDEP(K) = (HSTAR(K).GT.0.D0 .OR. F0(K).GT.0.D0 
     &              .OR. AIROSOL(K)) 
      ENDDO

to this:

!** If LDEP(K)=F, species does not deposit.
!** Deposition is applied only to species with LDEP=T.
      DO K = 1,NUMDEP

         ! Better test for depositing species K: We need both HSTAR and XMW
         ! to be nonzero, OR the value of AIROSOL to be true.  This should
         ! avoid any futher floating point invalid issues caused by putting
         ! a zero value in a denominator. (bmy, 8/29/13)
         ! (bmy, 8/29/13)
         IF ( ( HSTAR(K) > 0d0   .and. 
     &          XMW  (K) > 0d0 ) .or.  AIROSOL(K) ) THEN
            LDEP(K) = .TRUE.
         ELSE
            LDEP(K) = .FALSE.
         ENDIF
      ENDDO

A few species have a zero molecular weight in GeosCore/drydep_mod.F to denote they use the same deposition velocities as another species (i.e. PMN uses PAN's deposition velocity).

This error condition was caught by the GEOS-Chem unit tester when compiling with option FPE=yes.

Please also view the wiki post about this related error.

--Bob Y. 11:46, 12 September 2013 (EDT)

Error in definition of RIP and IEPOX drydep species

This fix was included in GEOS-Chem v9-02p.

NOTE: The code below was removed from GEOS-Chem v11-01d. Dry deposition parameters and molecular weights are now specified with the GEOS-Chem species database in v11-01d and higher versions.

In routine INIT_DRYDEP, we now assign a molecular weight of 118e-3 kg/mole to the species RIP and IEPOX. This corrects a typo, as described the conversation below.

Bob Yantosca wrote:

When compiling GEOS-Chem while checking for floating-point errors (using compiler flags -fpe0 –ftrapuv), I found an issue in GeosCore/drydep_mod.F.
Several species have XMW(K) = 0d0 (i.e. the molecular weight that goes into routine DEPVEL). This causes an error such as:
   DEPVEL: Deposition velocity is negative!
   Dep. Vel =  -1.080384596850893E-062
   Species  =           43
   IJLOOP   =            1
   RADIAT   =   0.000000000000000E+000
   TEMP     =    216.756195068359    
   SUNCOS   =  -0.377513942098922    
   USTAR    =   0.310841470956802    
   CZ1      =    48.4691696913914    
   OBK      =    45.7151465305601    
   CFRAC    =   0.671608686447144    
   ZH       =    190.063915293881    
   LRGERA   =  F
   ZO       =   1.998901367187500E-003
   STOP in depvel.f!

Species 43 is RIP. I then looked up RIP (and its neighbor IEPOX) in INIT_DRYDEP and I found:
           !RIP
           ELSE IF ( N == IDTRIP ) THEN
              NUMDEP          = NUMDEP + 1
              NTRAIND(NUMDEP) = IDTRIP
              NDVZIND(NUMDEP) = NUMDEP
              DEPNAME(NUMDEP) = 'RIP'
              HSTAR(NUMDEP)   = 1.70d6
              F0(NUMDEP)      = 1.0d0!Karl et al. 2010(science)
              XMW(NUMDEP)     = 0d0
              AIROSOL(NUMDEP) = .FALSE.

           !IEPOX
           ELSE IF ( N == IDTIEPOX ) THEN
              NUMDEP          = NUMDEP + 1
              NTRAIND(NUMDEP) = IDTIEPOX
              NDVZIND(NUMDEP) = NUMDEP !USE H2O2
              DEPNAME(NUMDEP) = 'IEPOX'
              HSTAR(NUMDEP)   = 1.30d8
              F0(NUMDEP)      = 1.0d0 !Karl et al. 2010(science)
              XMW(NUMDEP)     = 0d0
              AIROSOL(NUMDEP) = .FALSE.

It is OK for dry deposition species to have a zero value for XMW, but only if they are assigned the dry deposition velocity from another species (e.g. PMN uses the drydep velocity of PAN). But in this context, it is not appropriate for RIP and IEPOX to have a zero value for XMW molecular weight value. How should we define XMW for RIP and IEPOX?

Fabien Paulot wrote:

I think you got it right. This is probably an oversight. When I wrote the code, I didn't have any information regarding the Henry's constant for IEPOX and RIP and used H2O2 as a proxy. Henry's constants got updated later on. Though those are very uncertain, I would suggest to keep them and set XMW to the appropriate molecular weight for RIP and IEPOX.

Eloise Marais wrote:

Here are the Henry's constants I used for my work (with references):
   RIP    = 1.7d6 M atm^-1 (US EPA, 2011)
   IEPOX  = 1.3d8 M atm^-1 (Eddingsaas et al., 2010)

--Bob Y. 15:40, 6 September 2013 (EDT)

Add extra variables to PRIVATE declaration in parallel loop

This fix was validated by 1-month benchmark simulation GEOS-Chem v9-02p (approved 13 Sep 2013).

We have explicitly added variable IW to the !$OMP+PRIVATE list in routine DEPVEL in GeosCore/drydep_mod.F. This should already be held private because it is a loop index, but it doesn't hurt to add it to the list.

We also have added ZO as a PRIVATE variable when compiling with GCAP meteorology.

Old code:
!$OMP+PRIVATE( DTMP1,  DTMP2,  DTMP3,  DTMP4,  VDS,    CZH    )
!$OMP+PRIVATE( CKUSTR, REYNO,  CORR1,  CORR2,  Z0OBK,  RA     )
!$OMP+PRIVATE( DUMMY1, DUMMY2, DUMMY3, DUMMY4, DAIR,   RB     )
!$OMP+PRIVATE( C1X,    VK,     I,      J                      )

New code:
!$OMP+PRIVATE( DTMP1,  DTMP2,  DTMP3,  DTMP4,  VDS,    CZH    )
!$OMP+PRIVATE( CKUSTR, REYNO,  CORR1,  CORR2,  Z0OBK,  RA     )
!$OMP+PRIVATE( DUMMY1, DUMMY2, DUMMY3, DUMMY4, DAIR,   RB     )
!$OMP+PRIVATE( C1X,    VK,     I,      J,      IW             )
#if defined( GCAP )
!$OMP+PRIVATE( ZO                                             )
#endif

--Bob Y. 12:04, 12 September 2013 (EDT)

Add minor fixes for ND44 dry deposition diagnostic

This fix was tested with the 1-month benchmark simulation GEOS-Chem v9-02p and approved on 13 Sep 2013.

We corrected some minor issues for the dry deposition diagnostic in the following files:

  1. GeosCore/gamap_mod.F:
    • Remove inconsistencies in the ND44 tracerinfo.dat listing.
    • Save ISOPNB as NTRACERS+1 and MVKN as NTRACERS+2.
    • Replace confusing code with more clearly defined code.
  2. GeosCore/drydep_mod.F:
    • Improve printout of drydep species to stdout.
    • Also improve commenting in INIT_DRYDEP for species ISOPND, ISOPNB, MACRN, MVKN.
  3. GeosCore/diag3.F:
    • Now save ISOPNB with index NTRACERS+1 and MVKN with index NTRACERS+2 in the ND44 diagnostic. NNPAR+1 and NNPAR+2 were previously used.

--Bob Y. 11:59, 12 September 2013 (EDT)

Parallelization error in routine METERO

This update was tested in the 1-month benchmark simulation v9-02r and approved on 14 Nov 2013. This update is included in Adjoint v35j.

We discovered a minor bug in the parallel DO loop of routine METERO (in module GeosCore/drydep_mod.F):

      ! Loop over surface grid boxes
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, IJLOOP, THIK )
      DO J = 1, JJPAR
      DO I = 1, IIPAR

         ! 1-D grid box index
         IJLOOP         = ( (J-1) * IIPAR ) + I

         ! THIK = thickness of layer 1 [m]
         THIK           = State_Met%BXHEIGHT(I,J,1)

         ! Midpoint height of first model level [m]
         CZ1(IJLOOP)    = THIK / 2.0d0

         ! Local surface pressure [hPa] (mpayer, 1/10/12)
         SP             = GET_PEDGE(I,J,1)

         ! Convert from hPa to Pa for SFCPRESS
         PRESSU(IJLOOP) = SP * 1.d2

The SP variable (which is surface pressure) had not been held PRIVATE in the above parallel loop. This caused slightly different surface pressure values to be saved into the PRESSU array. Later on in the code, PRESSU is passed as an argument to function DIFFG, which computes the computes the molecular diffusivity of each drydep species. This in turn affects the value of the drydep velocity of each species. Therefore, as a result of this error, we would observe small numerical differences in the results of GEOS-Chem simulations when using OpenMP parallelization versus when OpenMP parallelization was turned off.

The fix is to add SP to the PRIVATE list:

       ! Loop over surface grid boxes
 !$OMP PARALLEL DO
 !$OMP+DEFAULT( SHARED )
 !$OMP+PRIVATE( I, J, IJLOOP, THIK, SP )

which then ensures that we obtain identical results whether or not OpenMP parallelization is turned on.

--Bob Y. 15:25, 14 November 2013 (EST)

Parallelization error in routine DRYFLX

This update was tested in the 1-month benchmark simulation v9-02r and approved on 14 Nov 2013. This update is included in Adjoint v35j.

In routine DRYFLX (in module GeosCore/drydep_mod.F), the variable AREA_CM2 needs to be added to the !$OMP PRIVATE list for the parallel DO loop. We changed this line of code:

       !$OMP+PRIVATE( I, J, L, N, NN, NK, JJ, JLOOP, TDRYFX )

to:

       !$OMP+PRIVATE( I, J, L, N, NN, NK, JJ, JLOOP, TDRYFX, AREA_CM2 )

This error was causing the ND44 drydep diagnostic to return slightly different results when using single vs. multiple processors, under the following conditions:

  1. For all simulations using GEOS-4 meteorology, or
  2. For simulations using the full boundary layer mixing scheme (aka TURBDAY).

--Bob Y. 15:26, 14 November 2013 (EST)

Remove obsolete ISN2 drydep species

This update was validated with the 1-month benchmark simulation v10-01i and approved on Approved 01 May 2015.

Katherine Travis wrote:

In dry deposition, we still have a dry dep option for ISN2. ISN2 was the nighttime isoprene nitrate in old versions of the chemistry, and ISN1 was the second generation product from this species. But now, we call ISN1 the isoprene nitrate species. So, we should replace IDTISN2 in dry dep and tracerid_mod with IDTISN1. That fix is also in the above directory.

In GeosCore/drydep_mod.F, we replaced these lines of code:

        ! ISN2 (uses same dep vel as HNO3)
        ELSE IF ( N == IDTISN2 ) THEN
           . . .
           NTRAIND(NUMDEP) = IDTISN2
           . . .
           DEPNAME(NUMDEP) = 'ISN2'

with these:

        ! ISN1 (uses same dep vel as HNO3)
        ! This should be ISN1, krt (3/31/15)
        ELSE IF ( N == IDTISN1 ) THEN
           . . .
           NTRAIND(NUMDEP) = IDTISN1
           . . .
           DEPNAME(NUMDEP) = 'ISN1'

and in GeosCore/tracerid_mod.F, we replaced these lines of code:

     INTEGER            :: IDISN2,  IDR4N2,   IDCH2O,  IDC2H6,  IDMP
     . . .
     INTEGER            :: IDTPPN,  IDTISN1,  IDTR4N2, IDTPRPE, IDTC3H8
     . . .
     IF ( NAMEGAS(I) == 'ISN2'   ) IDISN2   = I
     . . .
     IDISN2    = 0

with:

     INTEGER            :: IDISN1,  IDR4N2,   IDCH2O,  IDC2H6,  IDMP
     . . .
     INTEGER            :: IDTPPN,  IDTISN1,  IDTR4N2, IDTPRPE, IDTC3H8     
     . . .
     IF ( NAMEGAS(I) == 'ISN1'   ) IDISN1   = I
     . . .
     IDISN1    = 0

--Bob Y. 16:51, 1 April 2015 (EDT)

Fix for NO3 + ISOP product

This update is slated for inclusion in GEOS-Chem v11-02.

Katie Travis wrote:

The nitrate isoprene nitrate is now called ISN1. It must be added to GEOS-Chem as a transported tracer to undergo dry deposition. We can use the molecular weight given to other isoprene nitrates of 147 g/mol. In addition, the references to ISN2 in drydep_mod.F and planeflight_mod.F (where the ANs are added up) should be changed to ISN1.

NOTE: We will add this fix after FlexChem is integrated into GEOS-Chem. Because FlexChem will completly remove SMVGEAR from GEOS-Chem, we shall add this fix into the KPP master equation file (instead of the SMVGEAR globchem.dat file).

--Melissa Sulprizio (talk) 18:21, 17 June 2015 (UTC)
--Bob Yantosca (talk) 22:17, 8 February 2016 (UTC)

Unresolved issues

Bug in computation of aerodynamic resistance RA

NOTE: This update will be added to GEOS-Chem v11-02.

Brian Boys wrote:

In module GeosCore/drydep_mod.F, I discovered a bug in the subroutine DEPVEL. For the calculation of aerodynamic resistance RA under very stable atmospheric conditions (line 1811), the integration of the stability function phi_h (aka dimensionless vertical temperature gradient) from the roughness length to grid box center doesn't take into account the discontinuity occurring at z/L = 1 where phi_h switches from 1 + 5(z/L) to 5 + z/L (Holtslag, 1993). The result is too high of a value for the integral of phi_h and subsequently RA.

A straightforward solution is to calculate RA under stable conditions (L>0) using an integral form of the stability function, such as equation (12) of Holtslag, 1988. I implemented these changes in drydep_mod.F (from GCv 10-01, public release). To see the effect of the bug, namely a discontinuity in RA, refer to Figure 1, which plots RA as a function of stability (1/L).

Depvel figure 1.png

I've also attached two scatter plots comparing outputted RA (GCv10-01, geosfp, 2x2.5, 2012/11/01) from the base model and an updated simulation employing the bug-fix as described above--please refer to Figure 2 (linear comparison):

Depvel figure 2.png

and Figure 3 (logarithmic comparison):

Depvel figure 3.png

--Bob Yantosca (talk) 18:23, 17 February 2017 (UTC)