2.1 Cost function configuration

WRF-CO2 4D-Var is designed to optimize the CO₂ flux by assimilating CO₂ observations into an atmospheric chemistry transport model. The CO₂ flux is optimized through use of a linear scaling factor:

where $\stackrel{\mathrm{̃}}{E}$ is the CO₂ flux read from flux files, k_co2 is the flux scaling factor, and E is the effective CO₂ flux. It is the effective flux E that is used in WRF-Chem's emission driver to update CO₂ mixing ratio (q_co2). The flux scaling factor k_co2, its tangent linear variable g_k_co2, and its adjoint variable a_k_co2 are used in calculating model sensitivity and minimizing the cost function defined in Eq. (2).

The cost function J(x) of WRF-CO2 4D-Var follows the Bayes framework that is widely used in atmospheric chemistry and numerical weather prediction (NWP) data assimilations:

where the background cost function J_b(x) is defined as

and the observation cost function J_o(x) is defined as

In Eqs. (3)–(4), B is the background error covariance matrix, R is the observation error covariance matrix, M is the transport model, H is the observational operator, y_k is the observation vector, and x^b is the prior state vector. The superscript n indicates that xⁿ is the optimized state vector at the nth iteration. For a full list of variables used in this paper, please refer to Table 1. Throughout the paper, boldface lowercase characters represent vectors and boldface uppercase characters represent matrices.

Like other data assimilation systems, WRF-CO2 4D-Var is an optimization scheme. Its state vector x consists of the flux scaling factors k_co2 and lateral boundary condition scaling factors. The summation K in Eq. (4) indicates the entire assimilation time period is evenly split into K observation windows during which observational data are ingested into the assimilation system. Details about how observations are ingested in the variational optimization schemes are given in Sect. 2.2 and 2.3.

Table 1A list of symbols used in this article.

Download Print Version | Download XLSX

Two optimization schemes are implemented in WRF-CO2 4D-Var to minimize the cost function. The first scheme uses the L-BFGS-B algorithm (Byrd et al., 1995) and the second uses the Lanczos-CG (Lanczos, 1950) minimization algorithm. Both schemes are iterative processes, and they call on WRF-CO2 4D-Var model components (the forward, tangent linear, and adjoint models) to calculate the model sensitivity $\partial {\mathit{q}}_{\mathrm{co}\mathrm{2}}/\partial {\mathit{k}}_{\mathrm{co}\mathrm{2}}$ between the iterations. The two optimization schemes are described in Sect. 2.2 and 2.3, respectively, and the three model components are described in Sect. 2.4.

2.2 L-BFGS-B optimization

L-BFGS-B (Byrd et al., 1995) is a quasi-Newton method for nonlinear optimization with bound constraints. L-BFGS-B has been used in a number of atmospheric chemistry inverse modeling systems, including the GEOS-Chem adjoint model system (Henze et al., 2007) and the TM5 4D-Var system (Meirink et al., 2008). The diagram in Fig. 1 demonstrates the steps involved in the L-BFGS-B-based optimization scheme. The scheme is an iterative process which searches for the optimized k_co2 by minimizing the cost function defined in Eq. (2)–(4). Between its iterations, the minimization algorithm L-BFGS-B requires the values of the cost function and its gradient, which are supplied by the forward model and the adjoint model as indicated in Fig. 1.

The calculation of the cost function is carried out based on Eqs. (2)–(4). Starting with the prior estimate of k_co2, the forward model run generates the CO₂ mixing ratio q_co2, which is transformed from the WRF model space to the observation space by the forward observation operator H. This results in the H(M(xⁿ)) term in Eq. (4), which is then paired with the observation vector y_k to calculate the innovation vector ${\mathbf{d}}_k=H\left(M\right({\mathit{x}}^n\left)\right)-{\mathit{y}}_k$. Next, the innovation vector and observation error covariance R are used to calculate the observation cost function J_o(x) as expressed in Eq. (4). Finally, the background cost function J_b(x) is calculated according to Eq. (3) and combined with the observation cost function J_o(x) to form the total cost function J(x) according to Eq. (2).

L-BFGS-B requires the values of the cost function J(x) and its gradient ∇ J(x) in searching for the optimized k_co2. The gradient is calculated using Eq. (5).

Figure 1Diagram of L-BFGS-B-based optimization implemented for WRF-CO2 4D-Var.

Download

The first term on the right-hand side of Eq. (5) is the observation gradient and the second is the background gradient. The observation gradient is calculated in two steps. (1) The innovation vector is scaled by R⁻¹ and transformed to the WRF model space by the adjoint observation operator, resulting in ${\stackrel{\mathrm{̃}}{H}}^T{\mathbf{R}}^{-\mathrm{1}}\left(H\right(M\left({\mathit{x}}^n\right))-{\mathit{y}}_k)$, which is the adjoint forcing. (2) The adjoint forcing is then ingested by the WRF-CO2 adjoint model during its backward (in time) integration, which yields the observation gradient. Supplied with the values of the cost function and gradient, the L-BFGS-B algorithm finds a new value of k_co2, which is used for the next iteration. The iterative optimization process continues until a given convergence criterion is met. The L-BFGS-B-based optimization in WRF-CO2 4D-Var is implemented based on the Fortran code of Algorithm 788 version Lbfgsb.2.1 (Zhu et al., 1997). Version Lbfgsb.3.0 (Luis Morales and Nocedal, 2011) will be implemented in the next model update.

2.3 Incremental optimization

The second optimization scheme implemented for WRF-CO2 4D-Var is the incremental approach commonly used in NWP data assimilation systems, including ECMWF 4D-Var (Rabier et al., 2000) and WRFDA (Barker et al., 2012). A major difference between the L-BFGS-B-based optimization and the incremental optimization is that the former optimizes for the state vector while the latter optimizes for the state vector analysis increment. The incremental assimilation scheme uses a linear approximation to transform the observation cost function from what is defined in Eq. (4) to Eq. (6):

Compared to Eq. (4), Eq. (6) approximates the innovation vector by a sum of two parts. The first part, $H\left(M\right({\mathit{x}}^{n-\mathrm{1}}\left)\right)-{\mathit{y}}_k$, is the innovation vector from the previous iteration. The second part, $\stackrel{\mathrm{̃}}{H}\left(\stackrel{\mathrm{̃}}{M}\right({\mathit{x}}^n-{\mathit{x}}^{n-\mathrm{1}}\left)\right)$, is the state vector analysis increment $({\mathit{x}}^n-{\mathit{x}}^{n-\mathrm{1}})$ transformed by the tangent linear model $\stackrel{\mathrm{̃}}{M}$ and tangent linear observation operator $\stackrel{\mathrm{̃}}{H}$. With the linear approximation of the cost function the gradient is calculated by

In WRF-CO2 4D-Var, the incremental optimization is implemented as a double loop in which the outer loop calculates the first and second items on the right-hand side of Eq. (7), while the inner loop calculates the third and fourth items. The superscript n−1 indicates that xⁿ⁻¹ is the optimized state vector in the last outer loop, and superscript n indicates that xⁿ is the optimized state vector in the inner loop. The outer loop first calls the forward model M and adjoint model ${\stackrel{\mathrm{̃}}{M}}^T$ to calculate ${\stackrel{\mathrm{̃}}{M}}^T{\stackrel{\mathrm{̃}}{H}}^T{\mathbf{R}}^{-\mathrm{1}}\left(H\right(M\left({\mathit{x}}^{n-\mathrm{1}}\right)-{\mathit{y}}_k\left)\right)$ and $B^{-\mathrm{1}}({\mathit{x}}^{n-\mathrm{1}}-{\mathit{x}}^b)$, which remain unchanged during the subsequent inner loop calculation. The analysis increment $({\mathit{x}}^n-{\mathit{x}}^{n-\mathrm{1}})$ is optimized in the inner loop, which calls the tangent linear and adjoint models to calculate the third and fourth items of Eq. (7). Inner loop calculation is carried out by Lanczos-CG (Lanczos, 1950), which can optionally estimate eigenvalues of the cost function Hessian matrix (∇² J(x)). The diagram in Fig. 2 shows the structure of the Lanczos-CG-based incremental optimization implemented in WRF-CO2 4D-Var.

2.4 Forward, tangent linear, and adjoint models

WRFPLUS consists of three model components: the WRF model, its tangent linear model, and its adjoint model (Barker et al., 2012; Huang et al., 2009). The three models are used by WRFDA to optimize the initial meteorological condition in order to improve numerical weather prediction. Unlike WRFDA, WRF-CO2 4D-Var is designed to optimize the CO₂ flux, instead of the meteorological initial and boundary conditions. This difference means the physical and dynamical processes involved in the atmospheric CO₂ transport are needed in WRF-CO2 4D-Var's model components. To include these processes, the chemistry module was added to the forward model. The chemistry module includes chemistry, deposition, photolysis, advection, diffusion, and convective transport of chemistry species (Grell et al., 2005). These processes are included in different modules of WRF-Chem: ARW (Advanced Research WRF) dynamical core, physics driver, and chemistry driver. The GHG (greenhouse gas) tracer option was removed and CO₂ is treated as an inert tracer. In the emission driver, the CarbonTracker 2016 version (Peters et al., 2007) replaces the online biogenic CO₂ model Vegetation Photosynthesis and Respiration Model (VPRM) (Mahadevan et al., 2008). This change is made because WRF-CO2 4D-Var optimizes the CO₂ flux instead of online emission model parameters.

2.4.1 Variable dependence analysis

The tangent linear and adjoint models of WRFPLUS needed to be modified to include the physical and dynamical processes involved in the atmospheric transport of CO₂, so that they will be consistent with the forward model. A thorough variable dependence analysis was conducted and the results are summarized in Table 2, which groups WRF-Chem processes into three categories regarding CO₂ tracer transport. The first category includes the chemistry processes that do not apply to CO₂, including gas- and aqueous-phase chemistry, dry and wet deposition, and photolysis. These processes are simply excluded from the forward, tangent linear, and adjoint models in WRF-CO2 4D-Var.

Figure 2Diagram of Lanczos-CG-based incremental optimization implemented for WRF-CO2 4D-Var.

Download

The second category is comprised of the physical parameterizations that do not provide CO₂ tendency but provide meteorological tendency. This category includes radiation, surface, cumulus, and microphysics parameterizations. While the full physics schemes of surface, cumulus, planetary boundary layer (PBL), and microphysics are used in the forward model of WRFPLUS, simplified versions of these schemes are used in its tangent linear and adjoint models. In addition, WRFPLUS uses full radiation schemes (longwave and shortwave) in its forward model, but it excludes radiation schemes from its tangent linear and adjoint models. The difference in the physical parameterizations between the forward model and tangent linear/adjoint models in a 4D-Var system is a source of linearization error. For instance, Tremolet (2004) found linearization error in ECMWF 4D-Var larger than expected and recommended more accurate linear physics for higher resolution 4D-Var systems. Because WRF-CO2 4D-Var ignores the impacts of CO₂ mixing ratio variation on the meteorological fields, no tangent linear and adjoint variables for meteorological fields are needed in its tangent linear and adjoint models. Since this second category of processes is not directly involved in CO₂ transport, there is no need for their tangent linear and adjoint procedures in WRF-CO2 4D-Var. In WRFPLUS's tangent linear model, the tangent linear code of the simplified versions of the cumulus, surface, and microphysics schemes was removed and replaced with the code for the full schemes as used in the forward model. In WRFPLUS's adjoint model, the forward sweep updates the state variables and local variables just as in the forward model, but it also stores these variables' values for the subsequent backward sweep, which updates the adjoint variables of the state variables. The simplified versions of the cumulus, surface, and microphysics schemes used in the forward sweep of WRFPLUS's adjoint model were removed and replaced with the full schemes used in the forward model. Since these processes do not directly modify CO₂ mixing ratio, their corresponding adjoint code was removed from the backward sweep of the adjoint model, as indicted by the “X” in Table 2.

Table 2Summary of variable dependence analysis for developing WRF-CO2 4D-Var component models on top of WRFPLUS. In the table, an “F” means a full physics scheme is used in the forward model, tangent linear model, or the forward sweep of the adjoint model. An “X” means a process is not needed for CO₂ treatment. A “Dev” means a process does not exist in WRFPLUS and has been developed for WRF-CO2 4D-Var. An “Add” means a process for CO₂ is simply added using the existing WRFPLUS code for other tracers.

Download Print Version | Download XLSX

The third category includes advection, diffusion, emission, and turbulence mixing in the PBL, along with convective transport of CO₂. Because these processes directly modify CO₂ mixing ratio, their tangent linear code and adjoint code are needed for WRF-CO2 4D-Var. The modifications made for advection and diffusion are described in Sect. 2.4.2, and those for emission, turbulent mixing in the PBL, and convective transport of CO₂ are detailed in Sect 2.4.3.

3.1 Model setup

WRF-CO2 4D-Var is set up with a domain covering the continental United States with 48 km × 48 km grid spacing and 50 vertical levels (Fig. 3). The domain dimension is 110 points in the east–west and 66 points in the north–south directions. Model configuration includes the Rapid Radiative Transfer Model (RRTM) longwave radiation (Mlawer et al., 1997), Goddard shortwave radiation (Chou and Suarez, 1999), Pleim surface layer (Pleim, 2006), Pleim–Xiu land surface model (Pleim and Xiu, 2003), ACM2 PBL (Pleim, 2007), Grell–Freitas cumulus (Grell and Freitas, 2014), and Thompson microphysics (Thompson et al., 2008). Positive–definite transport is applied to the transport of scalars and CO₂.

CO₂ fluxes used for the simulations are from the CarbonTracker 2016 version (hereafter CT2016) (Peters et al., 2007). These fluxes are the optimized surface fluxes at a 3 h interval and at 1×1 degree spatial resolution. The four individual CO₂ fluxes (biosphere, fossil fuel, fire, and ocean) are spatially interpolated to the WRF grid and saved in chemistry input files. In the following sensitivity tests and inverse experiments, the emission scaling factor k_co2 is applied only to the biosphere flux. Daily mean biospheric fluxes are calculated as the arithmetic mean of the 3-hourly CT2016 fluxes at each surface grid cell, and the scaling factor k_co2 is applied as in Eq. (1). The daily mean biospheric flux used for the 24 h simulation is shown in Fig. 4. The model configuration and emission data used are summarized in Table 3. Although the daily mean biospheric flux was used for the forward and inverse modeling in this paper, the WRF-CO2 4D-Var implementation allows flexibility in configuring the prior fluxes. For instance, fluxes from respiration and photosynthesis can be estimated independently and at higher temporal resolution (Gourdji et al., 2012). When using these options with real observations, the balance between the degrees of freedom in the state vector and the constraints provided by the observations needs to be carefully considered.

Figure 3WRF 4D-Var simulation domain covering the continental United States with 48 km × 48 km grid spacing. The domain boundary is marked by the bold dark outline. Grid cells used for evaluating sensitivities are marked: red triangles are the 20 CO₂ tower sites used as receptor locations; blue stars are source locations. While receptors are placed at the first, fifth, and tenth vertical levels at each site; all sources are at the first level only.

Download

Figure 4Daily mean CarbonTracker biosphere CO₂ flux, calculated as the arithmetic mean of the 3-hourly flux between 2 June 2011 at 00:00:00 UTC and 3 June 2011 at 00:00:00 UTC.

Download

Model simulations span 24 h from 00:00 UTC on 2 June to 00:00 UTC on 3 June 2011. Meteorological initial and lateral boundary conditions are prepared using the National Centers for Environmental Prediction (NCEP) Climate Forecast System version 2 (CFSv2) 1×1 degree 6-hourly products (Saha et al., 2014). CO₂ initial and lateral boundary conditions are from the CT2016 global 3×2^∘ CO₂ mole fraction. A method similar to PREP-CHEM-SRC (Freitas et al., 2010) was used to horizontally and vertically interpolate CT2016 mole fraction data to the WRF grid.

Table 3WRF-CO2 4D-Var model configuration and CO₂ flux used in sensitivity and inverse modeling tests.

Download Print Version | Download XLSX

Figure 5Sea level pressure (hPa) and horizontal wind (m s⁻¹) at the model's lowest vertical level plotted at 6 h intervals during the 24 h simulation starting on 2 June 2011 at 00:00 UTC.

Download

First, the forward model (WRF-Chem) was run for 24 h with the CO₂ emission as described in the last section. Trajectory files that contain model state variables including both meteorology and CO₂ mixing ratio are saved at model dynamical time step intervals (120 s). These files are required for the subsequent tangent linear and adjoint model runs. Figure 5 shows the instantaneous values of sea level pressure (SLP) and horizontal wind at the model's lowest vertical level at every 6 h. The figure shows that a high pressure system was located off the west coast, causing a northerly surface wind off southern California, and a westerly wind for most of the Pacific Northwest. A low pressure system intensified over Montana and North Dakota during the 24 h, causing a strong southerly wind over the Midwest. In the northeast, as a low pressure system moved eastward out of the domain, the surface wind shifted from southwesterly to westerly.

In the model setup, the initial and boundary meteorological conditions are generated by downscaling the CFSv2 data. Downscaling coarse-resolution global reanalysis data could lead to poor WRF performance. Although this potential problem is not a concern for the present pseudo-observation-based inversion experiments, it must be properly treated in future applications with real observations. Error in the initial condition will lead to erroneous flux attribution, especially for inversions with a short assimilation window.

In order to be useful for applications which employ real observational data, WRF-CO2 4D-Var requires accurate simulations of the meteorological fields by the forward model. Because transport error can only be partially accounted for in the 4D-Var system through the observation error covariance matrix, it is imperative to minimize errors due to inaccurate simulation of meteorological processes as much as possible. Although the present paper uses pseudo-observation data (which have zero transport error by definition) in its inversion experiments, future applications with real observations will require vigorous evaluation of the model-simulated meteorology and associated transport error. In the following, the forward-model-simulated horizontal winds at the surface and 500 hPa constant pressure surface are evaluated using in situ measurements from weather stations and radiosondes.

For the surface level, WRF-simulated 10 m winds are compared against surface weather station measurements archived in the National Oceanic and Atmospheric Administration (NOAA) Integrated Surface Database (Smith et al., 2011). Hourly surface wind measurements from more than 2000 stations within the WRF domain are used for the evaluation. Comparisons of wind speed and wind direction are carried out at the top of each hour during the 24 h simulation period starting at 00:00 UTC on 2 June 2011. Excluding missing observations, this results in 31 745 valid data pairs, which are summarized in the histograms of Fig. 6. RMSE for the hourly wind speed is 2.16 m s⁻¹ and the mean difference in the hourly wind direction is 29.4^∘.

For the upper level, WRF-simulated 500 hPa horizontal winds were compared against radiosonde measurements from 90 stations obtained from the NOAA Earth System Research Laboratory (ESRL) radiosonde database (https://ruc.noaa.gov/raobs/, last access: 27 April 2018). Since most stations release balloons at 00:00 and 12:00 UTC, WRF winds were compared against the radiosonde measurements at a 12 h interval during the 24 h simulation period. The results are shown in Fig. 7: RMSEs of wind speed are 2.54, 4.0, and 5.11 m s⁻¹, on 2 June at 00:00 UTC, 2 June at 12:00 UTC, and 3 June at 00:00 UTC, respectively. Wind direction differences between WRF and radiosonde are 11.5, 16.4, and 19.1^∘ at the three times. Locations of the weather stations and radiosonde sites used in the evaluations can be found in the Supplement.

Figure 6Histograms of the 10 m wind speed difference (a) and wind direction difference (b) between the WRF simulation and surface meteorological station measurements.

Download

Figure 7Comparison of 500 hPa wind speed and wind direction between WRF simulation and radiosonde measurements. Panels (a) and (b) are the comparison on 2 June 2011 at 00:00 UTC; panels (c) and (d) are on 2 June 2011 at 12:00 UTC; and panels (e) and (f) are on 3 June 2011 at 00:00 UTC. RMSE and relative error (RE) for wind speed and mean difference in wind direction are shown in each figure.

Download

The above-described evaluations using in situ measurements indicate that the meteorological simulation is of adequate accuracy for the pseudo-observation-based inverse modeling tests conducted in this paper. When the 4D-Var system is applied with real observations, the error and bias must be considered. In WRF 4D-Var's cost function configuration, the observation error matrix R is a combination of three error sources: measurement error, aggregation error, and transport model error. The uncertainty of the CO₂ measurements is about 0.05 %, while the transport and aggregation errors are typically an order of magnitude larger (Bruhwiler et al., 2005). For real observation applications, the variance and covariance in R need to represent the transport error. Furthermore, Fig. 6 shows that WRF-simulated 10 m wind speed is biased high, which is likely to result in bias in the simulated atmospheric CO₂ mixing ratio. Because Bayes inversion framework assumes unbiased observation error, it may be imperative to correct the error for inversions. One approach is to nudge the meteorology fields toward the observations. For instance, Gupta et al. (1997) found that nudging the model-simulated winds in the boundary layer to the radar wind profile observations substantially improved estimates of plume dispersion. An alternative approach is to use a combined 4D-Var inversion of meteorology and CO₂ fluxes. For instance, Bocquet et al. (2015) discussed data assimilation using coupled chemistry meteorology models (CCMMs). If the CO₂ impact on meteorology is not considered, the current implementation of WRF-CO2 4D-Var can be extended to a joint meteorology and CO₂ assimilation system. Since the adjoint code for meteorology has been developed and tested in WRFPLUS and WRFDA (Zhang et al., 2013; Barker et al., 2012), the major modification would be in the optimization schemes where the combined state vector of meteorology and CO₂ is optimized. It should be noted that in such a joint assimilation framework, optimization of meteorology is an initial condition problem, whereas the CO₂ flux optimization is a boundary condition problem (bottom and lateral boundaries).

3.2 Accuracy of tangent linear and adjoint sensitivities

Next, the accuracy of the newly developed tangent linear and adjoint models was evaluated by comparing their sensitivity calculations against finite difference sensitivity calculated by the forward model. Grid cells involved in the sensitivity calculations are shown in Fig. 3, in which the 35 blue stars are the source cells, and the 20 red triangles are 20 tower sites where the receptors are placed (Table 4). All the 35 sources are placed at the grid's bottom vertical level. Receptors are placed at the first, fifth, and tenth vertical levels at each of the 20 tower sites, resulting in 60 receptor cells.

A tangent linear model run for a grid cell will calculate the tangent linear sensitivity $\partial {\mathit{q}}_{\mathrm{co}\mathrm{2}}/\partial {\mathit{k}}_{\mathrm{co}\mathrm{2}}$, which approximates a column vector of the forward model's Jacobian matrix and quantifies the influence of the cell's flux change on CO₂ mixing ratio of its receptor cells downwind. In comparison, an adjoint model run for a grid cell will calculate adjoint sensitivity $\partial {\mathit{q}}_{\mathrm{co}\mathrm{2}}/\partial {\mathit{k}}_{\mathrm{co}\mathrm{2}}$, which approximates a row vector of the forward model's Jacobian matrix and quantifies the influence on the cell's CO₂ mixing ratio by its source cells upwind. Because k_co2 multiplies emission in Eq. (1), the magnitude of the sensitivity is determined by both the magnitude of emission and meteorological transport.

To calculate tangent linear sensitivity at a grid cell, g_k_co2 is set to unity at the cell and zero at all other cells at the start of a tangent linear model run. Upon completion, the values of g_q_co2 are the tangent linear sensitivities $\partial {\mathbf{q}}_{\mathrm{co}\mathrm{2}}/\partial {\mathit{k}}_{\mathrm{co}\mathrm{2}}$. To calculate adjoint sensitivity at a cell, an adjoint model run starts with a_q_co2 set to unity at the cell and zero at all others, and the values of a_k_co2 at the end of the simulation are the adjoint sensitivities. The adjoint model running in this mode is analogous to using a Lagrangian particle transport model in backward trajectory mode to compute the footprint of a receptor, as shown in Fig. 4 of Gerbig et al. (2008).

The tangent linear sensitivity is first compared against the finite difference sensitivity. After confirming the accuracy of the tangent linear model, the adjoint sensitivity is compared against the tangent linear sensitivity. Finite difference sensitivities are calculated using the two-sided formula (Eq. 8).

The magnitude of Δx used in Eq. (8) is determined by comparing the result from a range of different values. The finite sensitivities were calculated at the 35 sites using Δx set to 0.01, 0.1, and 1.0, and the results show that the magnitude of all differences is less than 10⁻¹⁰ (results not shown) because WRF-CO2 is largely linear. For all subsequent calculations, Δx=0.1 is used for Eq. (8).

Since both finite difference and tangent linear sensitivities form columns of the Jacobian matrix, their values can be compared cell by cell for all receptor cells for a given site. Figure 8 shows the comparison between the finite difference and tangent linear sensitivities at 9 of the 35 source cells. The dark straight lines in the figures are the 1 : 1 line. The maximum and minimum of the difference between finite difference and tangent linear sensitivities are given for each source cell. Results at the rest of the sources are similar (not shown). All differences are less than 10⁻¹⁰, confirming that the tangent linear model is accurate.

Figure 8Comparison between $\partial {\mathit{q}}_{\mathrm{co}\mathrm{2}}/\partial {\mathit{k}}_{\mathrm{co}\mathrm{2}}$ calculated by finite difference (x axis) and tangent linear models (y axis) for nine source cell locations (see Fig. 3 for source locations).

Download

The adjoint model is next evaluated by comparing adjoint sensitivities against the tangent linear sensitivities. Because finite difference sensitivities form columns of the Jacobian matrix while adjoint sensitivities form rows of the Jacobian matrix, they can only be compared at the intersections of the rows and columns of the Jacobian matrix, meaning there are 2160 (35×60) pairs of comparison. We organized these 2160 pairs into three groups based on the vertical levels a receptor is placed at and the result is shown in Fig. 9. The minimum and maximum value of the difference between tangent linear and adjoint sensitivities in all three groups are no greater than 10⁻⁶, which is comparable to the accuracy tests from other adjoint model developments (Meirink, 2008; Henze, 2007), indicating that the adjoint model has been correctly implemented.

Figure 9Comparison between $\partial {\mathit{q}}_{\mathrm{co}\mathrm{2}}/\partial {\mathit{k}}_{\mathrm{co}\mathrm{2}}$ calculated by the tangent linear (x axis) and adjoint models (y axis).

Download

3.3 Spatial patterns of adjoint sensitivities

Adjoint sensitivity q_co2∕k_co2 quantifies how q_co2 of a given receptor is impacted by the flux scaling factor of all surface cells. It is similar to the receptor footprint typically calculated using LPDM, such as Fig. 4 of Gerbig et al. (2008) and Fig. 1 of Alden et al. (2016). But q_co2∕k_co2 differs from footprint in that the former contains the combined impact of tracer transport and flux magnitude, while the latter is determined by tracer transport alone. The spatial patterns of the adjoint sensitivity were examined to discern the impacts of tracer transport. Figure 10 shows q_co2∕k_co2 of Centerville, Iowa (top row), and WLEF (tower at Park Fall), Wisconsin (bottom row). At each tower site, q_co2∕k_co2 of the receptor placed at the first and tenth vertical levels is plotted.

Figure 10Adjoint sensitivities calculated by the WRF-CO2 adjoint model. The top panels show adjoint sensitivity of receptors placed at the first (a) and tenth (d) vertical levels at Centerville, Iowa. The bottom panels show adjoint sensitivity of receptors placed at the first (c) and tenth (d) vertical levels at WLEF, Wisconsin. The black cross in each figure marks the corresponding tower site.

Download

Figure 11Footprints calculated using Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) backward trajectories and CarbonTracker biospheric fluxes for the tower sites at Centerville, Iowa, and WLEF, Wisconsin. The receptor locations are the same as in Fig. 10. Each HYSPLIT footprint is plotted in the same color scale as its counterpart in Fig. 10 for comparison.

Download

Figure 12The first-guess biosphere CO₂ fluxes used in the two inverse modeling experiments. The x axis is the true daily mean CarbonTracker biosphere CO₂ value (as shown in Fig. 4), and the y axis is the first guess (background value). The solid line in each figure is the 1 : 1 line.

Download

The adjoint sensitivities of the Centerville tower site indicate its q_co2 results primarily from surface flux located immediately south of the site. This pattern agrees with the fact that low-level wind during the simulation period is predominantly southerly, transporting tracers northward. There is also a marked difference in the adjoint sensitivity of the same tower site when the receptor is placed at a different height. The figure in the top left panel of Fig. 10 shows that the highest magnitude of q_co2∕k_co2 is closest to the tower itself, indicating a large impact from local fluxes. In comparison, when the receptor is placed at the tenth vertical level, the peak magnitude of its adjoint sensitivity is located further south of the tower site. Results from the WLEF site show the adjoint sensitivity is located to the southeast of the site, matching the southeasterly wind patterns around Wisconsin during the simulation period. There are also clear difference between the receptors at the different vertical levels. Results from other sites all show a similar pattern of impacts of transport and receptor placement height (not shown).

To provide a comparative view of the source–receptor relations, backward trajectories of particles released from the Centerville and WLEF sites were also calculated using the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model (Stein et al., 2015). WRF-CO2 forward-model-simulated meteorology saved at 1 h intervals was used to drive the HYSPLIT trajectory calculations. To compare with the adjoint model result, two sets of simulations were carried out for each of the two tower sites. For each simulation, 30 000 particles were released from the location of the corresponding WRF grid box used in the adjoint sensitivity calculations. The starting locations of the particles were randomly distributed within the grid box. The resulting backward trajectories were combined with the biospheric CO₂ flux to calculate the footprint for the receptor locations. The HYSPLIT footprints were calculated on the same grid as used in the WRF-CO2 simulations to facilitate the comparisons between the two models.

Figure 11 shows the HYSPLIT-calculated footprints for the Centerville and WLEF sites at the two different vertical levels. The four figures in Fig. 11 are the HYSPLIT counterparts of the adjoint sensitivity figures in Fig. 10. A comparison between Figs. 10 and 11 shows that the results from HYSPLIT and the WRF-CO2 adjoint model compare well spatially. For instance, for the receptor placed at the first vertical level at Centerville, Iowa (Figs. 10a and 11a), the footprints from both models are primarily located in Missouri and northwestern Arkansas. Based on the horizontal wind fields at the first level, these areas were upwind of the receptor location during the simulation period. Overall, the WRF-CO2 adjoint sensitivities contain larger surface areas compared to their HYSPLIT footprint counterparts. This difference is likely caused by the more diffusive nature of tracer transport in WRF-CO2: its finite difference scheme for tracer advection contains numerical diffusion, and it also includes an explicit horizontal diffusion term in the tracer transport (Skamarock et al., 2008). A further comparison at individual grid points reveals magnitude differences between the footprints from HYSPLIT and the WRF-CO2 adjoint model. This is mainly caused by the different treatments of turbulent vertical mixing by the two models. In WRF-CO2, the PBL and convective schemes parameterize tracer vertical mixing (see Sect. 2.4.3). For vertical mixing, HYSPLIT either uses the PBL heights calculated by WRF or it calculates PBL heights independently by analyzing temperature profiles. The footprints shown in Fig. 11 were simulated by HYSPLIT using PBL heights from the WRF-CO2 ACM2 PBL scheme. In a separate set of HYSPLIT simulations with PBL heights calculated from the temperature profiles, only minor differences are observed in the resulting footprints (not shown).

Figure 13Results of the inverse modeling experiment (Case 1). Panel (a) shows the reduction of the cost function, represented by J(xⁿ)∕J(x^b). Panel (b) shows the reduction of the gradient norm, represented by $\Vert \mathrm{\nabla }J\left({\mathit{x}}^n\right)\Vert /\Vert \mathrm{\nabla }J\left({\mathit{x}}^b\right)\Vert$. Panel (c) shows the reduction of biospheric CO₂ flux RMSE.

Download

Figure 14Comparison between the true and optimized CO₂ fluxes by Lanczos-CG (left column) and L-BFGS-B (right column) in the inverse modeling experiment (Case 1). The comparison and RMSE after the 1st, 5th, 10th, 30th, 50th, and 70th iterations are shown in the figure. All iterations of Lanczos-CG are from one outer loop.

Download

3.4 Inverse modeling tests

3.4.2 Inverse modeling results

The results from inverse modeling experiments with Case 1 prior are shown in Figs. 13 and 14. Figure 13 shows the iterative reduction of the cost function J(x), gradient norm $\Vert \mathrm{\nabla }J\left(\mathit{x}\right)\Vert$, and RMSE. The iteration number for Lanczos-CG is all from its inner loop, and only one outer loop is used. The figures show both L-BFGS-B and Lanczos-CG reduce the cost function monotonically. In about the first 10 iterations, the cost function reduction is more or less similar for the two optimization schemes, but Lanczos-CG starts to gradually outperform L-BFGS-B after. In gradient norm reduction, both schemes feature periodic oscillations embedded in the large-scale downward trend. By comparison, Lanczos-CG has a smaller magnitude oscillation and steeper downward trend than L-BFGS-B. It should be noted that while L-BFGS-B calculates cost function and its gradient in each iteration, Lanczos-CG only approximates these values in its inner loop. The cost function and gradient norm from Lanczos-CG shown in Fig. 13 are calculated by extra calls to the forward and adjoint models in each inner iteration, which doubles the computation cost and is not needed in actual inversion applications. Figure 13c shows that both optimization schemes reduce RMSE of daily biosphere flux monotonically, and Lanczos-CG achieves better reduction after about the first 10 iterations. Figure 14 shows the snapshots of the optimized daily mean biosphere flux (obtained as the product of the prior flux and the optimized scaling factor) at a selected set of iterations. These figures depict the iterative process of priors converging to the true solution.

Figure 15Same as Fig. 13 but for the inverse modeling experiment (Case 2).

Download

Figure 16Same as Fig. 14 but for the inverse modeling experiment (Case 2).

Download

The results of inverse modeling experiments using Case 2 prior are shown in Figs. 15 and 16. The reductions of J(x), $\Vert \mathrm{\nabla }J\left(\mathit{x}\right)\Vert$, and RMSE are similar to Case 1 in that Lanczos-CG substantially outperforms L-BFGS-G after about the first 10 iterations. Table 5 summarizes the results from all four inverse modeling experiments described above. It must be pointed out that these inverse modeling results are obtained from a highly unphysical setup, and they are not the expected level of performance (in terms of cost function and RMSE reduction) that would be obtained in a inversion with real observations.