Sulfate aerosols

From Geos-chem
Revision as of 20:01, 19 September 2022 by Bmy (Talk | contribs) (Add bug fix in sulfate_mod.F involving reorder IF statements for Fahey and Seinfeld decision algorithm)

Jump to: navigation, search

On this page we provide information about the sulfate aerosol species in GEOS-Chem.

Note that in GEOS-Chem v10-01 and higher versions, all emissions of sulfate aerosols are now handled by the HEMCO emissions component.

Overview

Original formulation

From Park et al [2004]:

The sulfur simulation in GEOS-Chem is based on the Georgia Tech/Goddard Global Ozone Chemistry Aerosol Radiation and Transport (GOCART) model [Chin et al., 2000a], with a number of modifications described below. Our fossil fuel and industrial emission inventory is for 1999-2000 and is obtained by scaling the gridded, seasonally resolved inventory from the Global Emissions Inventory Activity (GEIA) for 1985 [Benkovitz et al., 1996] with updated national emission inventories and fuel use data [Bey et al., 2001a]. The emissions for the United States and Canada are from U.S. EPA [2001], and the emissions for European countries are from European Monitoring and Evaluation Programme (EMEP)/United Nations Economic Commission for Europe (UNECE). Asian sulfur emission in the model is 20 Tg S yr1, which can be compared to year 2000 estimates of 17 Tg S yr1 by Streets et al. [2003] and 25 Tg S yr1 by Intergovernmental Panel on Climate Change (hereinafter IPCC) [2001]. Anthropogenic sulfur is emitted as SO2 except for a small fraction as sulfate (5% in Europe and 3% elsewhere) [Chin et al., 2000a].
Other anthropogenic sources of SO2 in the model include gridded monthly aircraft emissions (0.07 Tg S yr1) taken from Chin et al. [2000a] and biofuel use. We use a global biofuel CO emission inventory with 1° x 1° spatial resolution from Yevich and Logan [2003] and apply an emission factor of 0.0015 mol SO2 per mole CO [Andreae and Merlet, 2001]. Seasonal variations in biofuel emissions are specified from the heating degree days approach [Park et al., 2003].
Natural sources of sulfur in the model include DMS from oceanic phytoplankton and SO2 from volcanoes and biomass burning. The oceanic emission of DMS is calculated as the product of local seawater DMS concentration and sea-to-air transfer velocity. The seawater DMS concentrations are gridded monthly averages from Kettle et al. [1999], and the transfer velocity of DMS is computed using an empirical formula from Liss and Merlivat [1986] as a function of the surface (10 m) wind speed. The GEOS surface winds used here assimilate remote sensing data from the Special Sensor Microwave Imager instrument. Volcanic emissions of SO2 from continuously active volcanoes are included from the database of Andres and Kasgnoc [1998]. Emissions from sporadically erupting volcanoes show large year-to-year variability and are not included in the model. No major volcanic eruptions occurred in 2001. Biomass burning emissions of SO2 are calculated using a gridded monthly biomass burning inventory of CO constrained from satellite observations in 2001 by Duncan et al. [2003] with an emission factor of 0.0026 mol SO2 per mole CO [Andreae and Merlet, 2001].
The gas-phase sulfur oxidation chemistry in the model includes DMS oxidation by OH to form SO2 and MSA, DMS oxidation by nitrate radicals (NO3) to form SO2, and SO2 oxidation by OH to form sulfate. Reaction rates are from DeMore et al. [1997] and the yields of SO2 and MSA from DMS oxidation are from Chatfield and Crutzen [1990]. Aqueous-phase oxidation of SO2 by O3 and H2O2 in clouds to form sulfate is included using kinetic data from Jacob [1986] and assuming a pH of 4.5 for the oxidation by O3. Cloud liquid water content is not available in the GEOS data, and we specify it instead in each cloudy grid box by using a temperature-dependent parameterization [Somerville and Remer, 1984]. The cloud volume fraction in a given grid box is specified as an empirical function of the relative humidity following Sundqvist et al. [1989].
Ammonia emissions in the model are based on annual data for 1990 from the 1° x 1° GEIA inventory of Bouwman et al. [1997]. Source categories in that inventory include domesticated animals, fertilizers, human bodies, industry, fossil fuels, oceans, crops, soils, and wild animals. We view the first five as anthropogenic and the last four as natural. Additional emissions from biomass burning and biofuel use are computed using the global inventories of Duncan et al. [2003] and Yevich and Logan [2003], with an emission factor of 1.3 g NH3 per kilogram dry mass burned [Andreae and Merlet, 2001].
Production of total inorganic nitrate (gas-phase nitric acid and aerosol nitrate) in the model is computed from the ozone-NOx-hydrocarbon chemical mechanism.

Important updates to the sulfate aerosol simulation

Notable additions since Park et al [2004]:

  1. Biomass emissions of SO2 and NH3 are now computed by the GFED inventory.
    • The most recent version (Oct 2015) is GFED4
    • You may still uses the older GFED2 or GFED3 inventories for research purposes.
  2. Incorporation of new Volcanic SO2 emissions from Aerocom
  3. Alkalinity computation for Sea salt aerosols
  4. Updates to regional and global anthropogenic emissions inventories
  5. Get liquid water content and cloud fraction directly from GEOS-5 met fields for SO2 chemistry (since GEOS-Chem v8-03-02)
  6. Other minor changes

Also, the following updates have been added (or are in the process of being added) since GEOS-Chem v9-02:

--Bob Y. (talk) 15:34, 26 October 2015 (UTC)

Cloud water pH for sulfate formation

This update was tested in the 1-month benchmark simulation v9-02p and approved on 13 Sep 2013.

Becky Alexander wrote:

Bulk cloud pH is calculated iteratively using concentrations of sulfate, total nitrate (HNO3 + NO3), total ammonia (NH3 + NH4), SO2, and CO2 = 390 ppmv based on their effective Henry's law constants and the local cloud LWC.
Over the oceans, the influence of cloud droplet heterogeneity in pH on in-cloud sulfate production rates is accounted for using the Yuen et al. (1996) parameterization. Based on isotopic evidence, this parameterization seems to work well over the oceans using sea salt aerosol as the course mode aerosol component, but tends to overestimate in-cloud sulfate production over land.

The reference for this work is:

Alexander, B., D.J. Allman, H.M. Amos, T.D. Fairlie, J. Dachs, D.A. Hegg and R.S. Sletten, Isotopic constraints on sulfate aerosol formation pathways in the marine boundary layer of the subtropical northeast Atlantic Ocean, J. Geophys. Res., 117, D06304, doi:10.1029/2011JD016773, 2012.

--Melissa Sulprizio 11:55, 5 September 2013 (EDT)

Update DMS climatology to Lana

This update was validated with 1-month benchmark simulation v11-01b and 1-year benchmark simulation v11-01b-Run0. This version was approved on 19 Aug 2015.

Monthly average DMS seawater concentrations at 1° x 1° resolution will be implemented in GEOS-Chem v11-01 (via the HEMCO emissions component). These data are described in Lana et al. (2011).

--Melissa Sulprizio (talk) 19:23, 20 July 2015 (UTC)

Sulfur oxidation by reactive halogens

This update was included in v11-02d (approved 12 Feb 2018).

From Chen et al. (2017):

Sulfur and reactive bromine (Bry) play important roles in tropospheric chemistry and the global radiation budget. The oxidation of dissolved SO2 (S(IV)) by HOBr increases sulfate aerosol abundance and may also impact the Bry budget, but is generally not included in global climate and chemistry models. In this study, we implement HOBr + S(IV) reactions into the GEOS-Chem global chemical transport model and evaluate the global impacts on both sulfur and Bry budgets. Modeled HOBr mixing ratios on the order of 0.1–1.0 parts per trillion (ppt) lead to HOBr + S(IV) contributing to 8% of global sulfate production and up to 45% over some tropical ocean regions with high HOBr mixing ratios (0.6–0.9 ppt). Inclusion of HOBr + S(IV) in the model leads to a global Bry decrease of 50%, initiated by the decrease in bromide recycling in cloud droplets. Observations of HOBr are necessary to better understand the role of HOBr + S(IV) in tropospheric sulfur and Bry cycles.

Text S2 in this supporting document describes the parameterization of HOBr + S(IV) reactions in GEOS-Chem.

Reference:

Chen, Q., J. A. Schmidt, V. Shah, L. Jaeglé, T. Sherwen, and B. Alexander, Sulfate production by reactive bromine: Implications for the global sulfur and reactive bromine budgets, Geophys. Res. Lett., 44, 7069–7078, doi:10.1002/2017GL073812, 2017.

Metal catalyzed oxidation of SO2

This update was included in v11-02e (approved 24 Mar 2018).

Becky Alexander wrote:

SO2 is oxidized in clouds by transition metals (Fe and Mn). Natural Fe and Mn atmospheric concentrations are scaled to dust, and anthropogenic are scaled to primary anthropogenic sulfate. It is assumed that 1% of natural Mn and Fe is soluble, for anthropogenic it is 10%. The oxidation state of Fe and Mn depends on sunlight. See Alexander et al. [2009] for more details.
We had discussed having a switch for this in input.geos, so people could easily turn it off if they wanted. The Aerosols WG decided to have it on by default.

Reference:

Alexander, B., Park, R.J., Jacob, D.J., and Gong, S., Transition metal catalyzed oxidation of atmospheric sulfur: Global implications for the sulfur budget, J. Geophys. Res., 114, D02309, 2009.

Viral Shah implemented the metal catalyzed in-cloud SO2 oxidation pathway originally described in Alexander et al. (2009) into GEOS-Chem v10-01.

Viral Shah wrote:

My method largely follows Becky's implementation. The main difference is that instead of using a tracer for primary sulfate to calculate anthropogenic Fe and Mn concentrations, I have added a tracer for anthropogenic Fe (pFe). pFe is emitted along with primary sulfate with an emissions ratio that equals the scaling factor used by Becky to calculate Fe concentrations from primary sulfate. This emission ratio is added as a scaling factor in HEMCO_Config and can be adjusted in the future. For wet and dry deposition, pFe is treated as an aerosol species. Anthropogenic Mn concentrations are calculated by scaling pFe concentrations. Note that Fe and Mn are also present in natural dust, and the GC dust species are used to calculate the natural Fe and Mn concentrations.

--Melissa Sulprizio (talk) 19:36, 1 February 2018 (UTC)

Tagged sulfate and nitrate simulation

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

Becky Alexander wrote:

Sulfate and nitrate each have several different tracers that correspond to their production pathway. This allows for an off-line calculation of their oxygen isotopes (O-17 excess), and allows one to determine the relative importance of each production pathway to the burden of sulfate and nitrate.

--Melissa Sulprizio (talk) 20:50, 15 June 2015 (UTC)

Ocean ammonia emission inventory

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

Fabien Paulot and others have created a new emission inventory for ammonia from oceans. The reference for this work is

Paulot, F., D.J. Jacob, M. Johnson, T.G. Bell, A.R. Baker, W.C. Keene, I.D. Lima, S.C. Doney, and C.A. Stock, Global oceanic emission of ammonia: constraints from seawater and atmospheric observations, Global Biogeochemical Cycles, in press, 2015. [ PDF ]

These emissions will be implemented into GEOS-Chem v11-02 via the HEMCO emissions component.

--Melissa Sulprizio (talk) 17:50, 4 August 2015 (UTC)

Source code and data

In GEOS-Chem v10-01 and newer versions, the sulfate emissions data files are read with the HEMCO emissions component. We have created data files (in COARDS-compliant netCDF format) for use with HEMCO. These data files are contained in the HEMCO data directory tree. For detailed instructions on how to download these data files to your disk server, please see our Downloading the HEMCO data directories wiki post.

--Melissa Sulprizio (talk) 19:17, 20 July 2015 (UTC)

Computing PM2.5 concentrations from GEOS-Chem output

For information on how to compute particulate matter (PM2.5) from GEOS-Chem diagnostic outputs, please see our Particulate matter in GEOS-Chem wiki page.

--Bob Yantosca (talk) 21:13, 10 February 2016 (UTC)

References

  1. Andreae, M. O., and P. Merlet (2001), Emission of trace gases and aerosols from biomass burning, Global Biogeochem. Cycles, 15(4), 95596
  2. Andres, R. J., and A. D. Kasgnoc, A time-averaged inventory of subaerial volcanic sulfur emissions, J. Geophys. Res., 103(D19), 25,251-25,261, 1998.
  3. Benkovitz, C. M., M. T. Scholtz, J. Pacyna, L. Tarrason, J. Dignon, E. C. Voldner, P. A. Spiro, J. A. Logan, and T. E. Graedel, Global gridded inventories of anthropogenic emissions of sulfur and nitrogen, J. Geophys. Res., 101(D22), 29,239-29,253, 1996.
  4. 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
  5. Bouwman, A. F., D. S. Lee, W. A. H. Asman, F. J. Dentener, K. W. VanderHoek, and J. G. J. Olivier, A global high-resolution emission inventory for ammonia, Global Biogeochem. Cycles, 11(4), 561-587, 1997.
  6. Chatfield, R. B., and P. J. Crutzen, Are there interactions of iodine and sulfur species in marine air photochemistry?, J. Geophys. Res., 95(D13), 22,319-22,341, 1990.
  7. DeMore, W. B., S. P. Sander, D. M. Golden, R. F. Hampson, M. J. Kurylo, C. J. Howard, A. R. Ravishankara, C. E. Kolb, and M. J. Molina, Chemical kinetics and photochemical data for use in stratospheric modeling, JPL Publ., 97-4, 1-278., 1997.
  8. Duncan, B. N., R. V. Martin, A. C. Staudt, R. Yevich, and J. A. Logan, Interannual and seasonal variability of biomass burning emissions constrained by satellite observations, J. Geophys. Res., 108(D2), 4100, doi:10.1029/2002JD002378, 2003. PDF
  9. Jacob, D. J., Chemistry of OH in remote clouds and its role in the production of formic acid and peroxymonosulfate,J. Geophys. Res., 91(D9), 9807-9826, 1986.
  10. Kettle, A. J., et al. A global database of sea surface dimethylsulfide (DMS) measurements and a procedure to predict sea surface DMS as a function of latitude, longitude, and month, Global Biogeochem. Cycles, 13(2), 399-444., 1999.
  11. Liss, P. S., and L. Merlivat (1986), Air-sea gas exchange rates: Introduction and synthesis, in The Role of Air-Sea Exchange in Geochemical Cycling, edited by P. Buat-Me´nard, pp. 113-127, D. Reidel, Norwell, Mass, 1986.
  12. Park, R. J., D. J. Jacob, B. D. Field, R. M. Yantosca, and M. Chin, Natural and transboundary pollution influences on sulfate-nitrate-ammonium aerosols in the United States: implications for policy, J. Geophys. Res., 109, D15204, 10.1029/2003JD004473, 2004. PDF
  13. Yevich, R., and J. A. Logan, An assessment of biofuel use and burning of agricultural waste in the developing world, Global Biogeochem. Cycles, 17(4), 1095, doi:10.1029/2002GB001952, 2003. PDF
  14. Somerville, R. C. J., and L. A. Remer, Cloud optical thickness feedbacks in the CO2 climate problem, J. Geophys. Res., 89(D6), 9668-9672, 1984.
  15. Streets, D. G., et al. An inventory of gaseous and primary aerosol emissions in Asia in the year 2000, J. Geophys. Res., 108(D21), 8809, doi:10.1029/2002JD003093, 2003.

--Bob Y. 14:36, 23 February 2010 (EST)

Previous issues that have now been resolved

The following bugs and/or technical issues in the sulfate aerosol module have since been resolved:

Give SO4s and NITs the same molecular weight as SALC

This update was validated with 1-month benchmark simulation v11-01e (approved 04 Jan 2016).

For several past GEOS-Chem versions, SO4s and NITs had the wrong molecular weights listed in the input.geos file.

Bob Yantosca wrote:

In GEOS-Chem v9-02, we changed the molecular weights of the sea salt species SALA and SALC in the input.geos file from 0.036 kg/mol to 0.0314 kg/mol. But somehow this change never got propagated to drydep_mod.F, where we were still using the old molecular weight of 0.036 kg/mol. This is now fixed in the GEOS-Chem v11-01e development code.

One question that I still have is this: What should we use for the MW’s of SO4s and NITs in drydep_mod.F? Both of these species currently use molecular weights of 0.036 kg/mol in drydep_mod.F. (This was the old sea salt molecular weight prior to v9-02.) But the actual molecular weight of SO4s is 0.096 kg/mol and that of NITs is 0.062 kg/mol. So for dry deposition purposes, should we use the MW of sea salt, or should we use the actual MW’s of SO4s and NITs? I looked back throughout the code but could not find any comments documenting this. This behavior has been in the code at least through v7-03-01, and probably earlier.

Becky Alexander replied:

The reason for using sea salt's molecular weight for SO4s and NITs is that these tracers are essentially internally mixed with coarse sea salt aerosol (SALC). As coarse sea salt aerosol likely dominates the mass of these aerosols, it is appropriate to use sea salt's MW.

Another explanation is that since SO4s and NITs are internally mixed with sea salt, they should be treated identically to SALC in the code for all processes.

Therefore, in GEOS-Chem v11-01e, we have changed the molecular weights of SO4s and NITs to 31.4 g/mol in the GEOS-Chem species database. These values will be used everywhere throughout GEOS-Chem.

NOTE: This update caused an error in the production of SO4s and NITS in routine SEASALT_CHEM. This will be fixed in v11-02. Please see this section for more information.

--Bob Yantosca (talk) 16:28, 27 January 2017 (UTC)

Bug fix for sea salt alkalinity in sulfate_mod.F

This fix was validated with the 1-month benchmark simulation v11-01c and approved on 14 Sept 2015.

Johan Schmidt wrote:

There is a bug in the GET_ALK subroutine (now part of sulfate_mod) that causes the computed rate of uptake of SO2 and HNO3 on to sea salt aerosol to be much too slow. GET_ALK only considers the surface area of freshly emitted sea salt (contained in N1 and N2) when calculating the rate, it should consider the surface area of all sea salt aerosols in the grid box.
The bug results in too low sulfate formation, and retained SSA alkalinity in regions that should have acidic SSA.
I would like to propose a fix where the uptake rates (KT1, KT2, KTN1, and KTN2) are calculated using the SSA surface area (from TAREA) and effective radius (from ERADIUS).
I have outlined the bug fix below.
I've added the following variable to GET_ALK
     JLOOP = JLOP(I,J,L)

     SA1= TAREA(JLOOP,4+NDUST) !in cm2/cm3
     SA2= TAREA(JLOOP,5+NDUST) !in cm2/cm3

     R1 = ERADIUS(JLOOP,4+NDUST) !in cm
     R2 = ERADIUS(JLOOP,5+NDUST) !in cm
I've then modified the section that calculates KT1, KT2, KTN1, KTN2 as below:
        ! calculate gas-to-particle rate constant for uptake of 
        ! SO2 onto fine sea-salt aerosols [Jacob, 2000] analytical solution
        CONST1 = 4.D0/(V*GAMMA_SO2)
        !A1     = (RAD1/DG)+CONST1
        !B1     = (RAD2/DG)+CONST1
        !TERM1A = ((B1**2)/2.0d0) - ((A1**2)/2.0d0)
        !TERM2A = 2.D0*CONST1*(B1-A1)
        !TERM3A = (CONST1**2)*LOG(B1/A1)
        !KT1    = 4.D0*PI*N1*(DG**3)*(TERM1A - TERM2A + TERM3A)
        KT1    = SA1/((R1/DG) + CONST1)
        
        !----------------------------------
        ! SO2 uptake onto coarse particles 
        !----------------------------------
        
        ! calculate gas-to-particle rate constant for uptake of 
        ! SO2 onto coarse sea-salt aerosols [Jacob, 2000] analytical solution
        CONST2 = 4.D0/(V*GAMMA_SO2)
        !A2     = (RAD2/DG)+CONST2
        !B2     = (RAD3/DG)+CONST2
        !TERM1B = ((B2**2)/2.0d0) - ((A2**2)/2.0d0)
        !TERM2B = 2.D0*CONST2*(B2-A2)
        !TERM3B = (CONST2**2)*LOG(B2/A2)
        !KT2    = 4.D0*PI*N2*(DG**3)*(TERM1B - TERM2B + TERM3B)
        KT2    = SA2/((R2/DG) + CONST2)
        KT     = KT1 + KT2
        
        !----------------------------------
        ! HNO3 uptake onto fine particles 
        !----------------------------------
        
        ! calculate gas-to-particle rate constant for uptake of 
        ! HNO3 onto fine sea-salt aerosols [Jacob, 2000] analytical solution
        CONST1N = 4.D0/(V*GAMMA_HNO3)
        !A1N     = (RAD1/DG)+CONST1N
        !B1N     = (RAD2/DG)+CONST1N
        !TERM1AN = ((B1N**2)/2.0d0) - ((A1N**2)/2.0d0)
        !TERM2AN = 2.D0*CONST1N*(B1N-A1N)
        !TERM3AN = (CONST1N**2)*LOG(B1N/A1N)
        !KT1N    = 4.D0*PI*N1*(DG**3)*(TERM1AN - TERM2AN + TERM3AN)
        KT1N    = SA1/((R1/DG) + CONST1N)
        
        !----------------------------------
        ! HNO3 uptake onto coarse particles 
        !----------------------------------
        
        ! calculate gas-to-particle rate constant for uptake of 
        ! HNO3 onto coarse sea-salt aerosols [Jacob, 2000] analytical solution
        CONST2N = 4.D0/(V*GAMMA_HNO3)
        !A2N     = (RAD2/DG)+CONST2N
        !B2N     = (RAD3/DG)+CONST2N
        !TERM1BN = ((B2N**2)/2.0d0) - ((A2N**2)/2.0d0)
        !TERM2BN = 2.D0*CONST2N*(B2N-A2N)
        !TERM3BN = (CONST2N**2)*LOG(B2N/A2N)
        !KT2N    = 4.D0*PI*N2*(DG**3)*(TERM1BN - TERM2BN + TERM3BN)
        KT2N    = SA2/((R2/DG) + CONST2N)

--Melissa Sulprizio (talk) 15:22, 3 August 2015 (UTC)

Updated THNO3.geos4.4x5 file

Lyatt Jaegle wrote:

I was trying to run GEOS-Chem (v8-01-03) in the offline aerosol mode with GEOS-4 met fields and ran into a problem: it seems that the offline file GEOS_4x5/sulfate_sim_200508/offline/THNO3.geos4.4x5 is in units of ppbv instead of v/v as expected in routine GET_HNO3_UGM3. This leads to issues in RPMARES which thinks that HNO3 is very large. I checked all the other files, and they are in v/v, as expected.

Bob Yantosca wrote:

Thanks Lyatt. I've downloaded the file and updated the README in the GEOS_4x5/sulfate_sim_200508/offline directory.

--Bob Y. 14:19, 24 April 2009 (EDT)

Fix for mass balance of HNO3 and NIT

NOTE: This fix is was standardized in GEOS-Chem v8-01-02.

Becky Alexander wrote:

We need to make a change in sulfate_mod in order to have mass balance for HNO3 and NIT. Duncan Fairlie noticed the bug. There is a simple change:
In routine SEASALT_CHEM in sulfate_mod.f: In order to have mass balance, you need to change:
   !HNO3 lost [eq/timestep] converted back to [v/v/timestep]
   HNO3_ss = TITR_HNO3 * 0.063 * TCVV(IDTHNO3)/AD(I,J,L)
to:
   !HNO3 lost [eq/timestep] converted back to [v/v/timestep]
   HNO3_ss = HNO3_SSC * 0.063 * TCVV(IDTHNO3)/AD(I,J,L)
In my original code where I added isorropia and the new tracers, NITs and SO4s, the line above:
   !HNO3 lost [eq/timestep] converted back to [v/v/timestep]
   HNO3_ss = TITR_HNO3 * 0.063 * TCVV(IDTHNO3)/AD(I,J,L)
is appropriate as long as you also have PNIT (analogous to PNITs). PNIT is in my original code where I did all my mass balance testing. PNIT got dropped when going to the standard version. I don't recall dropping this, but my guess is that I decided it was redundant to have it when isorropia would just repartition HNO3 and NIT anyway according to thermodynamic equilibrium. But when dropping PNIT, you have to change TITR_HNO3 to HNO3_SSC in the above equation in order to achieve mass balance.

--Bob Y. 16:36, 19 February 2010 (EST)

Fix for sulfate production in HET_DROP_CHEM

This fix was included in v11-02a and approved on 12 May 2017.

Qianjie Chen wrote:

I want to report a bug in the sulfate module of GEOS-Chem.

The bug is in the HET_DROP_CHEM subroutine. It results in too much sulfate production in the Yuen et al. 1996 parameterization (the SR), i.e. too much sulfate produced from ozone over the ocean. Below shows how we can fix it:

(1) The calculation of NDss should multiply by 1.0D-6, to make it with the unit of #/cm3 air.
     NDss = ((3.d0/4.d0) * CNss ) / 
    &       (PI * SS_DEN * RG_S**3.d0 * exp( (9.d0/2.d0) * 
    &       (LOG(SIG_S)) ** 2.d0 ) )*1.d-6
(2) The updraft velocity over the ocean "W" of 500 cm/s is too high. We should get it from the met field:
    W = -OMEGA(I,J,L)/(AIRDEN(L,I,J)*grav)*100

    Then add the following:

    REAL*8, PARAMETER   :: grav = 9.80665d0
    REAL*8, POINTER     :: OMEGA(:,:,:)
    AIRDEN => State_Met%AIRDEN
    OMEGA  => State_Met%OMEGA
    NULLIFY( OMEGA  )
    NULLIFY( AIRDEN )
(3) To be consistent with Yuen et al. (1996), we should compute SR only when air parcel rises, i.e. W>0. Thus,
     IF (W>0d0) THEN
      
        ! additional sulfate production that can be attributed to
        ! ozone [ug/m3/timestep]
        SR = DSVI - B

        ! Convert SR from [ug/m3/timestep] to [v/v/timestep]
        SR = SR *TCVV(IDTSO4) * 1.D-9 / AIRDEN(L,I,J)

        ! Don't allow SR to be negative
        SR = MAX( SR, 0.d0 )

        ! Don't produce more SO4 than SO2 available after AQCHEM_SO2
        SR = MIN( SR, SO2_sr )

     ELSE
       SR = 0d0
     ENDIF
(4) In order to get OMEGA, we make the following changes:
    ## a6_read_mod.f ###

    - get_a6_fields

    Add: (line 608)
        &                 OMEGA      = State_Met%OMEGA,      

    Add: (line 694)
        &                 OMEGA      = State_Met%OMEGA)         

    - read_a6

    change:
         SUBROUTINE READ_A6( NYMD,      NHMS,   
        &                    CLDF,      CLDMAS,    CLDTOPS,   CMFMC, 
        ......
        &                    ZMMD,      ZMMU,      T_FULLGRID)
    into:
         SUBROUTINE READ_A6( NYMD,      NHMS,   
        &                    CLDF,      CLDMAS,    CLDTOPS,   CMFMC, 
        ......
        &                    ZMMD,      ZMMU,   T_FULLGRID,OMEGA)

    Add: (line 985)
         REAL*8,  INTENT(OUT), OPTIONAL :: OMEGA(IIPAR,JJPAR,LLPAR)

    In  "CASE ( 'OMEGA' )"
    change:
                  IF ( CHECK_TIME( XYMD, XHMS, NYMD, NHMS ) ) THEN
                     NFOUND = NFOUND + 1 
    into:
                  IF ( CHECK_TIME( XYMD, XHMS, NYMD, NHMS ) ) THEN
                     IF ( PRESENT( OMEGA ) ) CALL TRANSFER_3D( D, OMEGA )     
                     NFOUND = NFOUND + 1 

    Add: (line 1811)
            ! GEOS-5 OMEGA  
            IF ( PRESENT( OMEGA ) ) THEN
               AD66(:,:,1:LD66,7) = AD66(:,:,1:LD66,7) + OMEGA(:,:,1:LD66)
            ENDIF       


    ### gigc_stat_met_mod.f ###
    Add: (line 170)
        REAL*8,  POINTER :: OMEGA     (:,:,:)  ! Updraft velocity [Pa/s]     

    Add: (line 691)
       ALLOCATE( State_Met%OMEGA     ( IM, JM, LM   ), STAT=RC )
       IF ( RC /= GIGC_SUCCESS ) RETURN
       State_Met%OMEGA     = 0.d0

    Add: (line 971)
       IF ( ASSOCIATED( State_Met%OMEGA      )) DEALLOCATE( State_Met%OMEGA      )     


    ### CMN_DIAG_mod.f ###
    change:
         INTEGER, PARAMETER :: PD66=6
    into:
         INTEGER, PARAMETER :: PD66=7                  


    ### gamap_mod.f ###
    change: (line 6721)
            NTRAC(66) = 6
    into:
            NTRAC(66) = 7     !(qjc, 11/13/15) add OMEGA, vertical velocity

    Add: (line 6755)
                  CASE( 7 )          ! (qjc, 11/13/15) add OMEGA, vertical velocity
                     NAME(T,66) = 'OMEGA'
                     UNIT(T,66) = 'Pa/s'

    ### diag3.f ###
    Add: (line 4130)
                 ! OMEGA    
                  CASE( 7 )
                     SCALEX = SCALE_ND66 
                     UNIT   = 'Pa/s'

--Melissa Sulprizio (talk) 16:08, 2 February 2017 (UTC)

Fix error in production of SO4s, NITs in SEASALT_CHEM routine

This fix was included in v11-02a and approved on 12 May 2017.

Prasad Kasibhatla reported an issue in the production of SO4 and NITs in routine SEASALT_CHEM (in GeosCore/sulfate_mod.F) due to a unit conversion error.

Prasad Kasibhatla wrote:

The problem is that in going from v/v to equivalents and back to v/v in rather than doing it the straightforward way as you say, 2 molecular weights are used each time. The way the conversion is done is: Before calculating the production of SO4s or NITs
    equivalents = v/v * (mol wt of SO2 or HNO3)
                      / (mol wt of air)
                      * mass of air in grid
                      / (mol wt of SO2 or HNO3)
                      * (2 for SO2 or 1 for HNO3 )
and then in calculating the production of SO4s or NITs
     v/v = equivalents * (mol wt of SO4s or NITs)  
                       * (mol wt of air)
                       / (mol wt of SO4s or NITs)
                       / mass of air in grid 
                       / (2 for SO2 or 1 for HNO3)
So the problem is the difference in the terms highlighted in RED, which correspond to the TCVV factors that have been used in GEOS-Chem unit conversions. Btw, I think this affects SO4 too because uptake of SO2 on fine mode seasalt is considered to produces SO4 (and not SO4s). And also, uptake of HNO3 on fine mode seasalt is not accounted for either as a loss of HNO3 or a gain of NIT or NITs.

Bob Yantosca replied:

Older GEOS-Chem versions used the TCVV factor (MW of air / MW of species) in the equations. In GEOS-Chem v11-01, we removed all instances of TCVV and replaced them with direct coding (e.g. TCVV(IDTSO2) became AIRMW / MW_SO2, etc.

Further complicating matters is that the equations in routine SEASALT_CHEM uses hardwired molecular weights (e.g. 0.064 kg/mole, 0.096 kg/mole) as well.

This bug appears to have arisen due to two factors:

  1. Removing TCVV, and
  2. Changing the molecular weight of SO4s and NITs species to 31.4 g/mole (in v11-01e), as per Becky Alexander's recommendation:
The reason to use the molecular weight of sea salt for NITs and SO4s is that these tracers are internally mixed with sea salt aerosol, so their deposition fluxes should be identical to course-mode sea salt aerosol.

To illustrate the issue, we list below the lines in the code where the production of SO4s and NITs are computed (skipping other lines for clarity) in v9-02, v10-01, and v11-01.

   #############################################################
   #### Production of NITs, computed at end of SEASALT_CHEM ####
   #### Other terms have been omitted for clarity           ####
   #############################################################

   In v9-02:

     !SO4F produced converted from [eq/timestep] to [v/v/timestep] 
     PSO4F        = SO4F * 0.096 * TCVV(IDTSO4S) /
    &               AD(I,J,L) / 2.0d0

     Correct, because the MW of SO4s = 96 g/mole.

   In v10-01:

     !SO4F produced converted from [eq/timestep] to [v/v/timestep] 
     PSO4F        = SO4F * 0.096 * TCVV(IDTSO4S) /
    &               AD(I,J,L) / 2.0e+0_fp

    Correct, because the MW of SO4s = 96 g/mole.

   In v11-01:

     !SO4F produced converted from [eq/timestep] to [v/v/timestep] 
     PSO4F        = SO4F * 0.096 * ( AIRMW / MW_SO4S ) /
    &               AD(I,J,L) / 2.0e+0_fp

     Incorrect, because the MW of SO4s = 31.4 g/mole, not 96 g/mole.

   #############################################################
   #### Production of NITs, computed at end of SEASALT_CHEM ####
   #### Other terms have been omitted for clarity           ####
   #############################################################

   In v9-02:

     ! NITS produced converted from [eq/timestep] to [v/v/timestep] 
     PNITs(I,J,L) = HNO3_SSC * 0.063 * TCVV(IDTNITS) / AD(I,J,L)

     Correct, because the MW of NITs = 63 g/mole.

   In v10-01:

     ! NITS produced converted from [eq/timestep] to [v/v/timestep] 
     PNITs(I,J,L) = HNO3_SSC * 0.063 * TCVV(IDTNITS) / AD(I,J,L)
 
     Correct, because the MW of NITs = 63 g/mole.

   In v11-01:

     ! NITS produced converted from [eq/timestep] to [v/v/timestep] 
     PNITs(I,J,L) = HNO3_SSC * 0.063 * ( AIRMW / MW_NITS ) / AD(I,J,L)

     Incorrect, because the MW of NITs = 31.4 g/mole, not 63 g/mole.
Becky Alexander had also written:
However, for this particular unit conversion going from molar equivalents to v/v, the calculation should use the molecular weight of HNO3 for mass balance. I see in the wiki that the unit conversion used to use TCVV(HNO3). In v10 it uses TCVV(NITS). I don't know when this change happened, but it should not have (though the difference is minor, only 1 g/mole). However, for v11, it is definitely a factor of 2 off if you use the molecular weight of seasalt to convert NITS from molar equivalents to v/v. This should be fixed in v11. This unit conversion should use TCVV(HNO3).

Therefore, to fix this bug, one could do the following:

  1. Use 0.096 * ( AIRMW / MW_SO4 ) in the equation for PSO4F
  2. Use 0.063 * ( AIRMW / MW_NIT ) in the equation for PNITs

But an even better solution (and the one we will follow) is to rewrite the various equations to remove all references to hardwired molecular weights, as shown below:

     SO2_eq       = ( ( 2.0_fp * SO2_cd * AD(I,J,L) ) / AIRMW  ) * 1000.0_fp
 
     HNO3_eq      = ( ( HNO3_vv * AD(I,J,L) ) / AIRMW ) * 1000.0_fp

     SO2_ss       = ( SO2_chem * AIRMW / AD(I,J,L) ) / 2000.0_fp
 
     PSO4E        = ( SO4E * AIRMW / AD(I,J,L) ) / 2000.0_fp

     PSO4F        = ( SO4F * AIRMW / AD(I,J,L) ) / 2000.0_fp

     HNO3_ss      = ( HNO3_SSC  * AIRMW / AD(I,J,L) ) / 1000.0_fp
      and   
     HNO3_ss      = ( TITR_HNO3 * AIRMW / AD(I,J,L) ) / 1000.0_fp

     PNITS(I,J,L) = ( HNO3_SSC * AIRMW / AD(I,J,L) ) / 1000.0_fp
Furthermore, we can take this opportunity to remove the hardwired molecular weights for SO4 and HNO3 by rewriting these equations as well:
     C_FLUX_A     = ( 2.0_fp * F_SO2_A * AD(I,J,L) / AIRMW ) * 1000.0_fp

     C_FLUX_C     = ( 2.0_fp * F_SO2_C * AD(I,J,L) / AIRMW ) * 1000.0_fp

     N_FLUX_A     = ( F_HNO3_A * AD(I,J,L) / AIRMW ) * 1000.0_fp

     N_FLUX_C     = ( F_HNO3_C * AD(I,J,L) / AIRMW ) * 1000.0_fp

--Bob Yantosca (talk) 19:33, 10 February 2017 (UTC)

Fix bug in CHEM_NIT routine

This fix was included in v11-02a and approved on 12 May 2017.

Prasad Kasibhatla wrote:

I am wondering if there is a small bug in CHEM_NIT. I think the nitrate uptake on fine sea salt (the line displayed in GREEN)
          ! Fraction of box (I,J,L) underneath the PBL top [unitless]
          F_UNDER_TOP = GET_FRAC_UNDER_PBLTOP( I, J, L )     
   
          ! Only apply drydep to boxes w/in the PBL
          IF ( F_UNDER_TOP > 0e+0_fp ) THEN  
 
             !===========================================================
             ! NITs chemistry
             !=========================================================== 
 
             ! NIT prod from HNO3 uptake on fine sea-salt [v/v/timestep]
             NITs = NIT0s + PNITs(I,J,L) 
 
             ! Store final concentration in Spc [v/v]
             Spc(I,J,L,id_NITs) = NITs
             
          ENDIF

should be done throughout the troposphere and not just in the PBL. The PBL check seems to be a legacy from when dry dep was done in CHEM_NIT. I am guessing that even if this is a bug, it is a very small bug because PNITs should be quite small outside the marine PBL.

--Bob Yantosca (talk) 19:43, 25 January 2017 (UTC)

Fix bugs in sulfate chemistry routines

These fixes were included in v11-02a and approved on 12 May 2017.

Viral Shah wrote:

I am running GEOS-Chem v10-01 in the full chemistry configuration with GEOS-FP met fields, and focusing on the sulfate aerosol simulation in Feb and March over the U.S. While diagnosing the cause for low modeled sulfate concentrations compared to surface observations, I found three bugs in sulfate_mod.F. These are present in the public release of GEOS-Chem v11-01 too, and the line numbers below refer to v11-01.

1. HNO3 effective Henry's law coefficient

Module: sulfate_mod.F
Function: kHNO3

Line Original code Recommended fix
5087 Hhno3eff = 3.2e6/HPLUS Hhno3eff = Hhno3_T*(1.0e+0_fp+(Kn1_T/HPLUS))

The original expression is valid for 298K (Seinfeld and Pandis 2006, pp:299-301), and Kn1 has a strong temperature dependence. The recommended fix follows Eq. 7.59 of Seinfeld and Pandis (2006, pp 301). This error was resulting in a low cloud water pH (< 3) over the Arctic and other NH mid-latitude cold regions in Feb-Mar, and affecting sulfate production rates.

2. ∆H/R parameter

Module: sulfate_mod.F
Functions: kCO21, kCO22, kSO21, kSO22, kNH3

Function Line Original code Recommended fix
kCO21 4809 REAL(fp), PARAMETER :: Dhco2 = -4.85/0.04 REAL(fp), PARAMETER :: Dhco2 = 4.85e+0_fp/1.986e-3_fp
kCO22 4875 REAL(fp), PARAMETER :: Dhco2 = -4.85/0.04 REAL(fp), PARAMETER :: Dhco2 = 4.85e+0_fp/1.986e-3_fp
kSO21 4940 REAL(fp), PARAMETER :: Dhso2 = -6.25/0.04 REAL(fp), PARAMETER :: Dhso2 = 6.25+0_fp/1.986e-3_fp
kSO22 5005 REAL(fp), PARAMETER :: Dhso2 = -6.25/0.04 REAL(fp), PARAMETER :: Dhso2 = 6.25+0_fp/1.986e-3_fp
kNH3 5190 REAL(fp), PARAMETER :: Dhnh3 = -8.17/0.04 REAL(fp), PARAMETER :: Dhnh3 = 8.17+0_fp/1.986e-3_fp

The above values are from Table 7.3 of Seinfeld and Pandis (2006, pp 289), but the values here should be positive for consistency with the way they are used in the respective functions. Secondly, the value of R, in units of kcal mol-1 K-1, is 1.986x10^-3, not 0.04.

3. LWC units

Module: sulfate_mod.F
Subroutine: CHEM_SO2

Line Original code Recommended fix
3244 ELSEIF( LWC < 0.1e+6_fp ) THEN
!10^6 coversion from m3/m3 --> g/m3
ELSEIF( LWC < 0.1e-6_fp ) THEN
!10^-6 coversion from g/m3 --> m3/m3
3250 IF ( LWC >= 0.3e+6_fp .and. IF ( LWC >= 0.3e-6_fp
3255 & LWC >= 0.5e+6_fp .and. & LWC >= 0.5e-6_fp .and.
3258 IF ( LWC >= 0.1e+6_fp .and. IF ( LWC >= 0.1e-6_fp .and.
3261 ELSEIF( LWC >= 0.1e+6_fp ) THEN ELSEIF( LWC >= 0.1e-6_fp ) THEN
3309 ELSEIF ( LWC >= 0.3d6 ) THEN ELSEIF ( LWC >= 0.3e-6_fp ) THEN
3345 ELSEIF ( LWC >= 0.5e+6_fp ) THEN ELSEIF ( LWC >= 0.5e-6_fp ) THEN
There was an error in unit conversion. We are comparing LWC (m3/m3) with X (g/m3) or X*10^-6 (m3/m3).

The impacts of these updates are summarized below.

Viral Shah wrote:

The effects of my proposed fixes are a bit more pronounced during the NH winter. The plots show significant increases in sulfate concentrations over large areas of the NH. The SO4 increases are accompanied by decreases in NIT in some areas. The increase in SO4 in the cold areas of NH is mostly because of the correction to HNO3 effective Henry's law coefficient that leads to higher cloud water pH (between 4 and 5) and higher production by in-cloud oxidation by O3. The lower SO4 concentrations in the SH are probably due to the fix to ∆H/R parameter. Also note that these plots do not include the effects of fix #3 (LWC units), but I expect those effects to be relatively small.

--Melissa Sulprizio (talk) 17:59, 22 March 2017 (UTC)

Unresolved issues

None at this time.