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/