Articles | Volume 17, issue 4
Development and technical paper
22 Feb 2024
Development and technical paper |  | 22 Feb 2024

Development of the tangent linear and adjoint models of the global online chemical transport model MPAS-CO2 v7.3

Tao Zheng, Sha Feng, Jeffrey Steward, Xiaoxu Tian, David Baker, and Martin Baxter

We describe the development of the tangent linear (TL) and adjoint models of the Model for Prediction Across Scales (MPAS)-CO2 transport model, which is a global online chemical transport model developed upon the non-hydrostatic Model for Prediction Across Scales – Atmosphere (MPAS-A). The primary goal is to make the model system a valuable research tool for investigating atmospheric carbon transport and inverse modeling. First, we develop the TL code, encompassing all CO2 transport processes within the MPAS-CO2 forward model. Then, we construct the adjoint model using a combined strategy involving re-calculation and storage of the essential meteorological variables needed for CO2 transport. This strategy allows the adjoint model to undertake a long-period integration with moderate memory demands. To ensure accuracy, the TL and adjoint models undergo vigorous verifications through a series of standard tests. The adjoint model, through backward-in-time integration, calculates the sensitivity of atmospheric CO2 observations to surface CO2 fluxes and the initial atmospheric CO2 mixing ratio. To demonstrate the utility of the newly developed adjoint model, we conduct simulations for two types of atmospheric CO2 observations, namely the tower-based in situ CO2 mixing ratio and satellite-derived column-averaged CO2 mixing ratio (XCO2). A comparison between the sensitivity to surface flux calculated by the MPAS-CO2 adjoint model with its counterpart from CarbonTracker–Lagrange (CT-L) reveals a spatial agreement but notable magnitude differences. These differences, particularly evident for XCO2, might be attributed to the two model systems' differences in the simulation configuration, spatial resolution, and treatment of vertical mixing processes. Moreover, this comparison highlights the substantial loss of information in the atmospheric CO2 observations due to CT-L's spatial domain limitation. Furthermore, the adjoint sensitivity analysis demonstrates that the sensitivities to both surface flux and initial CO2 conditions spread out throughout the entire Northern Hemisphere within a month. MPAS-CO2 forward, TL, and adjoint models stand out for their calculation efficiency and variable-resolution capability, making them competitive in computational cost. In conclusion, the successful development of the MPAS-CO2 TL and adjoint models, and their integration into the MPAS-CO2 system, establish the possibility of using MPAS's unique features in atmospheric CO2 transport sensitivity studies and in inverse modeling with advanced methods such as variational data assimilation.

1 Introduction

Estimating CO2 fluxes through inverse modeling, using atmospheric chemical transport models and atmospheric CO2 measurements, is an important approach for understanding the global carbon budget. Beyond providing seasonal flux estimates that are useful for understanding the magnitude and phase of photosynthesis and respiration, it provides annual mean flux estimates that shed light on the key processes driving the response to climate change. When these annual mean CO2 estimates are adjusted to account for lateral fluxes (e.g., due to rivers or the transport of crops and wood products), it gives an independent means of validating carbon stock change estimates from the terrestrial biogeochemical models and inventories (Byrne et al.2023). However, atmospheric transport models, which play a key role in inverse modeling, remain a significant source of uncertainty at both regional and global scales (Schuh et al.2019, 2022; Hurtt et al.2022).

Two classes of chemical transport models – online and offline – are commonly used for simulating atmospheric CO2 transport. Offline models, such as TM5 (Krol et al.2005; Meirink et al.2006), PCTM (Kawa et al.2004; Baker et al.2006), and GEOS-Chem (Kopacz et al.2009), solve the tracer continuity equation using winds and vertical mixing fields computed from an independent run of a meteorological model or from a meteorological analysis. Online models, such as WRF-Chem (Grell et al.2011), OLAM (Walko and Avissar2008; Schuh et al.2021), and MPAS-CO2 (Skamarock et al.2012; Zheng et al.2021), integrate chemistry, transport, and meteorology simultaneously. Although offline models typically have lower computational costs, the separation of chemistry and transport from meteorology leads to a loss of information regarding atmospheric processes occurring at timescales shorter than the meteorological model output frequency (Grell et al.2005). In comparison, online models, owing to their simultaneous integration of meteorology and chemistry, have the potential to improve transport accuracy, particularly for vertical transport of chemistry. Recent advances in computer power and parallelization have greatly reduced the computational cost of online transport models, making them increasingly more accessible and practical for atmospheric CO2 research.

A number of studies have demonstrated that transport model accuracy can be improved by increasing the model's horizontal resolution (Feng et al.2016; Agusti-Panareda et al.2019). Because global high-resolution CO2 transport simulations are computationally demanding, limited-area models (regional models) are often used instead (Pillai et al.2012; Lauvaux et al.2012; Zheng et al.2018). However, regional models introduce the lateral boundary condition, posing challenges for CO2 inverse modeling (Zheng et al.2019; Rayner et al.2019). As a global online chemical transport model, the Model for Prediction Across Scales (MPAS)-CO2 (Zheng et al.2021) avoids the lateral boundary condition problem. Like OLAM (Schuh et al.2021), MPAS-CO2 uses a global variable-resolution mesh to facilitate local grid refinement for high-resolution simulations in specific regions without incurring prohibitively high computational costs and avoiding the disadvantages of lateral boundary conditions.

The primary objective of this study is to develop the tangent linear (TL) and adjoint (AD) models associated with the global online transport model MPAS-CO2 (Zheng et al.2021). Adjoint model techniques have been widely used in both meteorological and atmospheric greenhouse gas research (Errico1997; Courtier et al.1994; Giering et al.2006; Meirink et al.2008; Henze et al.2007; Tian and Zou2021) and play critical roles in variational data assimilation and sensitivity analyses (Baker et al.2006; Zheng et al.2018; Tian and Zou2020).

The subsequent sections of this paper provide an overview of the MPAS-CO2 forward model developed in Zheng et al. (2021) (Sect. 2) and the development and verification of the TL and AD models based on the forward model (Sects. 3 and 4). The utility of the newly developed AD model is demonstrated with adjoint sensitivity analyses in Sect. 5. Finally, a summary and conclusions are given in Sect. 6.

2 MPAS-CO2 forward model

Zheng et al. (2021) documented the development of MPAS-CO2, verifying its mass conservation and assessing its accuracy. Hereafter, we refer to MPAS-CO2 as the forward model, whose TL and AD model counterparts we develop in the present paper. A brief description of the forward model is provided here; see Zheng et al. (2021) for comprehensive details. The forward model characterizes CO2 transport through the continuity equation, as follows:

(1) ρ ̃ q CO 2 t = - ρ ̃ q CO 2 V ζ + F bl + F cu ,

where qCO2 is CO2 dry-air mixing ratio, ρ̃=ρd/(ζ/z), ρd is dry-air density, ζ is the vertical coordinate, z is geometric height, t is time, and V=(u,v,w) is the velocity vector (u, v, and w are the zonal, meridional, and vertical wind components, respectively). The meteorological variables, such as wind velocity and dry-air density, are updated simultaneously with CO2 by the model's dynamical core and physics parameterizations. The left-hand side (LHS) of Eq. (1) is the total CO2 time tendency ((ρ̃qCO2)/t). The first term on the right-hand side (RHS) represents the contributions to CO2 time tendency from the advection. The second (Fbl) and third (Fcu) terms of RHS represent the contribution from the vertical mixing by the planetary boundary layer (PBL) and cumulus convective transport parameterizations, respectively. Advection of CO2 in MPAS-CO2 is handled in the model's dynamical core and can be expressed as Eq. (2), where the first two terms on the RHS represent the horizontal advection, and the third term represents the vertical advection:

(2) ρ ̃ q CO 2 V ζ = ( ρ ̃ u q CO 2 ) x + ( ρ ̃ v q CO 2 ) y ζ + ρ ̃ w q CO 2 ζ .

CO2 vertical mixing by the PBL parameterization is implemented based on the YSU scheme (Hong et al.2006) and can be expressed as

(3) q CO 2 t bl = z K h q CO 2 z - w q CO 2 h z h 3 ,

where z is the vertical distance to the surface, h is the boundary layer top height, and Kh is the vertical eddy diffusivity. The second term on the RHS in the square bracket of Eq. (3) represents the contribution from CO2 entrainment flux at the inversion layer. The term [qCO2/t]bl from Eq. (3) is coupled with dry-air density ρ̃ to form the term Fbl of Eq. (1).

Convective transport of CO2 is implemented based on the Kain–Fritsch convection scheme (Kain2004), and it can be expressed as Eq. (4):

(4) q CO 2 t cu = ( M u + M d ) ρ A q CO 2 z + M ud M q CO 2 u - q CO 2 + M dd M q CO 2 d - q CO 2 ,

where qCO2, qCO2u, and qCO2d are the CO2 mixing ratio in the environment, updraft, and downdraft, respectively; Mu and Md are the updraft and downdraft mass, respectively; ρ is the environment air density; A is the horizontal area of a cell; M=ρAδz is the mass of environmental air in a grid box; and Mud and Mdd are the detrainment from the updraft and downdraft, respectively. The term [qCO2/t]cu of Eq. (4) is coupled with dry-air density ρ̃ to form the term Fcu of Eq. (1).

3 Development of the MPAS-CO2 TL model

3.1 TL model development

The CO2 advective transport process described in Eq. (2) is implemented by two different numerical schemes in the forward model: (1) a monotonic scheme with hyperviscosity (β) set to 0.25; and (2) a non-monotonic scheme with β=1.0 (Skamarock et al.2012). The monotonicity in the first scheme is achieved by applying a flux limiter in the last step of the third-order Runge–Kutta solver (Wang et al.2009; Skamarock and Gassmann2011). While the second scheme is linear in CO2, the first scheme is nonlinear due to the application of the flux limiter. Because both the YSU PBL and Kain–Fritsch convection schemes are linear in CO2, using the linear advective scheme makes the forward model a linear model in CO2. In this paper, we develop the TL and adjoint models based on the linear version of the MPAS-CO2 forward model, which can be symbolically expressed as

(5) x t = M ( x 0 , e ) ,

where x0 and xt are the CO2 dry-air mixing ratio at the initial and forecast time (t), respectively. ℳ( ) represents the MPAS-CO2 forward model, and e represents a time series of CO2 fluxes between times 0 and t. While both x0 and xt are 3-dimensional in space, e is 2-dimensional in space, indicating that CO2 flux is applied only to the model's surface cells. Equation (5) indicates that CO2 mixing ratio at a forecast time (xt) is determined by the CO2 mixing ratio at an initial time (x0) and the CO2 flux (e) through the forward model.

The TL and adjoint models are designed to calculate the sensitivity of xt with respect to x0 and e. This is achieved by introducing the TL and adjoint variables of their counterparts in the forward model (Giles and Pierce2000). While the introduction of the TL and adjoint variables for the initial CO2 mixing ratio (x0) is straightforward, it is a bit more complex for the CO2 fluxes (e). This complexity arises from the fact that CO2 flux, at each surface cell of the model, varies with time throughout the model's entire simulation period. Depending on the underlying biosphere model and emission inventory used, CO2 flux varies at a certain temporal frequency, ranging from hourly to monthly. Although it is possible to introduce TL and adjoint variables for CO2 flux at the flux's temporal frequency, it is neither practical nor necessary to do so. Instead, a common approach is to introduce flux scaling factors (Henze et al.2007; Zheng et al.2018) as follows:

(6) e = S ( k ) e ̃ ,

where ẽ are time-variant CO2 fluxes, typically from a process model or inventory, and S(k) is a generic scaling function. Equation (6) means that at each surface cell, the magnitude of the CO2 flux (ẽ) is adjusted using a flux scaling factor before it is used to modify the cell's CO2 mixing ratio. We implemented Eq. (6) in the forward model in a way that allows the flexibility of choosing the temporal frequency of the flux scaling factor. For instance, for a 24 h forward model simulation forced by 3 h CO2 flux, one can choose to have eight scaling factors at each surface cell (one for each of the eight 3 h segments) or just one scaling factor for the entire time period. All the MPAS-CO2 model runs used in the remainder of this paper are conducted using a single scaling factor for each surface cell that is repeated for each flux time step in the entire simulation period. In this case, the scaling function S(k) in Eq. (6) is a function of a scaling vector k that has the same dimension as the model's surface mesh. The introduction of the flux scaling factors turns CO2 flux from active variables to parameters, and the impacts of their variation on CO2 mixing ratio are calculated through their corresponding scaling factors k. Accordingly, the MPAS-CO2 forward model can be symbolically expressed as

(7) x t = M ( x 0 , k ) .

Equation (7) shows that for a given set of CO2 flux (ẽ), the forecast time CO2 mixing ratio (xt) is a function of the initial time CO2 mixing ratio (x0) and the flux scaling factor (k).

Table 1Results of the correctness check for the newly developed MPAS-CO2 tangent linear model. The results are from 1 month integration (from 1 October 2018 at 00:00 UTC to 1 November 2018 at 00:00 UTC) of the forward and tangent linear models, using the 120–480 km global variable-resolution mesh (Fig. 1). The terms in the table refer to Eq. (9).

Download Print Version | Download XLSX

The TL counterpart of the MPAS-CO2 forward model represented by Eq. (7) can be symbolically expressed as the first derivative of the forward model:

(8) Δ x t = M ( Δ x 0 , Δ k ) ,

where M( ) represents the MPAS-CO2 TL model, Δx0 and Δxt are the TL variable of CO2 mixing ratio at the initial and forecast time, respectively, and Δk is the TL variable of the flux scaling factor k. In essence, Eq. (8) shows that the TL model computes the perturbation in the forecast time CO2 mixing ratio (Δxt), given the perturbation in the flux scaling factor (Δk) and/or perturbation in the initial time CO2 mixing ratio (Δx0).

Based on the source code of the forward model, we developed the TL code by differentiating each process relevant to CO2 flux and transport, including advection, vertical mixing by the YSU PBL scheme, convective transport by the Kain–Fritsch scheme, and the CO2 emission driver that implements Eq. (6). Automatic differentiation tools, such as Tapenade (Hascoet and Pascual2013) and Tangent and Adjoint Model Compiler (Giering and Kaminski1998), can be used to assist TL and adjoint code generation. However, the code these tools generate typically contains redundancies and is difficult to read, particularly for the adjoint code. To optimize the computation efficiency and facilitate future code upgrading, we manually developed the TL and adjoint code for MPAS-CO2 with some minor assistance from Tapenade.

3.2 TL model validation

After the TL model is completed, a thorough examination of its correctness was undertaken. As indicated in Eq. (8), the TL model can calculate the sensitivity of xt with respect to both x0 and k. The calculation of the sensitivity of xt with respect to x0 involves the TL code of all the CO2 transport processes, including advection, PBL, and convective transport. In comparison, the calculation of the sensitivity of xt with respect to the flux scaling factor k involves the TL code of the CO2 emission driver in addition to the TL code of all the CO2 transport processes. Because the calculation of sensitivity to k includes the TL code of all the processes in the TL model and because both the transport processes and emission driver are linear, the correctness of the entire MPAS-CO2 TL model can be verified by checking whether the following equation is satisfied (Errico1997; Tian and Zou2020):

(9) Φ ( α ) = M ( x 0 , ( 1 + α ) k ) - M ( x 0 , k ) M ( 0 , α k ) = 1 ,

where M() is the TL model, ℳ() is the forward model, and α is a scalar. The second item in the numerator of Eq. (9), ℳ(x0,k), is a forward model run. The first item in the numerator, M(x0,(1+α)k), is an identical forward model run, except that its flux scaling factor at each surface cell is adjusted by multiplying 1+α. In the denominator, M(0,αk) is a TL model run with its perturbation in the initial time CO2 mixing ratio set to zero (Δx0=0) and perturbation in flux scaling factor Δk=αk, which is the difference in the flux scaling factors between the two forward model runs.

Figure 1A MPAS-CO2 global variable resolution mesh ranging from  120 km over most of land regions to  480 km over oceans. Also shown in the figures are the ground tracks of two Orbiting Carbon Observatory-2 orbits, which are used for the adjoint sensitivity studies described in Sect. 5. The blue-colored ground track crosses North America from the Caribbean Sea northward between 18:31 UTC and 18:48 UTC on 30 June 2016. The red-colored ground track crosses from South America to North America between 17:36 UTC and 18:13 UTC on 23 August 2016.

If the TL model is correctly coded with regard to the forward model, Eq. (9) should be satisfied to the extent of machine accuracy until α is too small, so that the result is affected by rounding-off errors and drifts away from unity. To verify using Eq. (9), we ran a series of simulations using the forward model and newly developed tangent linear model with the scale factor α ranging from 1.0 × 103 to 1.0 × 10−4 (Table 1). All of the simulations start on 1 October 2018 at 00:00 UTC, run for 1 month, and end on 1 November 2018 at 00:00 UTC. The meteorological initial condition is from the ERA5 reanalysis (Hoffmann et al.2019), and the CO2 initial condition (x0) is from CarbonTracker (Jacobson et al.2020) v2022 (CT2022) posterior CO2 mole fraction at this time. The 3 h CO2 fluxes for the biogenic, fire, fossil fuel, and oceanic components from the CT2022 posterior are applied throughout the 1-month simulation period for each model run. Flux scaling factors of k=1 were used in all our simulations here, with 1 being a vector the same length as k, with ones in every element. The model simulations are conducted using the global variable-resolution (VR) mesh shown in Fig. 1. This VR mesh has a total of 15 898 cells, which range from 120 km over most of the land regions to 480 km over oceans. Table 1 shows that the magnitudes of both the numerator and the denominator in Eq. (9) decrease as α decreases. Moreover, the table also shows that the ratio remains close to unity until α decreases to 1.0 × 10−1, beyond which rounding-off errors lead to a deviation from unity. These results confirm that the MPAS-CO2 TL model has been correctly developed with regard to the forward model. In the next section, we proceed to develop the MPAS-CO2 adjoint model.

4 Development of the MPAS-CO2 adjoint model

4.1 Adjoint model development

An adjoint model is an essential component of a variational data assimilation system and is very useful for adjoint sensitivity analysis (Tian and Zou2021; Zheng et al.2018; Bosman and Krol2023). Symbolically, the MPAS-CO2 adjoint model can be expressed as

(10) Δ x ^ 0 , Δ k ^ = M T Δ x ^ t ,

where MT() is the MPAS-CO2 adjoint model, Δk^ is the adjoint variable of the flux scaling factor, and Δx^0 and Δx^t are the adjoint variables of CO2 mixing ratio at the initial and forecast time, respectively. Equation (10) shows that starting with Δx^t at the forecast time, the MPAS-CO2 adjoint model runs backward in time to the initial time, while ingesting CO2 observations along the way, resulting in the adjoint variable of CO2 mixing ratio at the initial time (Δx^0) and the adjoint variable of the flux scaling factor (Δk^).

Similar to its TL model counterpart, the development of the MPAS-CO2 adjoint model was carried out through manual implementation to avoid redundancy and optimize computational efficiency. The calculation of the CO2 transport needs access to the meteorological fields at each time step. Since the forward and TL model both run forward in time, this access is straightforward. However, because the adjoint model runs backward in time, accessing the meteorological fields is more challenging. One approach to this problem is saving meteorological fields in memory during the adjoint model's forward sweep, enabling accessing during the subsequent backward sweep (Guerrette and Henze2015; Zheng et al.2018). However, since the MPAS-CO2 adjoint model is intended for long simulations, this approach becomes impractical due to the excessive memory it demands. As an alternative strategy, we adopt an approach that combines both recalculation and storage of the meteorological fields. This strategy effectively divides a long simulation into segments, and the forward and backward sweeps are carried out sequentially for each segment, requiring internal memory only large enough to accommodate one segment's worth of meteorological fields. This internal manipulation is handled seamlessly by the adjoint model, enabling it to run as long as needed without overburdening memory resources. Another strategy we adopted for developing the adjoint code is to have the forward sweep, save some immediate variables that are needed by the subsequent backward sweep so that they do not need to be recalculated. For instance, the values of some variables related to mass fluxes in the Kain–Fritsch convection scheme Kain (2004) are saved by the forward sweep in the memory to speed up the subsequent backward sweep execution. This strategy not only increases the adjoint model efficiency but also simplifies some of its code development.

Table 2Results of the correctness check for the newly developed adjoint model of MPAS-CO2. All simulations are of the 120–480 km variable-resolution mesh (Fig. 1). The LHS and RHS in the table refer to Eq. (11).

Download Print Version | Download XLSX

4.2 Adjoint model validation

The correctness of the newly developed MPAS-CO2 adjoint model can be verified using the following equation (Tian and Zou2020):

(11) Δ x , M ( 0 , Δ k ) = M T ( Δ x ) , Δ k ,

where 〈 〉 represents the inner product operator, Δx is a perturbation of CO2 mixing ratio, and Δk is a perturbation of CO2 flux scaling factor. If the adjoint model is correctly coded with respect to the TL model, then Eq. (11) should be satisfied for any choice of Δx and Δk. M(0,Δk) on the LHS of the equation is the perturbation in forecast CO2 mixing ratio resulting from a TL model run whose perturbation in initial CO2 mixing ratio is set to zero and perturbation to flux scaling factor is set to Δk. The first item of the RHS, MTx), represents the adjoint variable of flux scaling factor, which is an output from the adjoint model integration from the forecast time backward to the initial time. The TL and adjoint model runs on the two sides of Eq. (11) have the same simulation time period, but the latter runs backward in time.

We conducted two sets of experiments using the TL and adjoint models following Eq. (11) to verify the correctness of the newly developed adjoint model. In the first set of experiments, we set Δk=10-11, and Δx=M(0,Δk). The experiments were carried out in two steps. First, the TL model was integrated 7 d from the initial time (1 October 2018 at 00:00 UTC) to the end time (8 October 2018 at 00:00 UTC), with Δk=10-11, resulting in M(0,Δk) which is the perturbation in forecast time CO2 mixing ratio. Second, the adjoint model is initialized on 8 October 2018 at 00:00 UTC, with its adjoint variable for CO2 mixing ratio set to M(0,Δk). The adjoint model is then integrated backward in time for 7 d to 1 October 2018 at 00:00 UTC, resulting in MTx). The LHS and RHS of Eq. (11) are then calculated using the above results (Table 2). The table shows that the agreement between the LHS and RHS of Eq. (11) is about 5.16 × 10−15. We note that this value is not exactly zero due to the machine rounding errors. This experiment is repeated with the same configuration, but the simulation length is increased to 31 d, ending on 1 November 2018 at 00:00 UTC. As expected, the magnitude of both the LHS and RHS increased, and they agreed to about 2.55 × 10−16. In the second set of experiments, Δk=10-11 (same as the first set of experiments) but Δx=M14d(x0,k), which is the CO2 mixing ratio at the end of 14 d forward model run (1 October 2018 at 00:00 UTC to 15 October 2018 at 00:00 UTC). We note that this forward model run uses x0 from CT2022 posterior CO2 mole fraction, and k=1; however, Eq. (11) should satisfy for any configuration and simulation period of the forward model. The resulting LHS and RHS of Eq. (11) from the second set of experiments are about 2 orders of magnitude larger than their counterpart of the first experiments. This is caused by the much larger Δx of the second set of experiments. The LHS and RHS agree to about 3.42 × 10−15 for the 7 d simulation and about 2.66 × 10−15 for the 31 d simulation (Table 2).

The results shown in Table 2 obtained from the experiments based on Eq. (11) confirm that the MPAS-CO2 adjoint model has been correctly developed with regard to the TL model. As the TL model has already been confirmed correct with respect to the forward model, it follows that both TL and adjoint models are correct with respect to the forward model of MPAS-CO2. This validation ensures the reliability and integrity of the entire MPAS-CO2 model system, since the forward model has already been validated in Zheng et al. (2021). It allows MPAS-CO2 to be used as the basis of a variational assimilation system for carbon flux estimation and as a platform for conducting sensitivity analyses in atmospheric carbon research.

Table 3 presents the computational cost of model simulations using the MPAS-CO2 system. Using the global 120–480 km VR mesh (Fig. 1; 15 898 cells), the 1-month forward model simulation completes in 20 min when using 128 processors. Both the TL and adjoint model simulations using the same configuration take approximately 10 % longer than the forward model. This extra computation time for the TL/adjoint model is incurred by the execution of the TL/adjoint code of the CO2 transport processes. Furthermore, we conducted another set of 1-month simulations using the models on a global quasi-uniform resolution (UR) mesh of about 120 km, consisting of a total of 40 962 cells. Table 3 demonstrates that the simulations with the VR mesh reduce the computational cost by over 50 % for all three models, primarily due to its substantially smaller number of cells. This reduction in computation cost, while preserving the high resolution over areas of interest, should prove advantageous when the models are applied in variational assimilation problems, which typically require many iterations of forward and adjoint model runs.

Table 3The computational costs for a 30 d simulation of MPAS-CO2 at  120 km quasi-uniform resolution and at a variable resolution ranging from  120 to  480 km. Computational costs are shown for the forward, tangent linear, and adjoint models. All simulations are conducted using 128 AMD Epyc 7H12 2.595 GHz processors running in parallel.

Download Print Version | Download XLSX

5 Adjoint sensitivity analysis

5.1 Comparison with CT-L footprints

In addition to forming a key component of variational assimilation systems (Baker et al.2006; Zheng et al.2018; Tian and Zou2021), adjoint models are powerful tools for sensitivity analysis (Errico and Vukicevic1992; Errico1997; Zou et al.1997; Tian and Zou2020). Studies focused on carbon flux estimation are often interested in exploring the sensitivity of atmospheric CO2 measurements to surface CO2 fluxes. This sensitivity is commonly referred to as observation influence functions or footprints (Cui et al.2022). The computation of observation footprints using forward models requires a large number of model runs, making it impractical, except for optimizing state vector elements at coarse horizontal resolutions. In contrast, adjoint models can calculate observation footprints much more efficiently. For point measurements, such as those from tower data, Lagrangian dispersion models offer an efficient alternative for obtaining footprints (Lin et al.2003; Stohl et al.2005). For an example of this, see the publicly available CO2 observation footprints from CarbonTracker–Lagrange (CT-L) (Hu et al.2019), which are generated using the Lagrangian particle dispersion model STILT (Lin et al.2003), driven by meteorology generated by the Weather Research and Forecast (WRF) model (Skamarock et al.2008). This approach involves releasing a certain number of particles from the observation location/height and tracing their backward transport in time. Note that CT-L is a regional modeling system that only provides observation footprints within the latitude range 10–80 N and longitude 0–180 W for up to 10 d backward in time.

In this section, we perform sensitivity analyses using the MPAS-CO2 adjoint model, which employs backward-in-time integration to calculate two quantities: (1) the sensitivity of atmospheric CO2 to the model's initial CO2 mixing ratio and (2) the sensitivity to the surface flux scaling factor. When a spatially uniform time-invariant surface flux of a unity value is used, and S(k)=k in Eq. (6), the sensitivity to the surface flux scaling factor calculated by the MPAS-CO2 adjoint model is the observation footprint. To compare with CT-L footprints, all MPAS-CO2 adjoint model simulations conducted in this section use a time-invariant CO2 surface fluxes of 1.0 µmol (m2 s)−1 for all surface cells, including both land and ocean cells. Because S(k)=k is used, Eq. (6) takes the form of e=kẽ. The units of the CO2 flux (µmolm-2s-1) and the multiplicative nature of the flux scaling factor k determine the units of the adjoint variable Δk^ (which represents observation footprints) to be in ppm(µmolm-2s-1)-1. Meteorological initial conditions for MPAS-CO2 model simulations conducted in this section are from the ERA5 reanalysis (Hoffmann et al.2019). Footprints calculated by the MPAS-CO2 adjoint model are of the 120–480 km variable resolution grid, while the CT-L footprints are of 1× 1.

Figure 2Sensitivity of CO2 mixing ratio at the WKT tower at 00:00 UTC on 31 March 2018 to the initial CO2 mixing ratio (a–d; units in ppm ppm−1) and the surface flux scaling factor (e–h; units in ppm(µmolm-2s-1)-1. The four rows from top to bottom show the sensitivities at 5, 10, 20, and 30 d before the observation. The sensitivities to the initial CO2 mixing ratio (a–d) are plotted as the column average. The WKT tower (31.3149 N, 97.3269 W) measurements used here are taken at 457 m above the ground level and labeled by the red color cross in the figures of the left column (a–d).

Figure 3The variation over time in the standard deviation (σ) of two quantities: the sensitivity to the initial CO2 mixing ratio and the sensitivity to the flux scaling factor (footprint). The standard deviations were calculated from MPAS-CO2 adjoint model simulations starting on 31 March 2018 at 00:00 UTC, running 30 d backward in time, and ending on 1 March 2018 at 00:00 UTC. Panel (a) is for the WKT tower (457 ma.g.l.) while panel (b) is for the WGC tower (483 ma.g.l.). In each figure, the triangles represent the standard deviation of sensitivity to the CO2 mixing ratio field (units in ppm ppm−1), and the circles represent the standard deviation of footprint (units in ppm(µmolm-2s-1)-1).


First, we conduct MPAS-CO2 adjoint model simulations for in situ CO2 observations at two towers in the United States, namely WKT, located at Moody, Texas (31.31 N, 97.33 W), and WGC, located at Walnut Grove, California (38.26 N, 121.49 W). For each tower, the adjoint model is initialized at 00:00 UTC on 31 March 2018. We add an adjoint forcing of 1 ppm CO2 at that time to the model grid cell closest to the tower location and the intake height (457 m at WKT and 483 m at WGC). The forcing is turned off for subsequent time steps, and the adjoint model is run backward in time for 30 d, ending at 00:00 UTC on 1 March 2018. The resulting sensitivity of CO2 at the WKT tower to the model's CO2 mixing ratio, which is 3-dimensional, is shown as a column average in the left column of Fig. 2. The right column of Fig. 2 shows the observation footprint (sensitivity to the surface flux scaling factor) at the corresponding times. The figure panels show that the sensitivity to the initial CO2 is highest and concentrated closest to the tower site at the time closest to the measurement, i.e., 5 d. With the increasing length of the backward-in-time integration, the sensitivity spreads over a larger area, and its magnitude decreases. After 30 d, the sensitivity to the initial condition has propagated across most of the Northern Hemisphere. Figure 2 also indicates that the spatial variation in the sensitivity magnitude decreases with time. To examine this, we calculated the standard deviation (σ) of sensitivity for each day of the 30 d (Fig. 3). The triangles in Fig. 3 show that the magnitude of the standard deviation of the sensitivity to the initial CO2 mixing ratio decreases rapidly with the increasing length of the adjoint model simulation for both towers. On the other hand, the sensitivity to the surface flux scaling factor (footprint) exhibits a different pattern from the sensitivity to the initial CO2 mixing ratio. As shown in Fig. 2, the footprint spreads spatially but the near-field to the tower maintains a much higher magnitude than the far-fields. By the end of 30 d, the footprint of the WKT tower covers almost the entire Northern Hemisphere, with the area north and northwest of the tower within the contiguous United States exhibiting a much higher magnitude than the more distant area. The circles in Fig. 3 indicate that the standard deviation of the footprint increases with time, but the rate of increase diminishes substantially after about 10 to 15 d. The finding suggests that extending the adjoint model integration further backward in time will still result in changes to the footprint but with a much-reduced change rate.

Figure 4The 10 d backward in time CO2 measurement footprint (units in ppm(µmolm-2s-1)-1) given by two tall towers, namely WKT and WGC. The figures in the top row are the footprint of the WKT tower calculated using the MPAS-CO2 adjoint model (a) and CT-L (c). The figures in the bottom row are the footprint of the WGC tower calculated by the MPAS-CO2 adjoint model (b) and CT-L (d). The location of the towers is marked by the black crosses in the figures in the right column (c, d).

For comparison, in Fig. 4, we plot the MPAS-CO2 adjoint model-calculated 10 d footprints in the CT-L geographic domain. The figure reveals that the MPAS-CO2 adjoint model-calculated WKT tower footprint spans most of the western and northwestern United States, with the highest sensitivity in Texas, Missouri, Iowa, Kansas, and Nebraska. Additionally, the footprint extends to a substantial area over the northeastern Pacific Ocean. The spatial pattern of the CT-L calculated footprint (Fig. 4c) is similar to that from the MPAS-CO2 adjoint model, but it is visibly less continuous. Figure 4 also shows that the MPAS-CO2 adjoint model-calculated footprint for the WGC tower covers northern California, Oregon, west Nevada, and a portion of the northeastern Pacific Ocean. The CT-L-calculated footprint exhibits a similar spatial pattern and magnitude. Overall, both the MPAS-CO2 adjoint model and CT-L provide valuable information on the sensitivity of atmospheric CO2 measurements to the surface flux; there are similar spatial patterns, although with some differences due to resolution and the Lagrangian/Eulerian framework difference.

In the second set of experiments, we compare CT-L and MPAS-CO2 adjoint model footprints for a swath of Orbiting Carbon Observatory-2 (OCO-2) XCO2 measurements. The ground track of the OCO-2 orbit used in the experiments is indicated by the blue line in Fig. 1. This orbit crosses North America from the Caribbean Sea to Canada's Northwest Territories in a northward direction between 18:31 UTC and 18:48 UTC on 30 June 2016. Since OCO-2 XCO2 represents the column average of atmospheric CO2, CT-L calculates XCO2 footprints at 14 discrete height levels, ranging from 50 to 14 000 m above the ground. For each height level, footprints are computed by placing a number of particles at that specific height. To ensure consistency with the CT-L approach, the MPAS-CO2 adjoint model is configured to apply the adjoint forcing at the corresponding vertical levels within the model. This is done by interpolating the CT-L's 14 height levels to MPAS-CO2 model's 55 vertical levels. This configuration allows for a direct comparison between the footprints calculated by the MPAS-CO2 adjoint model and the CT-L footprints.

Figure 5Comparison of footprints calculated by MPAS-CO2 adjoint model and CT-L of an OCO-2 sounding (17.82 N, 77.88 W) (red crosses; a–d) along the ground track shown in Fig. 1 (blue color). The footprints are calculated at four different heights: 500, 2000, 4500, and 10 000 m, and 10 d backward in time.

Figure 5 shows the footprints of a point located south of Jamaica in the Caribbean Sea (17.82 N, 77.88 W) at four different height levels: 500, 2000, 4500, and 10 000 m. The top three rows of Fig. 5 show that both the MPAS-CO2 adjoint model and CT-L footprints largely extend eastward over the Atlantic Ocean, indicating transport from the surface due to the influence of the easterly trade winds. Additionally, the MPAS-CO2 adjoint model-calculated footprint includes a branch that crosses the Equator and extends southeastward to the Southern Hemisphere between 30 and 40 W longitude. This feature is not shown in the CT-L footprint due to its limited area domain. The bottom row of Fig. 5 shows that the footprints of 10 000 m calculated by both the MPAS-CO2 adjoint and CT-L feature a primarily counterclockwise extension, covering the Gulf of Mexico and Texas. Moreover, there is a second segment extending westward from Texas toward the west coast. Upon closer examination, we observe that the CT-L-calculated footprint has a lower magnitude than the MPAS-CO2 adjoint model at the 500, 2000, and 4500 m height levels but not at the 10 000 m height level. The distinct patterns in both systems' footprints at different height levels indicate significant differences in horizontal and vertical transport patterns.

Figure 6Same as Fig. 5, except for a different OCO-2 sounding location (36.8 N, 82.9 W).

Figure 7Same as Fig. 5, except for a different OCO-2 sounding location (56.96 N, 91.89 W).

Figure 6 shows the corresponding footprints for an OCO-2 XCO2 sounding location in eastern Kentucky (36.8 N, 82.9 W) for the same four height levels as shown in Fig. 5. The figure shows that the footprints of 500 and 2000 m extend predominantly northward, covering the Great Lakes region and part of the Canadian Shield. In comparison, the footprint for 4500 m is mostly directed to the west. Another notable difference is that the highest-magnitude portion of the 500 and 2000 m footprints are in close proximity to the sounding location, while the 4500 m footprint is not in proximity at all. These differences among the height levels are evident in both the MPAS-CO2 adjoint model and CT-L calculated footprints. Figure 7 shows the footprints of an OCO-2 sounding location on the southwest coast of Hudson Bay, Canada (56.96 N, 91.89 W), for the four height levels. Both the MPAS-CO2 adjoint and CT-L footprints for 500 and 2000 m are largely either close to or north of the sounding location, indicating that surface fluxes from these regions have a significant influence on the atmospheric CO2 at the two height levels. In comparison, the footprint of 4500 m is located more than 2000 km northwestward, mostly covering Alaska; the particles move that far in a horizontal manner in the time it takes them to advect and mix 4500 m in the vertical. These findings emphasize the significant impact of vertical mixing on the spatial distribution of footprint at different altitudes, highlighting the unique patterns of horizontal and vertical transport in each case.

Figure 8The footprint of the OCO-2 ground track shown in Fig. 1 (with blue color) calculated by the MPAS-CO2 adjoint model (a–d) and by CT-L (e–h). The footprints are calculated by placing the adjoint forcing (for the MPAS-CO2 adjoint model) or releasing particles (for CT-L) at four different height levels above the ground: 500, 2000, 4500, and 10 000 m. The footprints are computed for 10 d backward in time.

In additional MPAS-CO2 adjoint model runs, we quantitatively compare the footprints of the entire OCO-2 track at each of the 14 height levels with the CT-L footprints. This comparison is conducted by performing a single MPAS-CO2 adjoint model run for each height level to calculate the footprint at the end of the 10 d backward-in-time integration. We then compare these resulting footprints with their CT-L counterparts. Figure 8 shows the comparison at four height levels of 500, 2000, 4500, and 10 000 m above the surface. At each height level, the value in the figure represents the average of the footprints of all the cells that are part of the OCO-2 track. The figure reveals that the footprints calculated by the two systems have similar spatial patterns within the limited-area domain of CT-L. However, it is important to note that a substantial portion of the footprints extends beyond the CT-L domain. For instance, the footprints of 2000 and 4500 m levels have significant coverage over Siberia in Russia, while the footprint of the 10 000 m level extends from the eastern Pacific Ocean to northeastern and western China, both of which are outside the CT-L model domain.

Figure 9Comparison of the OCO-2 ground track footprints from the MPAS-CO2 adjoint model and CT-L after 10 d backward-in-time integration. For each of the 14 height levels, the values of the footprints (units in ppm(µmolm-2s-1)-1) are extracted as the average value of 2× 3 boxes within the range of the CT-L spatial domain (10–80 N, 0–180 W). The solid line in each subfigure is the 1:1 line, and the dashed line is a linear fit with zero intercept. The correlation coefficient R2 of the linear fit is also labeled in each subfigure.


In order to compare the footprints from the two systems quantitatively, we aggregated the footprints onto a 2× 3 (lat × long) grid within the area covered by the CT-L model domain for each of the 14 height levels. Figure 9 shows the comparison for each of the 14. In the figure, the CT-L calculated footprints are on the x axis, and MPAS-CO2 adjoint model calculated footprints are on the y axis. The solid line in each subfigure of Fig. 9 is the 1:1 line, and the dashed line is a linear fit without intercept. The correlation coefficient R2 is labeled in each subfigure. The figure demonstrates that the agreement between the two systems is better for footprints at lower heights, particularly between 250 and 1500 m, with R2 all greater than 0.7. Footprints from the two systems agree to a much lesser degree at between 3500 and 14 000 m, where R2 is less than 0.5 in all cases. The linear fit lines (dashed lines) show that the MPAS-CO2 adjoint model calculated footprints are of greater magnitude in general than their CT-L counterparts at heights ranging from 50 to 1000 m. Between 1500 and 2500 m, the two sets of footprints are of similar magnitude on average, and at 3500 m and above, the CT-L footprints are of larger magnitude in general.

These differences in magnitude between the two systems could be attributed to various factors, including differences in model configurations, spatial resolution, and treatment of vertical mixing processes. Previous studies have shown that Lagrangian models, such as CT-L, can sometimes have different vertical mixing behavior compared to Eulerian models, especially at high altitudes (Karion et al.2019).

5.2 Influence of vertical distribution of OCO-2 soundings on footprints

In a final experiment, we use the MPAS-CO2 adjoint model to examine the impact of different vertical distributions on footprint calculation. Two adjoint model simulations were conducted for the OCO-2 orbit that crosses South America and North America between 17:36 UTC and 18:13 UTC on 23 August 2016 (the red color track in Fig. 1). Both simulations have the same adjoint forcing of 1 ppm XCO2 added to each MPAS-CO2 model cell along the orbital track at 18:00 UTC on 23 August 2016 and run backward in time for 30 d. The difference between the two simulations lies in the vertical distributions of the adjoint forcing. For the first simulation, we adopt profile 1, which is obtained by combining the XCO2 averaging kernel and pressure weight function (O'Dell et al.2018). In contrast, profile 2 prioritizes XCO2 information in the lower part of the troposphere (Fig. 10). The 20 pressure levels in the figure are interpolated to the MPAS-CO2 model's 55 vertical levels for the adjoint forcing placement. This experiment aims to highlight how these differences in vertical distribution impact the footprint calculation, leading to variations in flux estimation using variational assimilation. The results of this experiment will provide valuable insights to the importance of selecting appropriate vertical distribution when using the adjoint model for CO2 flux estimation.

Figure 10Two different profiles for vertically distributing a unit (ppm) of XCO2. Profile 1 is determined by OCO-2 XCO2 averaging kernel and pressure weight functions. Profile 2 is based on a redistribution of Profile 1 that gives more weight towards CO2 in the lower troposphere than in the upper part of the atmospheric column. The circles of profiles are on the 20 pressure levels of OCO-2 XCO2 pressure weight function. Both profiles integrate to unity.


Figure 11MPAS-CO2 adjoint model-calculated footprints (units in ppm(µmolm-2s-1)-1) obtained after 30 d of backward-in-time integration starting on 23 August 2016 at 18:00 UTC (the time of the OCO-2 measurement). The top panel (a) is obtained when using Profile 1 (Fig. 10) to vertically distribute 1 ppm of adjoint forcing. The middle panel (b) is the footprint using Profile 2. The bottom panel (c) is the difference in footprint between the two profiles.

The top two panels of Fig. 11 show the footprints resulting from MPAS-CO2 adjoint model simulations using the two distinct vertical distribution profiles for the adjoint forcing (Fig. 10). Although the two footprints may initially appear very similar, substantial differences become evident, as shown in the bottom panel in Fig. 11. Specifically, the footprint calculated using Profile 1 exhibits lower magnitudes compared to that obtained using Profile 2 in most extratropical regions in both the Northern Hemisphere and Southern Hemisphere. Conversely, over some of the tropical regions, particularly the tropical Pacific Ocean, the footprint calculated using Profile 1 shows slightly higher magnitudes than for Profile 2. Since the two adjoint model simulations have the same meteorology, these differences in the resulting footprints might be explained by how the convective transport of CO2 impacts the two distinctive vertical distribution profiles of the adjoint forcing. The prevalence of deep convection over the tropical Pacific Ocean can more effectively transport surface CO2 flux to the upper atmosphere than over the extratropics, where surface CO2 flux is more likely to be confined in the lower atmosphere. Thus, Profile 1's higher amount of adjoint forcing in the upper atmosphere results in its higher-magnitude footprint over the tropical Pacific Ocean but not over the extratropics, where its lower amount of adjoint forcing in the lower atmosphere leads to its lower-magnitude footprint. These findings underscore the critical importance of selecting an appropriate vertical distribution for the model–data difference when using an adjoint model during variational assimilation.

6 Conclusions

The MPAS-CO2 system consists of forward, TL, and adjoint models that are built upon the variable-resolution capability of the compressible non-hydrostatic MPAS-A model (Skamarock et al.2012). It promises to be a useful tool for carbon flux inverse modeling at the global and regional scales. The forward model of MPAS-CO2 is documented by Zheng et al. (2021). In this paper, we focus on the development of its tangent linear and adjoint models. Through rigorous testing, we have confirmed the correctness and accuracy of the newly developed MPAS-CO2 TL and adjoint models. A key challenge in developing the adjoint model was efficiently accessing meteorological variables during the model's backward-in-time integration. We have successfully implemented a strategy that combines recalculation and storage of meteorological variables. This approach significantly reduces the memory requirement, making the adjoint model feasible for long simulations, which are often necessary for CO2 inverse modeling.

The results of the sensitivity analysis using the newly developed MPAS-CO2 adjoint model provide valuable insights for designing CO2 data assimilation systems. The increasing homogeneity of the sensitivity to the initial atmosphere CO2 mixing ratio with longer integration length highlights the importance of selecting an appropriate assimilation window length. The comparison of the CO2 observation footprints between the MPAS-CO2 adjoint model and the NOAA CT-L system demonstrates good agreement, validating the accuracy of the adjoint model's footprint calculations. The comparison of OCO-2 XCO2 footprints reveals differences in sensitivity between the two systems at different altitudes. MPAS-CO2 adjoint model-calculated footprints tend to have higher magnitudes at low altitudes and lower magnitudes at high altitudes compared to CT-L. These differences in footprints could be caused by the differences in configuration, spatial resolution, and vertical mixing processes between the two model systems. Last, the sensitivity analysis using two different vertical distribution profiles for adjoint forcing highlights the importance of correctly mapping model–data difference in XCO2 to the transport model's vertical levels.

In addition to being a powerful tool for sensitivity analysis, the adjoint model plays a critical role in CO2 variational data assimilation (Bosman and Krol2023; Tian and Zou2021; Zheng et al.2018). Our future research efforts will focus on integrating the forward and adjoint models of MPAS-CO2 into such a system. This integration has the potential to bridge a significant gap by establishing an online Eulerian-transport-model-based global variational assimilation system for CO2 that targets high resolution in critical regions while at the same time avoiding the pitfalls associated with the lateral boundaries needed in regional domain inversions.

Code and data availability

The MPAS-CO2 forward, TL, and adjoint models v7.3 described in this paper can be downloaded from the CERN-based Zenodo archive at (Zheng2023a). This includes the model source code, instructions for compilation, and an example script for running models. Instructions for how to compile and run the models are provided in the package. The computation and plotting scripts used to produce the figures in this article can be downloaded from the CERN-based Zenodo archive at (Zheng2023b). CarbonTracker CO2 flux and posterior mixing ratio data can be obtained from the NOAA website at (Jacobson et al.2020). CT-L footprints can be obtained from the NOAA website at (NOAA Global Monitoring Laboratory2014).

Author contributions

TZ designed and developed the MPAS-CO2 TL and adjoint models. XT generated the 120–480 km VR mesh used in model simulations. TZ, SF, JS, and XT designed and carried out the model accuracy verification experiments. TZ, SF, DB, and MB designed the adjoint sensitivity analysis experiments. All authors contributed to writing the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


We thank the MPAS development team for making their code available to the public. We thank William Skamarock for his advice regarding the MPAS-A advection scheme. We thank Andrea Stenke, Peter Bosman, and the anonymous reviewer for their constructive and helpful reviews that greatly improved the paper. The model developments and simulations were carried out using Michigan State University High-Performance Computing Center (HPCC) and NASA's High-End Computing (HEC) facility. PNNL is operated by the U.S. Department of Energy by the Battelle Memorial Institute under contract no. DE-AC05-76RL01830. This is contribution 198 of the Central Michigan University Institute for Great Lakes Research.

Financial support

This research was supported by the National Aeronautics and Space Administration Carbon Monitoring System (grant no. 80HQTR21T0069).

Review statement

This paper was edited by Andrea Stenke and reviewed by Peter Bosman and one anonymous referee.


Agustí-Panareda, A., Diamantakis, M., Massart, S., Chevallier, F., Muñoz-Sabater, J., Barré, J., Curcoll, R., Engelen, R., Langerock, B., Law, R. M., Loh, Z., Morguí, J. A., Parrington, M., Peuch, V.-H., Ramonet, M., Roehl, C., Vermeulen, A. T., Warneke, T., and Wunch, D.: Modelling CO2 weather – why horizontal resolution matters, Atmos. Chem. Phys., 19, 7347–7376,, 2019. a

Baker, D. F., Doney, S. C., and Schimel, D. S.: Variational data assimilation for atmospheric CO2, Tellus B, 58, 359–365, 2006. a, b, c

Bosman, P. J. M. and Krol, M. C.: ICLASS 1.1, a variational Inverse modelling framework for the Chemistry Land-surface Atmosphere Soil Slab model: description, validation, and application, Geosci. Model Dev., 16, 47–74,, 2023. a, b

Byrne, B., Baker, D. F., Basu, S., Bertolacci, M., Bowman, K. W., Carroll, D., Chatterjee, A., Chevallier, F., Ciais, P., Cressie, N., Crisp, D., Crowell, S., Deng, F., Deng, Z., Deutscher, N. M., Dubey, M. K., Feng, S., García, O. E., Griffith, D. W. T., Herkommer, B., Hu, L., Jacobson, A. R., Janardanan, R., Jeong, S., Johnson, M. S., Jones, D. B. A., Kivi, R., Liu, J., Liu, Z., Maksyutov, S., Miller, J. B., Miller, S. M., Morino, I., Notholt, J., Oda, T., O'Dell, C. W., Oh, Y.-S., Ohyama, H., Patra, P. K., Peiro, H., Petri, C., Philip, S., Pollard, D. F., Poulter, B., Remaud, M., Schuh, A., Sha, M. K., Shiomi, K., Strong, K., Sweeney, C., Té, Y., Tian, H., Velazco, V. A., Vrekoussis, M., Warneke, T., Worden, J. R., Wunch, D., Yao, Y., Yun, J., Zammit-Mangion, A., and Zeng, N.: National CO2 budgets (2015–2020) inferred from atmospheric CO2 observations in support of the global stocktake, Earth Syst. Sci. Data, 15, 963–1004,, 2023. a

Courtier, P., Thepaut, J. N., and Hollingsworth, A.: A Strategy for Operational Implementation of 4d-Var, Using an Incremental Approach, Q. J. Roy. Meteor. Soc., 120, 1367–1387, b, 1994. a

Cui, Y. Y., Zhang, L., Jacobson, A. R., Johnson, M. S., Philip, S., Baker, D., Chevallier, F., Schuh, A. E., Liu, J., Crowell, S., Peiro, H. E., Deng, F., Basu, S., and Davis, K. J.: Evaluating Global Atmospheric Inversions of Terrestrial Net Ecosystem Exchange CO2 Over North America on Seasonal and Sub-Continental Scales, Geophys. Res. Lett., 49, e2022GL100147,, 2022. a

Errico, R. M.: What is an adjoint model?, B. Am. Meteorol. Soc., 78, 2577–2592, 1997. a, b, c

Errico, R. M. and Vukicevic, T.: Sensitivity analysis using an adjoint of the PSU-NCAR mesoseale model, Mon. Weather Rev., 120, 1644–1660, 1992. a

Feng, S., Lauvaux, T., Newman, S., Rao, P., Ahmadov, R., Deng, A., Díaz-Isaac, L. I., Duren, R. M., Fischer, M. L., Gerbig, C., Gurney, K. R., Huang, J., Jeong, S., Li, Z., Miller, C. E., O'Keeffe, D., Patarasuk, R., Sander, S. P., Song, Y., Wong, K. W., and Yung, Y. L.: Los Angeles megacity: a high-resolution land–atmosphere modelling system for urban CO2 emissions, Atmos. Chem. Phys., 16, 9019–9045,, 2016. a

Giering, R. and Kaminski, T.: Recipes for adjoint code construction, ACM T. Math. Software, 24, 437–474, 1998. a

Giering, R., Kaminski, T., Todling, R., Errico, R., Gelaro, R., and Winslow, N.: Tangent Linear and Adjoint Versions of NASA/GMAO's Fortran 90 Global Weather Forecast Model, in: Automatic Differentiation: Applications, Theory, and Implementations, edited by: Bücker, M., Corliss, G., Naumann, U., Hovland, P., and Norris, B., Springer Berlin Heidelberg, Berlin, Heidelberg, 275–284, 2006. a

Giles, M. B. and Pierce, N. A.: An introduction to the adjoint approach to design, Flow Turbul. Combust., 65, 393–415, 2000. a

Grell, G., Freitas, S. R., Stuefer, M., and Fast, J.: Inclusion of biomass burning in WRF-Chem: impact of wildfires on weather forecasts, Atmos. Chem. Phys., 11, 5289–5303,, 2011. a

Grell, G. A., Peckham, S. E., Schmitz, R., McKeen, S. A., Frost, G., Skamarock, W. C., and Eder, B.: Fully coupled online chemistry within the WRF model, Atmos. Environ., 39, 6957–6975, 2005. a

Guerrette, J. J. and Henze, D. K.: Development and application of the WRFPLUS-Chem online chemistry adjoint and WRFDA-Chem assimilation system, Geosci. Model Dev., 8, 1857–1876,, 2015. a

Hascoet, L. and Pascual, V.: The Tapenade Automatic Differentiation Tool: Principles, Model, and Specification, ACM T. Math. Software, 39, 20, 2013. a

Henze, D. K., Hakami, A., and Seinfeld, J. H.: Development of the adjoint of GEOS-Chem, Atmos. Chem. Phys., 7, 2413–2433,, 2007. a, b

Hoffmann, L., Günther, G., Li, D., Stein, O., Wu, X., Griessbach, S., Heng, Y., Konopka, P., Müller, R., Vogel, B., and Wright, J. S.: From ERA-Interim to ERA5: the considerable impact of ECMWF's next-generation reanalysis on Lagrangian transport simulations, Atmos. Chem. Phys., 19, 3097–3124,, 2019. a, b

Hong, S.-Y., Noh, Y., and Dudhia, J.: A new vertical diffusion package with an explicit treatment of entrainment processes, Mon. Weather Rev., 134, 2318–2341,, 2006. a

Hu, L., Andrews, A. E., Thoning, K. W., Sweeney, C., Miller, J. B., Michalak, A. M., Dlugokencky, E., Tans, P. P., Shiga, Y. P., Mountain, M., Nehrkorn, T., Montzka, S. A., McKain, K., Kofler, J., Trudeau, M., Michel, S. E., Biraud, S. C., Fischer, M. L., Worthy, D. E. J., Vaughn, B. H., White, J. W. C., Yadav, V., Basu, S., and van der Velde, I. R.: Enhanced North American carbon uptake associated with El Nino, Science Advances, 5,, 2019. a

Hurtt, G. C., Andrews, A., Bowman, K., Brown, M. E., Chatterjee, A., Escobar, V., Fatoyinbo, L., Griffith, P., Guy, M., Healey, S. P., Jacob, D. J., Kennedy, R., Lohrenz, S., McGroddy, M. E., Morales, V., Nehrkorn, T., Ott, L., Saatchi, S., Carlo, E. S. Serbin, S. P., and Tian, H.: The NASA carbon monitoring system phase 2 synthesis: scope, findings, gaps and recommended next steps, Environ. Res. Lett., 17, 063010,, 2022. a

Jacobson, A. R., Schuldt, K. N., Tans, P., Arlyn, A., Miller, J. B., Oda, T., Mund, J., Weir, B., Ott, L., Aalto, T., Abshire, J. B., Aikin, K., Aoki, S., Apadula, F., Arnold, S., Baier, B., Bartyzel, J., Beyersdorf, A., Biermann, T., Biraud, S. C., Boenisch, H., Brailsford, G., Brand, W. A., Chen, G., Chen, H., Chmura, L., Clark, S., Colomb, A., Commane, R., Conil, S., Couret, C., Cox, A., Cristofanelli, P., Cuevas, E., Curcoll, R., Daube, B., Davis, K. J., De Wekker, S., Coletta, J. D., Delmotte, M., DiGangi, E., DiGangi, J. P., Di Sarra, A. G., Dlugokencky, E., Elkins, J. W., Emmenegger, L., Fang, S., Fischer, M. L., Forster, G., Frumau, A., Galkowski, M., Gatti, L. V., Gehrlein, T., Gerbig, C., Gheusi, F., Gloor, E., Gomez-Trueba, V., Goto, D., Griffis, T., Hammer, S., Hanson, C., Haszpra, L., Hatakka, J., Heimann, M., Heliasz, M., Hensen, A., Hermansen, O., Hintsa, E., Holst, J., Ivakhov, V., Jaffe, D. A., Jordan, A., Joubert, W., Karion, A., Kawa, S. R., Kazan, V., Keeling, R. F., Keronen, P., Kneuer, T., Kolari, P., Komínková, K., Kort, E., Kozlova, E., Krummel, P., Kubistin, D., Labuschagne, C., Lam, D. H. Y., Lan, X., Langenfelds, R. L., Laurent, O., Laurila, T., Lauvaux, T., Lavric, J., Law, B. E., Lee, J., Lee, O. S. M., Lehner, I., Lehtinen, K., Leppert, R., Leskinen, A., Leuenberger, M., Levin, I., Levula, J., Lin, J., Lindauer, M., Loh, Z., Lopez, M., Luijkx, I. T., Lunder, C. R., Machida, T., Mammarella, I., Manca, G., Manning, A., Manning, A., Marek, M. V., Martin, M. Y., Matsueda, H., McKain, K., Meijer, H., Meinhardt, F., Merchant, L., Mihalopoulos, N., Miles, N. L., Miller, C. E., Mitchell, L., Mölder, M., Montzka, S., Moore, F., Moossen, H., Morgan, E., Morgui, J.-A., Morimoto, S., Müller-Williams, J., William Munger, J., Munro, D., Myhre, C. L., Nakaoka, S.-I., Necki, J., Newman, S., Nichol, S., Niwa, Y., Obersteiner, F., O'Doherty, S., Paplawsky, B., Peischl, J., Peltola, O., Piacentino, S., Pichon, J.-M., Pickers, P., Piper, S., Pitt, J., Plass-Dülmer, C., Platt, S. M., Prinzivalli, S., Ramonet, M., Ramos, R., Reyes-Sanchez, E., Richardson, S. J., Riris, H., Rivas, P. P., Ryerson, T., Saito, K., Sargent, M., Sasakawa, M., Scheeren, B., Schuck, T., Schumacher, M., Seifert, T., Sha, M. K., Shepson, P., Shook, M., Sloop, C. D., Smith, P., Stanley, K., Steinbacher, M., Stephens, B., Sweeney, C., Thoning, K., Timas, H., Torn, M., Tørseth, K., Trisolino, P., Turnbull, J., Van Den Bulk, P., Van Dinther, D., Vermeulen, A., Viner, B., Vitkova, G., Walker, S., Watson, A., Wofsy, S. C., Worsey, J. Worthy, D., Young, D., Zaehle, S., Zahn, A., and Zimnoch, M.: CarbonTracker CT2022, NOAA Global Monitoring Laboratory [data set], NOAA Global Monitoring Laboratory,, 2023. a, b

Kain, J. S.: The Kain–Fritsch Convective Parameterization: An Update, J. Appl. Meteorol., 43, 170–181,<0170:TKCPAU>2.0.CO;2, 2004. a, b

Karion, A., Lauvaux, T., Lopez Coto, I., Sweeney, C., Mueller, K., Gourdji, S., Angevine, W., Barkley, Z., Deng, A., Andrews, A., Stein, A., and Whetstone, J.: Intercomparison of atmospheric trace gas dispersion models: Barnett Shale case study, Atmos. Chem. Phys., 19, 2561–2576,, 2019. a

Kawa, S. R., Erickson, D. J., Pawson, S., and Zhu, Z.: Global CO2 transport simulations using meteorological data from the NASA data assimilation system, J. Geophys. Res.-Atmos., 109, 2493–2525,<2493:RATCPP>2.0.CO;2, 2004. a

Kopacz, M., Jacob, D. J., Henze, D. K., Heald, C. L., Streets, D. G., and Zhang, Q.: Comparison of adjoint and analytical Bayesian inversion methods for constraining Asian sources of carbon monoxide using satellite (MOPITT) measurements of CO columns, J. Geophys. Res.-Atmos., 114, d04305,, 2009. a

Krol, M., Houweling, S., Bregman, B., van den Broek, M., Segers, A., van Velthoven, P., Peters, W., Dentener, F., and Bergamaschi, P.: The two-way nested global chemistry-transport zoom model TM5: algorithm and applications, Atmos. Chem. Phys., 5, 417–432,, 2005. a

Lauvaux, T., Schuh, A. E., Uliasz, M., Richardson, S., Miles, N., Andrews, A. E., Sweeney, C., Diaz, L. I., Martins, D., Shepson, P. B., and Davis, K. J.: Constraining the CO2 budget of the corn belt: exploring uncertainties from the assumptions in a mesoscale inverse system, Atmos. Chem. Phys., 12, 337–354,, 2012. a

Lin, J. C., Gerbig, C., Wofsy, S. C., Andrews, A. E., Daube, B. C., Davis, K. J., and Grainger, C. A.: A near-field tool for simulating the upstream influence of atmospheric observations: The Stochastic Time-Inverted Lagrangian Transport (STILT) model, J. Geophys. Res.-Atmos., 108, 4493,, 2003. a, b

Meirink, J. F., Eskes, H. J., and Goede, A. P. H.: Sensitivity analysis of methane emissions derived from SCIAMACHY observations through inverse modelling, Atmos. Chem. Phys., 6, 1275–1292,, 2006. a

Meirink, J. F., Bergamaschi, P., Frankenberg, C., d'Amelio, M. T. S., Dlugokencky, E. J., Gatti, L. V., Houweling, S., Miller, J. B., Roeckmann, T., Villani, M. G., and Krol, M. C.: Four-dimensional variational data assimilation for inverse modeling of atmospheric methane emissions: Analysis of SCIAMACHY observations, J. Geophys. Res.-Atmos., 113, d17301,, 2008. a

NOAA Global Monitoring Laboratory: CarbonTraker-Lagrange Footprints, NOAA Global Monitoring Laboratory [data set],, last access: 20 February 2024. a

O'Dell, C. W., Eldering, A., Wennberg, P. O., Crisp, D., Gunson, M. R., Fisher, B., Frankenberg, C., Kiel, M., Lindqvist, H., Mandrake, L., Merrelli, A., Natraj, V., Nelson, R. R., Osterman, G. B., Payne, V. H., Taylor, T. E., Wunch, D., Drouin, B. J., Oyafuso, F., Chang, A., McDuffie, J., Smyth, M., Baker, D. F., Basu, S., Chevallier, F., Crowell, S. M. R., Feng, L., Palmer, P. I., Dubey, M., García, O. E., Griffith, D. W. T., Hase, F., Iraci, L. T., Kivi, R., Morino, I., Notholt, J., Ohyama, H., Petri, C., Roehl, C. M., Sha, M. K., Strong, K., Sussmann, R., Te, Y., Uchino, O., and Velazco, V. A.: Improved retrievals of carbon dioxide from Orbiting Carbon Observatory-2 with the version 8 ACOS algorithm, Atmos. Meas. Tech., 11, 6539–6576,, 2018. a

Pillai, D., Gerbig, C., Kretschmer, R., Beck, V., Karstens, U., Neininger, B., and Heimann, M.: Comparing Lagrangian and Eulerian models for CO2 transport – a step towards Bayesian inverse modeling using WRF/STILT-VPRM, Atmos. Chem. Phys., 12, 8979–8991,, 2012. a

Rayner, P. J., Michalak, A. M., and Chevallier, F.: Fundamentals of data assimilation applied to biogeochemistry, Atmos. Chem. Phys., 19, 13911–13932,, 2019. a

Schuh, A. E., Jacobson, A. R., Basu, S., Weir, B., Baker, D., Bowman, K., Chevallier, F., Crowell, S., Davis, K. J., Deng, F., Denning, S., Feng, L., Jones, D., Liu, J., and Palmer, P. I.: Quantifying the impact of atmospheric transport uncertainty on CO2 surface flux estimates, Global Biogeochem. Cy., 33, 484–500, 2019. a

Schuh, A. E., Otte, M., Lauvaux, T., and Oda, T.: Far-field biogenic and anthropogenic emissions as a dominant source of variability in local urban carbon budgets: A global high-resolution model study with implications for satellite remote sensing, Remote Sens. Environ., 262, 112473,, 2021. a, b

Schuh, A. E., Byrne, B., Jacobson, A. R., Crowell, S. M., Deng, F., Baker, D. F., Johnson, M. S., Philip, S., and Weir, B.: On the role of atmospheric model transport uncertainty in estimating the Chinese land carbon sink, Nature, 603, E13–E14, 2022. a

Skamarock, W., Klemp, J., Dudhia, J., Gill, D., Barker, D., Duda, M., Huang, X., Wang, W., and Powers, J.: A description of the Advanced Research WRF version 3, NCAR Tech Note NCAR/TN-475+STR, University Corporation for Atmospheric Research, 2008. a

Skamarock, W. C. and Gassmann, A.: Conservative Transport Schemes for Spherical Geodesic Grids: High-Order Flux Operators for ODE-Based Time Integration, Mon. Weather Rev., 139, 2962–2975,, 2011. a

Skamarock, W. C., Klemp, J. B., Duda, M. G., Fowler, L. D., Park, S.-H., and Ringler, T. D.: A Multiscale Nonhydrostatic Atmospheric Model Using Centroidal Voronoi Tesselations and C-Grid Staggering, Mon. Weather Rev., 140, 3090–3105,, 2012. a, b, c

Stohl, A., Forster, C., Frank, A., Seibert, P., and Wotawa, G.: Technical note: The Lagrangian particle dispersion model FLEXPART version 6.2, Atmos. Chem. Phys., 5, 2461–2474,, 2005. a

Tian, X. and Zou, X.: Development of the tangent linear and adjoint models of the MPAS-Atmosphere dynamic core and applications in adjoint relative sensitivity studies, Tellus A, 72, 1–17,, 2020. a, b, c, d

Tian, X. and Zou, X.: Validation of a prototype global 4D-Var data assimilation system for the MPAS-atmosphere model, Mon. Weather Rev., 149, 2803–2817, 2021. a, b, c, d

Walko, R. L. and Avissar, R.: The Ocean-Land-Atmosphere Model (OLAM). Part I: Shallow-Water Tests, Mon. Weather Rev., 136, 4033–4044,, 2008. a

Wang, H., Skamarock, W. C., and Feingold, G.: Evaluation of scalar advection schemes in the Advanced Research WRF model using large-eddy simulations of aerosol–cloud interactions, Mon. Weather Rev., 137, 2547–2558, 2009. a

Zheng, T.: The forward, tangent linear, and adjoint models of MPAS-CO2 V7.3 global online chemical transport model system (v7.3), Zenodo [code],, 2023a. a

Zheng, T.: Computation and plotting scripts for the development of the tangent linear and adjoint models of the global online chemical transport model MPAS-CO2 v7.3, Zenodo [code],, 2023b. a

Zheng, T., French, N. H. F., and Baxter, M.: Development of the WRF-CO2 4D-Var assimilation system v1.0, Geosci. Model Dev., 11, 1725–1752,, 2018. a, b, c, d, e, f, g

Zheng, T., Nassar, R., and Baxter, M.: Estimating power plant CO2 emission using OCO-2 XCO2 and high resolution WRF-Chem simulations, Environ. Res. Lett., 14, 085001,, 2019. a

Zheng, T., Feng, S., Davis, K. J., Pal, S., and Morguí, J.-A.: Development and evaluation of CO2 transport in MPAS-A v6.3, Geosci. Model Dev., 14, 3037–3066,, 2021. a, b, c, d, e, f, g, h

Zou, X., Vandenberghe, F., Pondeca, M., and Kuo, Y.-H.: Introduction to adjoint techniques and the MM5 adjoint modeling system, NCAR Technical note, University Corporation for Atmospheric Research, 1997. a

Short summary
The tangent linear and adjoint models have been successfully implemented in the MPAS-CO2 system, which has undergone rigorous accuracy testing. This development lays the groundwork for a global carbon flux data assimilation system, which offers the flexibility of high-resolution focus on specific areas, while maintaining a coarser resolution elsewhere. This approach significantly reduces computational costs and is thus perfectly suited for future CO2 geostationery and imager satellites.