Development of the adjoint of GEOS-Chem

21 Pages·2007·7.65 MB·English

Atmos. Chem. Phys.,7,2413–2433,2007
Atmospheric
www.atmos-chem-phys.net/7/2413/2007/
Chemistry
©Author(s)2007. Thisworkislicensed
underaCreativeCommonsLicense. and Physics
Development of the adjoint of GEOS-Chem
D.K.Henze,A.Hakami,andJ.H.Seinfeld
CaliforniaInstituteofTechnology,Pasadena,CA,USA
Received: 4October2006–PublishedinAtmos. Chem. Phys.Discuss.: 19October2006
Revised: 6February2007–Accepted: 17April2007–Published: 11May2007
Abstract. We present the adjoint of the global chemical the substantial uncertainty that remains in many aspects of
transportmodelGEOS-Chem,focusingonthechemicaland detailed aerosol simulations, it is critical to further exam-
thermodynamic relationships between sulfate – ammonium inehowthenumerousparametersinsuchmodelssteertheir
– nitrate aerosols and their gas-phase precursors. The ad- predictions,especiallyestimatesofemissionsinventoriesfor
joint model is constructed from a combination of manually aerosolsandtheirprecursors. Thecomplexityofthethermo-
andautomaticallyderiveddiscreteadjointalgorithmsandnu- dynamicandphotochemicalprocessesthatgovernsecondary
merical solutions to continuous adjoint equations. Explicit formationofaerosolsprecludessimpleassessmentofthede-
inclusion of the processes that govern secondary formation pendenceofmodelpredictionsonsuchparameters. Working
of inorganic aerosol is shown to afford efficient calculation toarriveatCTMsthatmorereliablyreproduceobservations,
ofmodelsensitivitiessuchasthedependenceofsulfateand adjoint modeling is often employed as a method for deter-
nitrate aerosol concentrations on emissions of SO , NO , mining the sensitivity of model predictions to input param-
x x
and NH . The accuracy of the adjoint model is extensively eters and for optimizing these parameters to enforce agree-
3
verified by comparing adjoint to finite difference sensitivi- ment between the model predictions and an observational
ties,whichareshowntoagreewithinacceptabletolerances. dataset.
We explore the robustness of these results, noting how dis- Severalinversemodelingstudieshaveanalyzedsourcesof
continuities in the advection routine hinder, but do not en- aerosols and aerosol precursors on regional scales. As of
tirelypreclude,theuseofsuchcomparisonsforvalidationof yet,moststudieshavebeenfairlycoarse,limitedtooptimiza-
theadjointmodel. Thepotentialforinversemodelingusing tionofafewscalingfactorsforemissionsinventoriesspan-
the adjoint of GEOS-Chem is assessed in a data assimila- ning large domains. Park et al. (2003) used multiple linear
tionframeworkusingsimulatedobservations,demonstrating regression to estimate annual mean sources of seven types
thefeasibilityofexploitinggas-andaerosol-phasemeasure- of primary carbonaceous aerosol over the United States.
mentsforoptimizingemissioninventoriesofaerosolprecur- A Kalman filter approach was used to estimate improved
sors. monthly emissions scaling factors for NH emissions over
3
theUnitedStatesusingobservationsofammoniumwetdepo-
sitioninworksbyGillilandandAbbitt(2001)andGilliland
etal.(2003,2006).Mendoza-DominguezandRussell(2000,
1 Introduction
2001) optimized domain-wide emissions scaling factors for
eight species over the eastern Unites States using observa-
Chemical transport models (CTMs) enhance our ability to
tions of gas-phase inorganic and organic species and speci-
understand the chemical state of the atmosphere and allow
ated fine particles. Source apportionment models have also
detailedanalysisofissuesrangingfromintercontinentalpol-
beenrefinedusinginversemodeling(Knippingetal., 2006;
lution transport to the coupling of anthropogenic processes,
Schichteletal.,2006).
regional pollution and climate change. Of particular inter-
Data from satellite observations offer tremendous poten-
est in these realms is explicit consideration of the role of
tial for inverse modeling of aerosols (Collins et al., 2001;
aerosols,theimportanceofwhichiswelldocumented.Given
Kahnetal.,2004). Inordertobestexploitthese, andother,
Correspondenceto: D.K.Henze large data sets, it is desired to extend inverse analysis of
([email protected]) aerosol models to global scales and to finer decomposition
PublishedbyCopernicusGmbHonbehalfoftheEuropeanGeosciencesUnion.
2414 D.K.Henzeetal.: AdjointofGEOS-Chem
oftheemissionsdomains. Suchgoalsrequireconsideration 2 Forwardandinversemodels
ofinversemodelingmethodsdesignedforlargesetsofvari-
ableparameters. Theadjointmethodisknowntobeaneffi- The GEOS-Chem model is used to simulate global aerosol
cientmeansofcalculatingmodelsensitivitiesthataffordex- distributions(version6.02.05withahorizontalresolutionof
amination of numerous parameters, where these values can 4◦×5◦ and 30 layers up to 0.01hPa, GEOS-3 meteorologi-
subsequently be used in tandem with an observational data calfields). Thisversionofthemodelincludesdetailedgas-
setfordataassimilation. Firstappearinginthefieldofatmo- phasechemistrycoupledwithheterogeneousreactions,inor-
spheric science in the early 1970s (Marchuk, 1974; Lamb ganic aerosol thermodynamics, and oxidative aging of car-
et al., 1975), the method later came to be applied exten- bonaceousaerosols(Parketal.,2004). Afewofthespecific
sively in meteorology, e.g., Talagrand and Courtier (1987); equationsforvariousmodelprocessesaregiveninSect.3.3,
Errico and Vukicevic (1992). In the last decade, the ad- along with their corresponding adjoints. We note here that
joint approach has expanded to include ever more detailed gaseous SO and primary sulfate are co-emitted in GEOS-
2
CTMs, beginning with the abbreviated Lagrangian strato- Chemusingasingleemissionsinventory,referredtoasSO ,
x
sphericmodelofFisherandLary(1995)andtheLagrangian which is partitioned between the two species on a regional
tropospheric model of Elbern et al. (1997). Vukic´evic´ and basis, with sulfate comprising 5% of SO emissions in Eu-
x
Hess(2000)usedtheadjointmethodtoperformasensitivity rope,1.7%inNorthAmerica,and3%elsewhere(Chinetal.,
studyofaninertgas-phasetraceroverthePacific,whileEl- 2000).
bernandSchmidt(1999)presentedthefirstadjointofa3-D Thestandardmodelhasbeenmodifiedtofacilitatethespe-
EulerianCTMtoincludechemistry.Theseinitialworkshave cific inverse modeling goals of the present study. We ne-
beenfollowedmorerecentlybysimilardevelopmentandap- glect stratospheric chemistry, which over the course of the
plication of adjoint models of several CTMs: CHIMERE shortsimulationsconsideredhereshouldnothaveasubstan-
(Vautardetal.,2000;Menutetal.,2000;SchmidtandMar- tial impact. The standard GEOS-Chem tropospheric chem-
tin,2003),IMAGES(MullerandStavrakou,2005;Stavrakou ical mechanism comprises 87 species and 307 reactions in-
andMuller,2006),Polair(MalletandSportisse,2004,2006), tegrated using the SMVGEARII solver of Jacobson (1995).
TM4(Meirinketal.,2006),theCaliforniaInstituteofTech- We retain this standard chemical mechanism; however, we
nologyurban-scalemodel(Martienetal.,2006;Martienand implement a different numerical solver. The details of this
Harley, 2006), and DRAIS (Nester and Panitz, 2006). The aregiveninAppendixA.Tosummarize,weimplementa3rd
adjointoftheregionalmodelSTEMalsohasbeendeveloped order Rosenbrock solver that not only facilitates construc-
(Sandu et al., 2005a) and deployed (Hakami et al., 2005, tion of the adjoint model, but also improves forward model
2006;Chaietal.,2006). efficiency. We also consider using offline concentrations of
Ofalltheprevious3-Dadjointmodelingstudies,nonein- sulfateaerosolforcalculationofphotolysisratesandhetero-
cludesdetailedtreatmentofaerosols,likelyowingtothedif- geneousreactionprobabilities,seeSect.3.5.
ficult prospect of deriving the adjoint of the model routines
dealingwithaerosolthermodynamics. ThestudyofHakami 2.1 Inversemodeling
etal.(2005)dealsonlywithinertcarbonaceousaerosols,and
the work of Dubovik et al. (2004), though global in scale, An adjoint model is used to calculate the gradient of a cost
does not include full chemistry or aerosol thermodynam- function, J, with respect to a set of model parameters, p,
ics. Detailed adjoint modeling of aerosols began with the ∇pJ. Fordataassimilationapplications,thecostfunctionis
theoretical investigations of Henze et al. (2004) and Sandu definedtobe
et al. (2005b). However, these are preliminary studies per-
formedonidealizedboxmodelsystems. Inthecurrentwork J=
we present the first adjoint of a global CTM that includes 1X(c−c )TS−1(c−c )+1γ (p−p )TS−1(p−p )(1)
dynamics,fulltroposphericchemistry,heterogeneouschem- 2 obs obs obs 2 r a p a
c∈(cid:127)
istry, andaerosolthermodynamics. Wedemonstratethepo-
tentialvalueofthistoolforquantifyingandconstrainingfac- where c is the vector of species concentrations mapped to
tors that govern global secondary inorganic aerosol forma- the observation space, cobs is the vector of species observa-
tion. In addition, we note the general usefulness of the ad- tions,Sobs istheobservationerrorcovariancematrix,p isa
joint model of GEOS-Chem for a wide variety of applica- vectorofactivemodelparametersthroughoutthemodeldo-
tions,suchasconstrainingCOemissionsusingsatellitedata main,pa istheinitialestimateoftheseparameters,Sp isthe
(Kopaczetal.,20071). errorcovarianceestimateoftheseparameters,γr isaregular-
ization parameter, and (cid:127) is the domain (in time and space)
1Kopacz,M.,Jacob,D.,Henze,D.K.,Heald,C.L.,Streets,D. overwhichobservationsandmodelpredictionsareavailable.
G.,andZhang,Q.:AcomparisonofanalyticalandadjointBayesian Wewillsometimesusethenotationcandptorepresentsin-
inversionmethodsforconstrainingAsiansourcesofCOusingsatel- gle elements of the vectors c and p. Using the variational
lite(MOPITT)measurementsofCOcolumns,submitted,2007. approach, the gradient ∇pJ is supplied to an optimization
Atmos. Chem. Phys.,7,2413–2433,2007 www.atmos-chem-phys.net/7/2413/2007/
D.K.Henzeetal.: AdjointofGEOS-Chem 2415
routineandtheminimumofthecostfunctionissoughtiter- tionofthediscreteforwardmodelwhichadvancesthemodel
atively. At each iteration, improved estimates of the model statevectorfromstepntostepn+1.
parametersareimplementedandtheforwardmodelsolution Forsimplicity,weconsideracostfunctionevaluatedonly
is recalculated. In this study, the magnitude of each vari- at the final time step N with no penalty term. We wish to
ableparameterisadjustedusingascalingfactor,σ,suchthat calculatethegradientofthecostfunctionwithrespecttothe
p=σp . We use the L-BFGS-B optimization routine (Byrd modelstatevectoratanystepinthemodel,
a
etal.,1995;Zhuetal.,1994), whichaffordsboundedmini-
∂J(cN)
miAzalttieornn,ateinvseulyri,nfgorpsoesnistiivtievivtayluanesalfyosrist,htehseccaolisntgfufnaccttioorns.can ∇cnJ = ∂cn (4)
bedefinedassimplyasetofmodelpredictions, WedefinethelocalJacobianaroundanygivenstepas
X
J = g(c) (2) ∂cn+1 ∂F(cn)
= =Fn (5)
g∈(cid:127)s ∂cn ∂cn c
where (cid:127)s is the set of times at which the cost function is Using the chain rule, we can expand the right hand side of
evaluated. Thedesiredgradientvaluesarethesensitivitiesof Eq.(4)toexplicitlyshowthecalculationofcN fromcn,
thissetofmodelpredictionstothemodelparameters.
∂J(cN)
2.2 Adjointmodeling ∇cnJ =(Fcn)T(Fcn+1)T···(FcN−1)T ∂cN (6)
Evaluatingtheaboveequationfromlefttorightcorresponds
Equationsforcalculatingthedesiredgradientsusingthead-
to a forward sensitivity calculation, while evaluating from
joint method can be derived from the equations governing
righttoleftcorrespondstoanadjointcalculation.WhenK is
the forward model or from the forward model code. The
largerthanthedimensionofJ,whichinthiscaseisascalar,
prior approach leads to the continuous adjoint, while the
the adjoint calculation is much more efficient (Giering and
latter leads to the discrete adjoint (Giles and Pierce, 2000).
Kaminski,1998).
The continuous adjoint equations for CTMs have been de-
Fortheadjointcalculation,wedefinetheadjointstatevari-
rived previously, using methods based upon the Lagrange
ableλn,
duality condition (Vukic´evic´ and Hess, 2000; Pudykiewicz, c
1998; Schmidt and Martin, 2003) or Lagrange multipliers ∂J(cN)
λn = . (7)
(Elbernetal.,1997). Continuousadjointgradientsmaydif- c ∂cn
ferfromtheactualnumericalgradientsofJ,andcontinuous
Thiscanalsobeexpanded,
adjointequations(andrequisiteboundary/initialconditions)
forsomesystemsarenotalwaysreadilyderivable;however, " #T
∂cn+1 ∂J(cN)
solutionstocontinuousadjointequationscanbemoreuseful λn = (8)
c ∂cn ∂cn+1
forinterpretingthesignificanceoftheadjointvalues. Many
previous studies have also described the derivation of dis- ∂J(cN)
crete adjoints of such systems (Sandu et al., 2005a; Muller =(Fcn)T ∂cn+1 . (9)
and Stavrakou, 2005). An advantage of the discrete adjoint
The equation above suggests how to solve for the adjoint
model is that the resulting gradients of the numerical cost
variableiteratively. Initializingtheadjointvariableatthefi-
functionareexact,evenfornonlinearoriterativealgorithms,
naltimestep
makingthemeasiertovalidate. Furthermore,portionsofthe
discreteadjointcodecanoftenbegenerateddirectlyfromthe ∂J(cN)
forwardcodewiththeaidofautomaticdifferentiationtools. λN = (10)
c ∂cN
Here we present a brief description of the discrete adjoint
methodforthesakeofdefiningaself-consistentsetofnota- wesolvethefollowingequationiterativelyfromn=N,...,1,
tionforthisparticularpaper;wereferthereadertothecited λn−1 =(Fn)Tλn (11)
works for further derivations and discussions of continuous c c c
anddiscreteadjoints. The value of λ0 is then the sensitivity of the cost function
c
TheGEOS-Chemmodelcanbeviewedasanumericalop- withrespecttothemodelinitialconditions,
erator,F,actingonastatevector,c
λ0 =∇ J (12)
c c0
cn+1 =F(cn) (3)
The scheme above shows why calculating the adjoint vari-
where c is the vector of all K tracer concentrations, able is often referred to as “reverse integration” of the for-
cn=[cn,...,cn,...,cn]T atstepn. Inpractice,F comprises wardmodel,aswestepfromthefinaltimetotheinitialtime.
1 k K
manyindividualoperatorsrepresentingvariousphysicalpro- Thisshouldnotbeconfusedwithsimplyintegratingthefor-
cesses. ForthemomentwewillsimplyletF representapor- wardmodelequationsbackwardsintime.
www.atmos-chem-phys.net/7/2413/2007/ Atmos. Chem. Phys.,7,2413–2433,2007
2416 D.K.Henzeetal.: AdjointofGEOS-Chem
In order to calculate the sensitivity of J with respect to model domain, a much better approach to revealing poten-
othermodelparameters, suchasemissions, similaranalysis tial errors than performing validation checks in only a few
(see,forexample,Sanduetal.,2003)showsthatthegradient locations. Furthermore, as GEOS-Chem has many routines
ofthecostfunctionwithrespecttotheseparameters, commontoothermodels,itbehoovesustoconsiderthead-
jointoftheseroutinesseparately.
λ0p =∇pJ (13) Forward model sensitivities, 3, are calculated using the
finite difference (brute force) method. For component-wise
canbefoundbyiterativelysolvingthefollowingequation,
tests of nonlinear routines, 3 is calculated using the two-
λn−1 =(Fn)Tλn+λn (14) sidedformula,
p p c p
J(σ +δσ)−J(σ −δσ)
wherethesubscriptscandpindicatesensitivitywithrespect 3= (17)
2δσ
tocandp,respectively,and
whilefortestingthefullmodel, themoreapproximateone-
∂Fn
Fn = (15) sidedfinitedifferenceequation,
p ∂p
J(σ +δσ)−J(σ)
Whenapenaltytermisincludedinthecostfunction,thegra- 3= (18)
δσ
dientbecomes
isusedinordertominimizethenumberofrequiredforward
∇ J =λ0 +γ S−1(p−p ) (16) model function evaluations. The latter method is also ade-
p p r p a
quate for testing linear components of the model. We use
δσ=0.1–0.01formosttests,whichexperienceshowedtobe
3 Constructing and validating the adjoint of GEOS- an optimal balance between truncation and roundoff error.
Chem Formostofthesevalidationtests,itsufficestouseasimpli-
fiedcostfunctionthatdoesnotdependonanyobservational
HerewepresentthederivationoftheadjointofGEOS-Chem. data set, as in Eq. (2), defining g to be a predicted tracer
Whiletheadjointoftheadvectionschemeisbaseduponthe mass,eithergas-oraerosol-phase,inasinglegridcell,orthe
continuous approach, the remainder of the adjoint model is totalmassburdenoveralargerspatialdomain.
baseduponthediscreteformulation,usingautomaticdiffer-
entiation tools for assistance. We use the Tangent and Ad- 3.1 Aerosolthermodynamics
jointModelCompiler(TAMC,GieringandKaminski,1998),
afreewaremultipurposeprogram,andtheKineticPrePoces- The equilibrium thermodynamic model MARS-A
sor(KPP,Sanduetal.,2003;Damianetal.,2002;Daescuet (Binkowski and Roselle, 2003) is used to calculate the
al., 2003), a public domain numerical library for construct- partitioning of total ammonia and nitric acid between
ing the adjoint of chemical mechanisms. Always some, if aerosol and gas phases. While it is a relatively simple
notsignificant,manualmanipulationofthecodeisrequired treatement compared to others such as SCAPE (Kim et al.,
tousesuchtools. Weoftencombineautomaticallygenerated 1993) or ISORROPIA (Nenes et al., 1998), the MARS-A
adjoint code with manually derived discrete adjoint code to model is still fairly complex. It uses an iterative algorithm
improveefficiencyandtransparencyoftheadjointmodel. tofindequilibriumconcentrations, consideringtwoprimary
Validationoftheadjointmodelisanimportantpartofin- regimes defined by the ionic ratio of ammonium to sulfate
troducinganadjointmodelofthissizeandcomplexity. Dis- and several sub-regimes defined by conditions such as
crete portions of the adjoint code have the advantage of be- relativehumidity.
ing easily validated via comparison of adjoint gradients to Several factors have historically prevented rigorous treat-
forwardmodelsensitivitiescalculatedusingthefinitediffer- ment of aerosol thermodynamics from inclusion in adjoint
enceapproximation. Thehybridapproachadoptedhere(dis- modelingstudiesofCTMs,orevenadjointstudiesofaerosol
creteandcontinuous)requiresdetailedinspectionofthead- dynamics(Henzeetal.,2004;Sanduetal.,2005b). Division
joint gradients on a component-wise basis as discrepancies of the possible thermodynamic states into distinct regimes
owing to the continuous portion are anticipated to obscure causes many discontinuities in the derivatives, precluding
suchcomparisonsforthemodelasawhole. Additionalmo- easy derivation of continuous adjoint equations and raising
tivations exit for checking the gradients of subprocesses in doubts to the value of such sensitivities. Furthermore, sev-
themodelseparatelyandcollectively. ForlargeCTMs,itis eral coding tactics often employed in these types of models
notfeasibletocompareadjointandfinitedifferencegradients renderthemintractablefordirecttreatmentusingautomatic
foreachcontrolparameter,asthefinitedifferencecalculation differentiationtools.
requiresanadditionalforwardmodelevaluationperparame- WedeveloptheadjointofMARS-Ainpieces, separating
ter. However,component-wiseanalysisaffordssimultaneous the model into several subprograms, the adjoints of which
examinationoflargenumbersofsensitivitiesthroughoutthe arethencreatedusingTAMC.Trackingvariablesareadded
Atmos. Chem. Phys.,7,2413–2433,2007 www.atmos-chem-phys.net/7/2413/2007/
D.K.Henzeetal.: AdjointofGEOS-Chem 2417
x 104
s 5
e
sitiviti 4 R2 =0.999
n
e
S 3 m =1.00
e
c
en 2
er
Diff 1
e
Finit 00 1 2 3 4 5
-106 -105 -104 -103 -102 -101 101 102 103 104 105 106 x 104
Adjoint Sensitivities [kg / grid cell] Adjoint Sensitivities [kg / grid cell]
s
vitie 3000 R2 =0.999
siti 2000
n
e m =1.00
e S 1000
c
n
e 0
er
Diff-1000
e
nit-2000
Fi -2000-1000 0 1000 2000 3000
-106 -105 -104 -103 -102 -101 101 102 103 104 105 106
Adjoint Sensitivities [kg / grid cell] Adjoint Sensitivities [kg / grid cell]
Fig.1. Thermodynamicadjointvalidation.Intheleftcolumnaretheadjointsensitivitiesofnitrateaerosolmassatthesurfacewithrespectto
anthropogenicNH3andSOxemissionsscalingfactors.Intherightcolumnaretheadjointgradientscomparedtofinitedifferencegradients.
Thecostfunctionisevaluatedonceattheendofaweek-longsimulationthatincludesonlyaerosolthermodynamicsandemissionsofSOx
andNH3.
totheforwardmodelroutinetoindicatewhichofthesesub- adjoint of the tropospheric chemistry solver, which calcu-
routinestocallduringtheadjointcalculation. Initialunequi- lates gradients with respect to the initial species concentra-
libratedconcentrationsatthebeginningofeachexternaltime tions. We are also interested in the gradient with respect
steparesavedincheckpointfilesduringtheforwardcalcula- to the emission rates for those species whose emissions are
tion. Intermediatevaluesarerecalculatedfromtheseduring incorporated into the chemical mechanism itself, such as
theadjointintegration. Thistypeoftwo-levelcheckpointing NO , (as opposed to those that are simply injected into the
x
strategyhasbeenshowntooptimallybalancestorage,mem- modelgridcellsatintermediatetimes,suchasSO ).Thead-
x
ory and CPU requirements (Griewank and Walther, 2000; ditional equations for calculating discrete adjoint gradients
Sanduetal.,2005a). with respect to reaction rate constants are derived in Ap-
The accuracy of the resulting adjoint code is tested by pendix B. Though these equations have not been presented
comparingadjointgradientstofinitedifferencegradientscal- previously, KPP does provide the necessary subroutines for
culatedusingEq.(17)withδσ=0.1. Thesecomparisonscan solvingthem.
be made directly throughout the entire model domain by
Toassesstheaccuracyoftheadjointsofthechemistryrou-
turning off all transport processes. Figure 1 shows compar-
tine,wecalculatethesensitivityofthespeciesconcentrations
isonsforthesensitivityofsurfacelevelnitrateaerosolmass
at the end of a single chemistry time step (1h) with respect
withrespecttoscalingfactorsforemissionsofsurfacelevel
to the emissions of NO (emitted as NO) in a box model
x
anthropogenic SO and NH after a week-long simulation.
x 3 test. Forthistest,thechemicalenvironmentisthatofapol-
The gradients agree quite well, confirming the accuracy of
luted, urban grid cell in the afternoon. Figure 2 shows the
the thermodynamic adjoint code. Discussion of values of
ratio λ /3 for three separate cases. Using a two-
modelsensitivitiesisgiveninSect.4. sidedfiEnNitOexdiffEeNreOnxcecalculation(Eq.17)withδσ =0.1
ENOx
leads to agreement within a few percent. The dependence
3.2 Chemistry oftheinternaltimesteponspeciesconcentrationsisafeed-
backnotaccountedforintheadjointalgorithm; hence, also
KPP(v2.2)(Sanduetal.,2003;Damianetal.,2002;Daescu holdingtheinternaltimestepfixedat60sresultsinratiosof
et al., 2003) is used to automatically generate code for the nearly1.000forallspecies. Forcomparison,theratioswhen
www.atmos-chem-phys.net/7/2413/2007/ Atmos. Chem. Phys.,7,2413–2433,2007
2418 D.K.Henzeetal.: AdjointofGEOS-Chem
1.1 ThecodegeneratedbyKPPallowscomputationofeither
one−sided L the continuous or discrete adjoints of the chemical mecha-
1.08 two−sided L
nism. The continuous adjoint equation can be solved faster
two−sided L , constant step size
1.06 than the discrete adjoint equation at a given tolerance level,
1.04 as calculation of the latter requires recalculation of inter-
Ox mediate values from the forward integration and computa-
N 1.02
E tion of the Hessian during the adjoint integration, see Ap-
L / Ox 1 pendix B. At tight tolerance levels (i.e. very small internal
N 0.98 time steps), the results of these methods should converge.
E
l However, for tolerance levels appropriate for global model-
0.96
ing, the continuous adjoint is only approximate, as λ+δλ,
0.94 where||δλ||<C·Tol. Giventhatthecomputationalexpense
0.92 of the Rosenbrock solver increases substantially for tighter
tolerancelevels(seeAppendixA),itismoreefficienttouse
0.9
10 20 30 40 50 60 70 80 the discrete adjoint, even though this requires an additional
Species Index
forward integration. This is in contrast to the approach of
Errera and Fonteyn (2001), who chose to approximate the
Fig. 2. Chemistry adjoint validation. The ratios of the adjoint
necessaryintermediatevaluesbylinearlyinterpolatingfrom
tofinitedifferencesensitivitiesofeachspecieswithrespecttoNOx
emissionsarecalculatedfora1hboxmodelsimulation.Resultsare values stored at each external time step, an approach likely
shownforaone-sidedfinitedifferencecalculation,δσ =0.1(blue more appropriate for their stratospheric chemistry applica-
(cid:5)’s),atwo-sidedfinitedifferencecalculation(i.e.averageofδσ=0.1 tion.
and−0.1,redx’s)andatwo-sidedfinitedifferencecalculationwith GEOS-Chemaccountsfortheeffectofaerosolconcentra-
afixedinternaltimestepof60s(greeno’s). tions on the radiation available for photolysis reactions and
ontheavailablesurfaceareafortheheterogeneousreactions
Eq.(18)isusedfor3ENOx arealsoshown,whichcandiffer included in the main chemical mechanism. The influence
asmuchas8%fromunity,demonstratingthenonlinearityof of the concentration of sulfate-ammonium-nitrate aerosols
suchchemicalsystems. on such rates is not currently accounted for in the adjoint
The above test was reassuring, yet limited in scope for model. We assume such an effect is less than 5% (Liao
a global CTM. To test our adjoint model over a wide vari- etal.,1999;Martinetal.,2003),especiallyastheabsorbing
etyofchemicalconditions,wealsocomparetheaccuracyof aerosols(blackcarbon,mineraldust)arenotactivevariables
theadjointderivativesofthechemicalmechanisminglobal duringthesetests. Thegeneralagreementbetweenλand3,
simulations over much longer time scales. We turn off all onlythelatterofwhichaccountsforthiseffect,indicatesthis
transportrelatedprocessesinthemodelandcalculatethead- assumptionisadequate,atleastforsimulationsofthislength.
jointandfinitedifferencesensitivitiesofsurfaceleveltracer Further tests indicate that this assumption is valid for most,
masseswithrespecttoNOx emissionsineachlocationafter thoughnotall,cases,seeSect.3.5.
a week-long simulation. As lack of transport leads to unre-
alisticallyextremeconcentrations,emissionsarereducedby 3.3 Convection,turbulentmixing,andwetremoval
anorderofmagnitudetopreventthechemicalsystemsfrom
Wet removal of tracers in GEOS-Chem is generally treated
becomingtoostiff. Manychemicalchangesassociatedwith
as a first-order process, leading to discrete forward model
aerosols are treated separately from the main tropospheric
equationsoftheform,
chemistry mechanism in GEOS-Chem, such as aqueous re-
actions,drydeposition,chemicalaging,andemissionofSOx cn+1 =cne−rw,k4t (19)
k k
andNH (Parketal.,2004). Theadjointsoftheseprocesses
3
are constructed separately (manually and with TAMC) and Sincethelossraterw,k formostspeciesdoesnotdependon
includedinthefollowingtests. any active variables (Jacob et al., 2000), the corresponding
Figure 3 shows the adjoint and finite difference sensitivi- adjointissimply
ties of several species with respect to surface level, anthro- λn =λn+1e−rw,k4t (20)
pogenicNO emissionsscalingfactors. Wechoosetoshow k k
x
sensitivities of species such as acetone and methacrolein to The adjoints of these routines are generated using hand-
NO emissionstoalsohighlightthepotentialvalueofthead- created code, retaining efficiency and legibility. However,
x
jointmodelforanalysisofnon-aerosolspecies. Weseefrom the in-cloud formation and cycling of sulfate aerosol from
these, andsimilartestsforotheractivespecies(notshown), SO is decidedly nonlinear, as the soluble fraction of SO
2 2
thatthesensitivitiescalculatedusingtheadjointmodelcon- islimitedbyavailabilityofH O ,andafractionoftheSO
2 2 2
sistentlyagreewiththoseusingthefinitedifferencemethod is reintroduced into the gas phase as sulfate when droplets
overawiderangeofconditions. evaporate (Park et al., 2004). Such nonlinearities that span
Atmos. Chem. Phys.,7,2413–2433,2007 www.atmos-chem-phys.net/7/2413/2007/
D.K.Henzeetal.: AdjointofGEOS-Chem 2419
s x 104
e
viti 8
ensiti 6 R2 =0.999
S
e m =0.995
nc 4
e
er
Diff 2
e
nit 0
Fi 0 2 4 6 8
-106 -105 -104 -103 -102 -101 101 102 103 104 105 106 x 104
Adjoint Sensitivities [kg / grid cell] Adjoint Sensitivities
[kg / grid cell]
s x 104
e
sitiviti 0 R2 =0.999
n
Se m =0.993
e -(cid:21)1
c
n
e
er -(cid:21)2
Diff
e
nit -(cid:21)3
Fi -3(cid:21) -2(cid:21) -(cid:21)1 0
-106 -105 -104 -103 -102 -101 101 102 103 104 105 106 x 104
Adjoint Sensitivities [kg / grid cell] Adjoint Sensitivities
[kg / grid cell]
s
e
viti 14000
siti R2 =0.999
n
e 10000
S(cid:21)
e m =0.992
c
n(cid:21)
e
er 5000
Diff
e
nit 0
Fi 0 5000 10000
-106 -105 -104 -103 -102 -101 101 102 103 104 105 106
Adjoint Sensitivities [kg / grid cell] Adjoint Sensitivities
[kg / grid cell]
Fig.3. Chemistryadjointvalidation. Intheleftcolumnaretheadjointsensitivitiesofsulfate(SO4),methacrolein(MACR),andacetone
(ACET)atthesurfacewithrespecttosurfacelevelanthropogenicNOxemissionsscalingfactors.Intherightcolumnaretheadjointgradients
comparedtofinitedifferencegradients. Thecostfunctionisevaluatedonceattheendofaweek-longsimulationwithonlychemistryand
emissions×0.1.
multiple program modules are treated both manually and matrixform,thisequationreads,
with the help of TAMC, requiring additional recalculation µ n+1  m1 ··· mL  µ n
andcheckpointingofintermediatevalues. k,1 mT mT k,1
Turbulent mixing in the boundary layer in the forward  ...  = ... ... ... · ...  (22)
model is calculated according to a mass-weighted mixing µ m1 ··· mL µ
k,L mT mT k,L
algorithm applied every dynamic time step (30min for our
Direct application of Eq. (11) yields the corresponding ad-
case),
jointequation,
µkn,+j1 = PLl=1mmTlµnk,l (21) λµ...k,1n =mm...T1 ·.·..· mm...T1 ·λµ...k,1n+1 (23)
λ mL ··· mL λ
µk,L mT mT µk,L
where µ is the mixing ratio (c/ρ, ρ is the density of air)
k,j whichcanbesimplywrittenas,
of tracer k in layer j, m is the air mass in a single layer l,
l
mT isthetotalairmassintheboundarylayercolumn,andL λn = mjPLl=1λnµ+k,1l (24)
is the number of layers in the boundary layer. Rewritten in µk,j m
T
www.atmos-chem-phys.net/7/2413/2007/ Atmos. Chem. Phys.,7,2413–2433,2007
2420 D.K.Henzeetal.: AdjointofGEOS-Chem
Deep convection is calculated in the forward model using where λ is the adjoint of the mixing ratio. Note that we
µ
cumulus cloud fluxes and an RAS type algorithm, see Ap- haveassumedthatthewinds(oranyothermetfields)arenot
pendixAofAllenetal.(1996). Wecalculatethediscretead- active variables; taking the adjoint with respect to the me-
jointofthisschemeusingTAMC,notingthatTAMCinitially teorologyisanothertaskinitself(see,forexample,Giering
generatescodethatisaccurate,yetseveralordersofmagni- etal.,2005). Applyingthesimpletransformλˆ =λ /ρ,and
µ µ
tude slower than necessary due to several superfluous loops substitutingthisintoEq.(28),wearriveatthefollowingad-
thathavetoberemovedmanually.Thenumericalschemefor jointequation,
theforwardcalculationiterativelysolvesasetofessentially ˆ ˆ
∂(ρλ ) ∂(ρλ u)
linear equations, with an internal time step of five minutes. − µ = µ (29)
∂t ∂x
Ifweneglectasingleconditionalstatementthatchecksonly
whichissimilarinformtoEq.(26). Ifweassumethatρ is
for rare floating point exceptions, then storage or recalcula-
relatively constant over a single dynamic time step and that
tionoftheintermediatevaluesisnotrequiredfortheadjoint
theadvectionislinear,thenwecansimplysolveEq.(29)us-
calculation.
ingthesamenumericalcodethatwasusedtosolveEq.(26)
Theadjointmodelperformanceforasimulationincluding
intheforwardmodel,scalingtheadjointby1/ρ beforeand
convection, turbulent mixing, and wet deposition is tested
re-scaling by ρ afterwards, which is equivalent to solving
by comparison of finite difference sensitivities to the ad-
Eq.(28).
joint sensitivities of concentrations of a soluble tracer with
While the continuous approach was in part adopted for
respect to its initial concentrations in a location exhibiting
reasons of practicality (the discrete advection algorithm in
strongconvection,deposition,andmixing. Horizontaltrans-
theforwardmodelnotbeingdirectlyamenableforusewith
port, chemistry, and aerosol thermodynamics are turned off
automaticdifferentiationtools),subsequentinvestigationin-
forthesetests. Weuseaperturbationofonepercentforthe
dicatesthatthecontinuousapproachissuitable,ifnotprefer-
finitedifferencecalculation. Theratioλ /3 forsimulations
c c
able. This is not surprising, as it is well documented that
that are 6h, 1d and 3d in length are 0.9998, 1.0002 and
discreteadjointsofsignpreservingandmonotonic(i.e.non-
1.0003,fromwhichweseeconsistentsatisfactoryagreement
linearanddiscontinuous)advectionschemesarenotwellbe-
between the two methods. Performance is similar in other
havedandcancontainundesirablenumericalartifacts,seefor
testedlocations.
exampleThuburnandHaine(2001),Vukic´evic´ etal.(2001),
3.4 Advection andLiuandSandu(2006)2.
To illustrate the benefits of the continuous adjoint ap-
Weimplementtheadjointofthecontinuousadvectionequa- proach for our system, the following numerical test is per-
tions. GEOS-Chem nominally employs a monotonic piece- formed. The sensitivity of aerosol concentrations with re-
wiseparabolic(PPM)advectionroutine(ColellaandWood- specttoconcentrationsinaneighboringcellsixhoursearlier
ward, 1984; Lin and Rood, 1996). Below we briefly show arecalculatedforameridionalcrosssectionofthenorthern
howthisschemecanbeusedtosolvethecontinuousadjoint hemisphere. Toaffordsimultaneouscalculationoffinitedif-
advection equations and afterwards address some of the is- ferenceandadjointsensitivitiesthroughoutthisdomain,only
suesweddedtothisapproach. Weconsiderthe1-Dexample horizontaladvectionintheE/Wdirectionisincludedinthese
oftheadvectionequationforatracerinmassconcentration tests. Figure4showsfinitedifferencesensitivitiescalculated
units, usingEq.(18)forseveralvaluesofδσ aswellastheadjoint
∂c ∂(uc) gradients.Theundesirablenatureofthefinitedifferencesen-
=− (25)
sitivities is indicated by negative sensitivities that have no
∂t ∂x
physicalmeaning. Thatnegativevaluesbecomemorepreva-
whereuisthewindvelocityinthex-direction. Theforward
lent as δσ→0 indicates such values are caused by disconti-
numericalmodelactuallysolvesthefluxformofEq.(25)in
nuitiesinthediscretealgorithm(ThuburnandHaine,2001).
termsofthemixingratio(LinandRood,1996),
Wecanexpectthatadjointsensitivitiesofthediscreteadvec-
∂(ρµ) ∂(ρµu)
=− (26) tionalgorithmwouldcontainsimilarfeatures,which,despite
∂t ∂x beingnumericallyprecisegradientsofthecostfunction,can
Assumingthatthecontinuityequationforρ issatisfied,this resultinconvergencetoundesirablelocalminimumsfordata
canberewrittenintheadvectionform, assimilation (Vukic´evic´ et al., 2001). Given the importance
∂µ ∂µ oftransportforanalysisofaerosols,useofthecontinuousap-
=−u (27)
∂t ∂x proachisdeemedpreferabletoimplementingalineartrans-
port scheme with well-behaved discrete adjoints at the cost
Applying the adjoint variable as a Lagrange multiplier and
offorwardmodelperformance.
integratingbyparts(see,forexample,AppendixAofSandu
etal.,2005a),thecontinuousadjointofEq.(27)is 2Liu,Z.andSandu,A.:AnalysisofDiscreteAdjointsofNumer-
∂λ ∂(λ u) icalMethodsfortheAdvectionEquation,Int.J.Numer.Meth.Fl.,
− µ = µ (28)
submitted,2006.
∂t ∂x
Atmos. Chem. Phys.,7,2413–2433,2007 www.atmos-chem-phys.net/7/2413/2007/
D.K.Henzeetal.: AdjointofGEOS-Chem 2421
(a) Continuous adjoint sensitivities (b) Finite difference sensitivities,
(c) Finite difference sensitivities, (d) Finite difference sensitivities,
Fig.4. Sensitivitiesofaerosolconcentrationswithrespecttoconcentrationsinadjacentcells6hearlierconsideringonlyE/Wadvection.
Sensitivitiesarecalculatedusing: (a)continuousadjointequationand(b)–(d)one-sidedfinitedifferencemethodwithperturbationsofδσ.
The finite difference sensitivities contain more extreme values, including physically meaningless negative sensitivities that become more
prevalentasδσ→0.
3.5 Combinedperformance
Againwecomparethegradientscalculatedusingtheadjoint
modeltothosecalculatedusingthefinitedifferencemethod,
thistimeincludingallmodelprocesses.Wecalculatethesen-
sitivityofglobalaerosoldistributionsofsulfate,ammonium,
andnitratetosurfaceemissionsofanthropogenicSO ,NO
x x
andNH inselectlocations. Asnotedpreviously,suchcom-
3
parisons are quite time consuming to perform on a global
scale owing to the expense of the finite difference calcula-
tions. Attemptingtocoverawiderangeofconditions,while
keeping the number of required calculations within reason, Fig.5. Selectpointsforaccuracytests. Blacklocationsusedfor
wechoosetoanalyzetenlocationsforeachsetofemissions anthropogenic emissions of SOx and NOx, grey points for NH3,
withoneoverlappingpairinEurope.
considered,seeFig.5.Thesimulationsareonedayinlength,
andthecostfunction(Eq.2)isevaluatedonlyonceattheend
oftheday. Weuseaperturbationofδσ=0.1andEq.(18)for
thefinite-differencecalculations. forthelargervalues. Asthegradientsinagivensetusually
Figure 6 shows the adjoint gradients compared to the fi- spanseveralordersofmagnitude,manyoftheslopesarebi-
nitedifferencegradientsforeachofninerelationships. From asedbyafewsuchlargervaluesandarenotrepresentativeof
visualinspectionofthescatterplots,itisclearthattheagree- theoverallfit.However,accountingforsuchheteroscedastic-
ment is generally within reason given the fact that using a itybyre-scalingthegradientsby1/porperformingweighted
continuous adjoint for advection is expected to cause some regressionsthatplacelessemphasisonthelargervaluesstill
amountofdiscrepancy. Regressionlines,slopes,andR2val- leadstothesamegeneralresults. Pickingtwiceasmanytest
uesaregivenforeachsetofcomparisons. Theabsolutedif- cells, differenttestcells, oradifferentvalueofδσ alsowas
ference between the two methods is often more substantial notfoundtosubstantiallyaltertheoverallcomparisons.
www.atmos-chem-phys.net/7/2413/2007/ Atmos. Chem. Phys.,7,2413–2433,2007
2422 D.K.Henzeetal.: AdjointofGEOS-Chem
Sulfate Nitrate Ammonium
x 104 x 104 x 104
0
12 2
R2 =0.962 R2 =0.997 R2 =0.964
8 -(cid:21)1
m =0.88 m =1.23 50 1 m =1.00 ESOx
s 4 -(cid:21)2 0
e
viti 4 8 12 -2 -1-5-050 00 50 00 1 2
nsiti x 105 x 104 x 105 x 104 x 104 x 104
Se 1 R2 =0.994 4 R2 =0.991 8 R2 =0.998 ENOx
ce 0 m =1.19 2 m =1.20 4 m =1.10 8000
n 4000
e (cid:21)-1 0
er 0 0 -20-0200000 4000 8000
ff (cid:21)-1 0 1 0 2 4 0 4 8
Di x 105 x 105 x 104
x 105 x 105
e
nit4000 R2 =0.963 15 R2 =0.927 4 R2 =0.974
Fi2000 m =1.30 -500 10 m =1.09 2 m =1.18 ENH3
-100 5
0 -15-0150-100-50 0 0 0
0 2000 4000 0 5 10 15 0 2 4
x 105 x 105
Adjoint Sensitivities [kg]
Fig.6. Fullmodelperformance. Comparisonofsensitivitiesofglobalaerosolburdens(kg)toanthropogenicprecursoremissionsscaling
factorscalculatedusingtheadjointmethodvs.thefinitedifferencemethod. Afewoftheplotscontaininsetswithmagnifiedviewsofa
clusterofpoints.
Initialcomparison(notshown)ofgradientsforfiveofthe numericaldiffusion,withspatialoscillationsofthesensitivi-
90 tests showed underestimation of adjoint sensitivities by tiesindicativeoferrorsduetotransport.
more than an order of magnitude. Four of these tests were Inourtests,transportdoesnotdrasticallydegradethecon-
for the sensitivity of sulfate with respect to NH3 emissions sistency of the correlation between the two approaches; all
while one was for the sensitivity of nitrate with respect to of the R2 are near unity. There is, however, some amount
SOx emissions. Using offline concentrations for calcula- of bias in the comparisons, as indicated by slopes ranging
tionofthecontributionofsulfateaerosoltophotolysisrates from0.8to1.3, andthisdoesappeartobearesultoftrans-
and heterogeneous reaction probabilities in the main tropo- port. Figure 7 contains scatter plots of the sensitivities of
sphericchemicalmechanismforthesetestsalleviatedthedis- sulfatewithrespecttoNO emissionsforseveraladditional
x
crepancy,demonstratingthatwhilethisfeedbackisgenerally tests. Panel (a) shows the results when advection is turned
negligible, it is occasionally quite strong. Future work will off. Thisleadstoimprovedagreement,m=1.03,comparedto
extendtheadjointmodeltoaccountforthisfeedback. thecenterleftpanelofFig.6; hence, thesourceofthisbias
Napelenoketal.(2006)performedacomplementaryanal- ispresumablyadvection.AsshowninFig.4,theadjointgra-
ysis on a regional scale, calculating the sensitivities of lo- dients are likely smoother and more physically meaningful
calaerosoldistributionswithrespecttodomain-wideprecur- thanthefinitedifferencesensitivities.
soremissionsovertheUnitedStateswithaforwardsensitiv- Toassesstheextenttowhichusingthecontinuousadjoint
itymethod(DDM-3D),usingfinite-differencecalculationsto of advection hinders this approach to validating the adjoint
check their results. While they found similarly good agree- modelasawhole,weperformadditionaltests,theresultsof
mentforthemoredirectrelationships(suchassensitivityof which are shown in Fig. 7. Including advection, but evalu-
sulfate with respect to SO emissions, or ammonium with ating the cost function only in a single location, rather than
2
respect to NH emissions), they had difficulty verifying the globally, leadstoaveryunsmoothadjointfieldandtriggers
3
variability in the sensitivities of some of the more indirect many nonlinear and discontinuous aspects of the numeri-
relationships(suchasthesensitivityofsulfatetoNH emis- cal scheme in a manner inconsistent with advection of the
3
sions or nitrate to SO emissions). Granted, they used the relatively smooth concentration field in the forward model;
2
more complex and rigorous thermodynamic model ISOR- hence,agreementbetweenadjointandfinitedifferencegradi-
ROPIA;they suggested thatsuchdiscrepancieswere dueto entsundertheseconditionsisworse,seepanel(b). Allofthe
Atmos. Chem. Phys.,7,2413–2433,2007 www.atmos-chem-phys.net/7/2413/2007/