Adjoint of the global Eulerian–Lagrangian coupled atmospheric transport model (A-GELCA v1.0): development and validation

We present the development of the Adjoint of the Global Eulerian–Lagrangian Coupled Atmospheric (AGELCA) model that consists of the National Institute for Environmental Studies (NIES) model as an Eulerian threedimensional transport model (TM), and FLEXPART (FLEXible PARTicle dispersion model) as the Lagrangian Particle Dispersion Model (LPDM). The forward tangent linear and adjoint components of the Eulerian model were constructed directly from the original NIES TM code using an automatic differentiation tool known as TAF (Transformation of Algorithms in Fortran; http://www.FastOpt.com), with additional manual preand post-processing aimed at improving transparency and clarity of the code and optimizing the performance of the computing, including MPI (Message Passing Interface). The Lagrangian component did not require any code modification, as LPDMs are self-adjoint and track a significant number of particles backward in time in order to calculate the sensitivity of the observations to the neighboring emission areas. The constructed Eulerian adjoint was coupled with the Lagrangian component at a time boundary in the global domain. The simulations presented in this work were performed using the A-GELCA model in forward and adjoint modes. The forward simulation shows that the coupled model improves reproduction of the seasonal cycle and short-term variability of CO2. Mean bias and standard deviation for five of the six Siberian sites considered decrease roughly by 1 ppm when using the coupled model. The adjoint of the Eulerian model was shown, through several numerical tests, to be very accurate (within machine epsilon with mismatch around to ±6 e) compared to direct forward sensitivity calculations. The developed adjoint of the coupled model combines the flux conservation and stability of an Eulerian discrete adjoint formulation with the flexibility, accuracy, and high resolution of a Lagrangian backward trajectory formulation. A-GELCA will be incorporated into a variational inversion system designed to optimize surface fluxes of greenhouse gases.

The tangent and adjoint components of the Eulerian model were constructed directly from the original NIES TM code using an automatic differentiation tool known as TAF (Transformation of Algorithms in Fortran; http://www.FastOpt.com), with additional manual pre-and post-processing aimed at improving the performance of the 10 computing, including MPI (Message Passing Interface). As results, the adjoint of Eulerian model is discrete. Construction of the adjoint of the Lagrangian component did not require any code modification, as LPDMs are able to track a significant number of particles back in time and thereby calculate the sensitivity of observations to the neighboring emissions areas. Eulerian and Lagrangian adjoint components were cou-

Introduction
Forecasts of CO 2 levels in the atmosphere and predictions of future climate depend on our scientific understanding of the natural carbon cycle (IPCC, 2007;Peters et al., 2007). To estimate the spatial and temporal distribution of carbon sources and sinks, inverse methods are used to infer carbon fluxes from geographically sparse observa- 5 tions of the atmospheric CO 2 mixing ratio (Tans et al., 1989). The first comprehensive efforts in atmospheric CO 2 inversions date back to the late 1980s and early 1990s (Enting and Mansbridge, 1989;Tans et al., 1989). With the increase in spatial coverage of CO 2 observations and the development of 3-D tracer transport models, a variety of numerical experiments and projects have been performed by members of the so-called "TransCom" community of inverse modelers (e.g., Law et al., 1996Law et al., , 2008Denning et al., 1999;Gurney et al., 2002Gurney et al., , 2004Baker et al., 2006;Patra et al., 2011). A number of studies have proposed improvements to the inverse methods of atmospheric transport (Kaminski et al., 1999b;Rödenbeck et al., 2003;Peters et al., 2005;Peylin et al., 2005;Chevallier et al., 2005;Meirink et al., 2008;Maki et al., 2010). Despite 15 progress in atmospheric CO 2 inversions, a recent intercomparison (Peylin et al., 2013) demonstrated the need for further refinement.
In recent decades, a density of observational network established to monitor greenhouse gases in the atmosphere has been increased, and measurements taken onboard ships and aircraft are becoming available (Bovensmann et al., 1999). However, 20 on a global scale CO 2 observation are not existing for many remote regions not covered by networks. This lack of data is one of the main limitations of atmospheric inversions, which can be filled by monitoring from space (Rayner and O'Brien, 2001). The satellite observation data from current (GOSAT, Kuze et al., 2009;Yokota et al., 2009;OCO-2, Crisp et al., 2004)  To link surface fluxes of CO 2 to observed atmospheric concentrations, an accurate model of atmospheric transport and an inverse modeling technique are needed. Generally, there are the Eulerian and the Lagrangian method of modelling the atmospheric constituents transport. The Eulerian method treats the atmospheric tracers as a continuum on a control volume basis, so it is more effective in reproducing of long-term 5 patterns, i.e. the seasonal cycle or interhemispheric gradient. The Lagrangian Particle Dispersion Models (LPDMs) consider atmospheric tracers as a discrete phase and tracks each individual particle, therefore LPDMs are better for resolving synoptic and hourly variations.
To relate fluxes and concentrations of a long-lived species like CO 2 a transport model 10 must cover a long simulation period (e.g., Bruhwiler et al., 2005). Therefore, computing time is a critical issue and minimization of the computational cost is essential. If tracer is a chemically inert, the transport can be represented by a model's Jacobian matrix, because the simulated concentration at observational sites is a linear function of the flux sets. To compute such a matrix a transport model is running multiple times with set 15 of prescribed surface fluxes. The adjoint of the transport model is an efficient way to evaluate derivatives of concentration of the simulated tracer at observational locations towards to the sources and sinks of tracer (Kaminski et al., 1999). Marchuk (1974) first applied the adjoint approach in atmospheric science. After that, this method became widely used in meteorology. In the 1990s the approach was ex-Introduction (CTMs) with high-resolution grids in inversions. It would be beneficial to increase the model resolution close to observation points, where small uncertainties in the transport can seriously improve optimization of the resulting emission fluxes. LPDM running in the backward mode can explicitly estimate a source-receptor sensitivity matrix by solving the adjoint equations of atmospheric transport (Stohl et al.,5 2009), which is mathematically a Jacobian expressing the sensitivity of concentration at observational locations. Marchuk (1995), and Hourdin and Talagrand (2006) discussed the equivalence of the adjoint of forward transport models to backward transport models.
To utilize of the strongest sides of both methods, Lagrangian and Eulerian chemical 10 transport models can be coupled to develop the adjoint, which is suitable for the simultaneous estimation of global and regional emissions. Coupling can be performed in several ways; e.g., a regional-scale LPDM can be coupled to a global Eulerian model at the domain boundary (Rödenbeck et al., 2009;Rigby et al., 2011) Thompson and Stohl, 2014). One goal of this study is to present the development and evaluation of an Adjoint of the Global Eulerian-Lagrangian Coupled Atmospheric model (A-GELCA), which consists of an Eulerian National Institute for Environmental Studies global Transport Model (NIES-TM; Maksyutov et al., 2008;Belikov et al., , 2013a) and a Lagrangian 20 particle dispersion model (FLEXPART; Stohl et al., 2005). This approach utilizes the accurate transport of the LPDM to calculate the signal near to the receptors, and rapid calculation of background responses using the adjoint of the Eulerian global transport model. In contrast to previous works (Rödenbeck et al., 2009;Rigby et al., 2011;Thompson and Stohl, 2014), in which the regional models were coupled at the spatial 25 boundary of the domain, we implemented a coupling at the time boundary in the global model domain.
The remainder of this paper is organized as follows. An overview of the coupled model is provided in Sect. 2, and in Sect. 3 we describe the variational inversion pro- cess. In Sect. 4 we address several problems regarding the coupled model that have not been covered previously (Ganshin et al., 2012). In Sect. 5 we describe the formulation and evaluation of the adjoint model. The computational efficiency of the adjoint model is analyzed in Sect. 6, and finally the conclusions are presented in Sect. 7.
2 Model and method 5

Global coupled Eulerian-Lagrangian model
In the paper we use a global Eulerian-Lagrangian coupled model, the principles of which are described by Ganshin et al. (2012). In this section we provide the formula in its discrete form, as implemented in the model for the case of surface fluxes: 10 where i , j , and k are the indices that characterize the position of the particle in the cell; l is the time index; p is the particle index; F l i j are the surface fluxes in kg m −2 s −1 ; C B i j k are the background concentrations in the Eulerian model; f n i j k equals unity if the particle is within cell i , j , k, otherwise it equals zero; T is the duration of the trajectory; L is the number of steps in time; N is the total number of particles; h is the height 15 up to which the effect of the surface fluxes is considered significant; ρ is the average air density below height h; and m air and m CO 2 are the molar masses of air and carbon dioxide, respectively. The FLEXPART model starts at the observation point and calculates seven days' worth of backward trajectories for 1000 air particles, which are dispersed under the influence of turbulent diffusion. The background grid values of the h (500 m). The value of the first term is proportional to the flux in each cell along the trajectory, and to the time during which the air particle is inside this cell (Ganshin et al., 2013). We implemented a coupling at the time boundary in the global domain. The coupled model consists of FLEXPART (version 8.0; run in backward mode) as the Lagrangian particle dispersion model, and NIES TM (version NIES-08.1i) as the 5 Eulerian off-line global transport model to calculate the background CO 2 values.

NIES transport model
Since the first publication of the GELCA model in 2012, the NIES transport model has undergone significant updates. We provide a brief outline of the major features of the current model. NIES TM is a global three-dimensional CTM that simulates the global sivity we follows the method proposed by Hack et al. (1993), with a separate evaluation of transport processes in free troposphere and the planetary boundary layer (PBL). The PBL heights are provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis. The modified Kuo-type parameterization scheme is used for cumulus convection (Belikov et al., 2013a). 20 Inverse modeling assumes that the model reasonably well reproduces the relationship between atmospheric mixing ration and surface fluxes, assuming that the biases between the simulated and observed concentrations are mostly due to the emission inventories errors. To ensure that this is the case, the NIES TM model has been evaluated extensively, and it consistently performs well in intercomparisons against SF 6 and Introduction

FLEXPART
FLEXPART similar to other LPDMs consider atmospheric tracers as a discrete phase and tracks pathways of each individual particle. The advantage of this approach is direct estimation the sensitivity of measurements to the neighboring sink and sources by running the particles back in time. Usually it is enough to simulate for a limited 5 number of days (2-10) to determine, where particles intercept the surface layer.

Meteorological data
To run both models we use reanalysis which combines the Japanese 25 yr Reanalysis (JRA-25) and the Japanese Meteorological Agency Climate Data Assimilation System (JCDAS) dataset (Onogi et al., 2007). The JRA-25/JCDAS dataset is distributed on a Gaussian grid T106 with horizontal resolution 1.25 • ×1.25 • , 40 sigma-pressure levels and in 6 h time steps. The use of JRA-25/JCDAS data for Eulerian and Lagrangian models provides a consistency in the calculated fields; however, some features of FLEX-PART and NIES TM require different methods for processing the meteorological data. 15 Isolation of the transport equations is an effective way to save a significant amount of CPU time during tracer transport simulation. At the preprocessing stage, the NIES TM core produced a static archive of advective, diffusive, and convective mass fluxes with time step similar to the one of original JRA-25/JCDAS data (6 h). After that the archive is used by an "offline" model specially designed only for passive transport of tracer. 20 Intermediate fluxes are derived by interpolation.

Meteorological data processing for NIES TM
Besides the mass fluxes, the static archives contain fields of temperature, pressure, humidity, vertical grid parameters (variation of the sigma-isentropic vertical coordinate over time), and others. The pre-calculated and stored data field can be used directly for any of the inert tracers. It is also possible to simulate chemically active tracers Introduction Approximately 20 3-dimensional and 1-dimensional arrays are written to a hard disk for every record. This comprises around 10 GB of data per month for the model's standard resolution of 2.5 • × 2.5 • . 5 Originally, FLEXPART was driving by ECMWF reanalysis dataset distributed on a grid with regular latitude-longitude horizontal structure and sigma-pressure vertical coordinate. Current version of the model was adapted to use JRA-25/JCDAS data, by horizontal bilinear interpolation of the required parameters from a Gaussian grid to a regular 1.25 × 1.25 grid. The vertical structure and temporal resolution of of JRA-25/JCDAS data were used without modification. Given large differences in structure, resolution and parameter estimation method used in different reanalysis dataset the use of the same meteorology for both Eulerian and Lagrangian models is a significant benefit.

Meteorological data processing for FLEXPART
3 Inverse modeling for the flux optimization problem 15 Although the variational inversion method theory for minimizing the discrepancy between modeled and observed mixing ratios has been well described and published (i.e. Chevallier et al., 2005), we summarize it here.
The aim of the inverse problem is to find the value of a state vector x with n elements that minimizes a cost function J(x) using a least-squares method: where y is a vector of observations with m elements, and the matrix H represents the forward model simulation mapping the state vector x to the observation space. Here, R is the covariance matrix (size m × m) for observational error, which includes instrument and representation errors. The matrix R also includes errors of the forward model H. B is the covariance matrix (size n × n) of error for prior information of the state vector x b . Use of the cost function in the form of Eq.
(2) assumes that all errors must have Gaussian statistics and be unbiased (Rodgers, 2000).

Assessment of the coupled model
The effect of different horizontal resolutions on Eulerian models is discussed in detail by Patra et al. (2008). In general, higher resolution helps to resolve a more detailed 15 distribution of the tracer. However, the use of a more detailed grid leads to additional computational cost, which is not always justified by the resulting model output. This is largely due to the limited availability of high-resolution meteorology and tracer emission datasets. Ganshin et al. (2012) in various test showed that the coupled model surpasses the 20 Eulerian model in 4-month simulations. The advantage of GELCA in reproducing the high-concentration spikes and short-term variations caused mainly by anthropogenic emissions is more vivid with use of high resolution (1 km × 1 km) surface fluxes compared to standard low-resolution (1 • × 1 • ) fluxes.
We (JR-STATION) as described in Table 2 (Sasakawa et al., 2010). Siberia is assumed to be a substantial source and sink of CO 2 emissions, with high uncertainties in the fluxes describing them (McGuire et al., 2009;Hayes et al., 2011;Saeki et al., 2013). As a result, CTMs tend to reproduce the interseasonal variability of CO 2 quite poorly.  putation costs between NIES TM for low-and high-resolution grids (i.e. a difference by a factor of ∼ 15 between grids with resolution 10.0 and 2.5 • ), the advantage of the GELCA model is clear. Performance is important, as the setup considered here is almost identical to the case used in the inverse modeling of CO 2 . This case shows that the coupled model is effective even in the case of multiple 5 limiting factors, such as high uncertainty of fluxes, a small number of observations, and the low resolution of the Eulerian model. We recognize that the use of the concentrations simulated from the highly uncertain surface fluxes to judge the quality of different model configurations is quite problematic. Nevertheless, we cannot improve our analysis, because we do not have concentration measurements for tracers with 10 more accurate fluxes, like SF 6 .

Construction
In this section, we present the development of the adjoint of the coupled model. Construction of the adjoint to the Lagrangian part does not require any modification to 15 the code, as LPDMs are able to track a significant number of particles backwards in time and thereby calculate the sensitivity of observations to the neighboring emissions areas.
The development of the adjoint to the Eulerian part is more complicated. We decided to develop a discrete adjoint of NIES TM in order to make it consistent with the forward 20 model. An alternative approach is a construction of continuous adjoint derived from the leading equations of the forward model (Giles and Pierce, 2000). The main advantage of the discrete adjoint model is that the resulting gradients of the numerical cost function are exact, even for nonlinear or iterative algorithms, making them easier to validate, as validation of the adjoint model is an essential and complicated task. Introduction The tangent linear and adjoint models for NIES TM were created using the Transformation of Algorithms in Fortran (TAF) software (http://www.FastOpt.com). Use of this tool required some manual treatment of the code. We often manually redesign and optimize the automatically generated adjoint code to optimize the efficiency and improve readability and clarity of the adjoint model. 5 The advantages of our coupled adjoint model are as follows.

Validation of the NIES TM adjoint
The discrete adjoint obtained through automatic differentiation can be easily validated by comparing the adjoint sensitivities with forward model gradients calculated using the finite difference approximation (Henze et al., 2007). Forward model sensitivity, λ, is calculated using the one-or two-sided finite difference where M denotes the tangent linear model. A range of ε = 0.1-0.01 proved in most cases to give an optimal balance between truncation and roundoff error (Henze at al., 10

2007).
In the first test, forward simulations were carried out with an initial CO 2 distribution and zero surface flux for 2 days using a horizontal grid with resolution 2.5 • ×2.5 • . Adjoint simulations were then performed with CO 2 distribution perturbed by 1 ppm per grid cell. The adjoint gradient was then compared with that from the finite difference calculated 15 using Eq. (3) in order to save CPU time by minimizing the number of forward model function calculation for the case ε = 0.01. To quantify the difference between the two calculations of sensitivity λ we define the local relative error

20
where the subscripts A and F refer to adjoint and finite difference respectively, lon, lat -longitude and latitude correspondently. Figure 3c shows E (lon, lat) with a logarithmic color scale. The sensitivities obtained by the adjoint have maximum relative error 5996 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | of order 10 −16 , indicating that transport in the NIES TM adjoint is correct over short timescales. The overall comparisons did not seriously changed if we select a different grid cells or use various values of ε.
The definition of the adjoint model M * requires that for an inner product , and two random vectors u and v , the following expression should be valid: For practical use the identity in Eq. (6) is reworded as follows (Wilson et al., 2014): We use Eq. (7) to test the adjoint model initialized using several different random setups. For all cases, Eq. (7) compares well with machine epsilon.

Real case simulation
The next series of calculations was made for real measurements. As in the first part of the article, we used data from the Siberian observation network (Table 2) for the period 01-04 January 2010. The NIES adjoint was simulated with a horizontal resolution of 2.5 • × 2.5 • , and the Lagrangian response was simulated with a horizontal resolution of The FLEXPART model sensitivity shows more irregular distributions, and higher values closer to the observational sites, thereby reflecting the model's ability to monitor small-scale changes.
During coupling, the sensitivity is aligned due to the crosslinking of components (Fig. 6). Thus, intensity has a maximum near the stations and smoothly decreases 5 with increasing distance. The Eulerian and Lagrangian models employ different approaches and grid resolutions for the modeling of atmospheric tracers, and can thus resolve processes with different time and spatial scales, and underlying physics. Figure 7 shows the sensitivity calculated for the same setup as for Fig. 6, but using NIES TM with a 10.0 • resolution. By changing the Eulerian model resolution, it is pos-10 sible to change size of the footprint. This system can utilize responses calculated at higher resolutions, such as 0.5 or 0.1 • , but these setups require more accurate driving data and regular observations available for smaller time steps.

Computational efficiency
We tested several different methods to reduce the computational burden of the adjoint 15 model. First, the Eulerian part of the adjoint model was driven by static archives of meteorological parameters, as described in Sect. 2.4.1. Second, the forward NIES model was altered so that at each model timestep it saved any variables that would also be needed by the adjoint model. These variables therefore did not have to be recalculated for use in the adjoint model. (This was possible because we used a discrete version of 20 the adjoint, which was fully compatible with the forward model.) Third, the Lagrangian part of the adjoint model made use of pre-calculated response functions, as described in Sect. 2.4.2.
To run the adjoint model we used a Linux workstation with 8 Intel ( is about 1 GB. Henze et al. (2007) reports that the ratio between simulation time in backward and forward modes for adjoint models derived for other CTMs, as follows: GEOS-Chem: 1. 5,STEM: 1.5,and CIT: 11.75. Thus, the adjoint of the developed coupled model is quite efficient. To achieve this level of efficiency, a substantial amount of manual programming effort is required, 5 despite the automatic code generated by TAF. The main disadvantage of TAF is that many redundant recomputations are often generated by the compiler. A crucial optimization of the adjoint code is required to eliminate these extra recomputations.

Summary
In this papers we have presented the construction and evaluation of an adjoint of the 10 global Eulerian-Lagrangian coupled model that will be integrated into a variational inverse system designed to optimizing surface fluxes. The coupled model combines the NIES three-dimensional transport model as its Eulerian part and the FLEXPART plume diffusion model as its Lagrangian component. The model was originally developed to study the carbon dioxide and methane atmospheric distribution.

15
FLEXPART is tracking a significant number of particles back in time, and thereby calculates the sensitivity of observations to the neighboring emissions areas. Thus, construction of the adjoint to the Lagrangian part does not require any modification to the code.
For Eulerian part the discrete adjoint was constructed directly from the original NIES 20 TM code, in contrast to a construction of continuous adjoint derived from the forward model basic equations. The tangent and adjoint models of the Eulerian model were derived using the automatic differentiation software TAF (http://www.FastOpt.com), which significantly accelerated the development. However, considerable manual processing of forward and adjoint model codes was necessary to improve transparency of the 25 model and to optimize the performance of computing, including MPI. The TAF code used here (version 1.5) did not fully support MPI routines. The Eulerian and Lagrangian adjoints were coupled at the time boundary in the global domain. The main benefit of the developed discrete adjoint is accurate calculation of the numerical cost function gradients, even if the algorithms are nonlinear. The overall advantages of the developed model also include relatively simple construction of the adjoint to the Lagrangian part and computation using the Lagrangian component, 5 scalability of sensitivity calculation depending on distance to monitoring sites, thereby reducing aggregation errors, and computational efficiency even for high-resolution simulations.

GMDD
The accuracy of the transport scheme of the forward coupled model was investigated using simulation of distribution of the atmospheric CO 2 . The GELCA components and the model itself had previously been validated in various tests and by comparison with both measurements and other transport models for CO 2 and other tracers. Performed in the paper analyses showed, that GELCA is effective in capturing the seasonal variability of atmospheric tracer at observation sites. Decreasing of the Eulerian model resolution are not able to significantly distort the transport model performance, how-15 ever, running the coupled model using NIES TM with low resolution grid can maximize simulation speed and use of data storage.
The Eulerian and Lagrangian components of the adjoint model were validated using various tests in which the adjoint gradients were compared to gradients calculated with numerical finite difference. We evaluated each individual routine of discrete adjoint of 20 Eulerian model and the adjoint gradients of the cost function. The obtained precision of the results in considered numerical experiments demonstrates proper construction of the adjoint.
The CPU time of the adjoint model is comparable with that of other models, as we used several methods to reduce the computational load. The forward NIES model was 25 altered so that at each model timestep it saved any variables that would also be needed by the adjoint model. These variables therefore did not have to be recalculated for use in the adjoint model. In addition, the adjoint simulation was isolated from the recalculation of NIES TM meteorological parameters and Lagrangian response functions.