Regional atmospheric CO

While rising atmospheric CO

Like ensemble Kalman smoothers, 4D-Var inversions do not require
precalculation of the influence function, but they do require the adjoint of
a CTM to calculate the observation cost function. 4D-Var inversions use a CTM
and prior CO

4D-Var systems have been widely used for CO

Regional 4D-Var inversion studies have been driven in part by the need to
resolve biosphere–atmosphere carbon exchange at smaller scales

In this paper, a regional CO

In the WRF system, CO

Compared with offline regional inversion systems, WRF-CO2 4D-Var has an
advantage provided by the close one-way coupling between meteorological
processes and chemistry transport. For example, adequately representing
CO

The remainder of this paper is organized as follows: Sect. 2 details the
implementation of the two variational optimization schemes for cost function
minimization, and the modification to the tangent linear and adjoint models.
Section 3 examines the accuracy of sensitivity calculated by the tangent
linear and adjoint models, and the system's effectiveness in inverse
modeling. Section 4 describes the treatment of CO

This section describes the WRF-CO2 4D-Var cost function configuration and the associated minimization schemes, followed by a description of the forward, tangent linear, and adjoint models.

WRF-CO2 4D-Var is designed to optimize the CO

The cost function

Like other data assimilation systems, WRF-CO2 4D-Var is an optimization
scheme. Its state vector

A list of symbols used in this article.

Two optimization schemes are implemented in WRF-CO2 4D-Var to minimize the
cost function. The first scheme uses the L-BFGS-B
algorithm

L-BFGS-B

The calculation of the cost function is carried out based on Eqs. (2)–(4).
Starting with the prior estimate of

L-BFGS-B requires the values of the cost function

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

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

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

Compared to Eq. (4), Eq. (6) approximates the innovation vector by a sum of
two parts. The first part,

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

WRFPLUS consists of three model components: the WRF model, its tangent linear
model, and its adjoint model

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

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

The second category is comprised of the physical parameterizations that do
not provide CO

Summary of variable dependence analysis for developing WRF-CO2 4D-Var
component models on top of WRFPLUS. In the table, an “

The third category includes advection, diffusion, emission, and turbulence
mixing in the PBL, along with convective transport of CO

WRF includes the advection and diffusion of inert tracers along with other
scalars in its ARW dynamical core. The tangent linear and adjoint code for
these processes has been implemented in WRFPLUS. It should be noted that the
variables for these inert tracers are part of WRF, instead of WRF-Chem.
WRF-Chem uses a separate array for its chemistry species. Since WRF-Chem is
used as the forward model of WRF-CO2 4D-Var, the CO

An accurate representation of vertical mixing is important for inversion
accuracy, because misrepresentation causes transport error, which manifests
itself in the innovation vector and causes error in posterior estimation

In WRF-Chem, vertical mixing of chemical species is treated in three separate parts: in the vertical diffusion (subgrid-scale filter) in the dynamical core, in the PBL scheme in the physics driver, and in the convective transport in the chemistry driver. The subgrid-scale filter in the dynamical core treats both horizontal and vertical diffusion, but vertical diffusion is turned off if a PBL scheme is used.

For PBL parameterization, the asymmetrical convective model version 2 (ACM2)

Convective transport of chemistry species in WRF-Chem is not treated by the
cumulus scheme in the physics driver but by a separate convective transport
module (module_ctrans_grell) in the chemistry driver

Because the ACM2 PBL and chemistry convective transport are not included in
WRFPLUS, their tangent linear and adjoint code was developed for WRF-CO2
4D-Var. First, the automatic differentiation tool Tapenade

This section presents an accuracy assessment of the newly developed WRF-CO2 4D-Var system. First, the simulation model setup is described; then, the sensitivity tests and inverse modeling experiments are presented.

WRF-CO2 4D-Var is set up with a domain covering the continental United States
with 48 km

CO

WRF 4D-Var simulation domain covering the continental United States
with 48 km

Daily mean CarbonTracker biosphere CO

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)

WRF-CO2 4D-Var model configuration and CO

Sea level pressure (hPa) and horizontal wind (m s

First, the forward model (WRF-Chem) was run for 24 h with the CO

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

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 (

Histograms of the 10 m wind speed difference

Comparison of 500 hPa wind speed and wind direction between WRF
simulation and radiosonde measurements. Panels

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

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

To calculate tangent linear sensitivity at a grid cell,

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

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

Comparison between

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 (

Comparison between

Adjoint sensitivity

Adjoint sensitivities calculated by the WRF-CO2 adjoint model.
The top panels show adjoint sensitivity of receptors placed at the
first

Footprints 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.

The first-guess biosphere CO

The adjoint sensitivities of the Centerville tower site indicate its

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

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

Results of the inverse modeling experiment (Case 1). Panel

Comparison between the true and optimized CO

The sensitivity tests in Sect. 3.2 have confirmed that the tangent linear and
adjoint models of WRF-CO2 4D-Var are correctly implemented. In this section,
inverse modeling tests are conducted to confirm that the two optimization
schemes described in Sect. 2.2 and 2.3 are also correctly implemented. The
inverse modeling tests here are designed following the approach used in

Run the WRF-CO2 forward model for 24 h using the daily mean biospheric CO

Generate pseudo-observations by saving the model-simulated atmospheric CO

Generate a set of first-guess biospheric CO

Set the background error to infinity (

Run the L-BFGS-B and incremental optimizations with the first-guess biospheric CO

Case 1: set flux scaling factor

Case 2: set flux scaling factor

A very simple error configuration (

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

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

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

The results of inverse modeling experiments using Case 2 prior are shown in
Figs. 15 and 16. The reductions of

The lateral tracer boundary condition is necessary to connect regional tracer
simulations to the global background tracer distribution

Summary of CO

Summary of inverse modeling experiment results. The reductions of
cost function

In the WRF-Chem dynamical core, chemistry mixing ratios are updated at each
time step by the advection and diffusion tendencies. Then, chemistry mixing
ratios at the lateral boundaries are updated with the chemistry boundary
condition using the flow-dependent method, which uses the horizontal wind
direction to determine whether the chemistry mixing ratio at a boundary grid
point should be updated by the lateral boundary. If the horizontal wind
direction indicates tracer inflow at a boundary grid, Eq. (9) will be applied
to the grid point:

Most code development for tracer lateral boundary conditions is in the
input_chem_data module of the chemistry directory, along with some
additional code modification to enable the lateral boundary condition
variables to be passed forward (

When applied with real observations, whether and how to aggregate lateral
boundary scaling factors is not trivial

WRF-CO2 4D-Var was developed as a data assimilation system designed to
constrain surface CO

WRF-CO2 4D-Var's tangent linear and adjoint models were tested by comparing their sensitivities' spatial patterns with the dominant wind patterns. The results make physical sense given the meteorological transport. The accuracy of tangent linear and adjoint models was evaluated by comparing their sensitivity against finite difference sensitivity calculated by the forward model. The results show that both tangent linear and adjoint sensitivities agree well with finite difference sensitivity. Finally, the system was tested in inverse modeling with pseudo-observations, and the results show that both optimization schemes successfully recovered the true values with reasonable accuracy and computation cost.

While Lanczos-CG performs better than L-BFGS-B in the inverse modeling tests,
it must be pointed out that the tests are very limited. Although a
comprehensive comparison between the two optimization schemes is beyond the
scope of the present paper, it is important to point out some of their
differences as implemented in WRF-CO2 4D-Var. First, the Lanczos-CG calls the
tangent linear model in each inner loop iteration, while L-BFGS-B calls the
forward model. For a tracer transport system like WRF-CO2 4D-Var, the tangent
linear model can skip some of the costly physics parameterizations, such as
the radiation scheme. This difference means that typically the tangent linear
model is faster than the forward model, and as a result Lanczos-CG runs
faster than L-BFGS-B. In our inversion modeling experiments (24 h simulation
with

Second, provided with the cost function and its gradient, each iteration of
L-BFGS-B calculates an updated state vector from its previous iteration. In
WRF-CO2 4D-Var, this calculation is carried out on only the root core and
broadcasted to the other process cores. In comparison, Lanczos-CG calculates
the state vector increment based on the cost function gradient alone (without
the need for

Another consideration for memory requirements is related to I/O time cost. WRFPLUS saves its entire trajectory in memory to avoid expensive I/O operations. This is not a practical solution for WRF-CO2 4D-Var, which is designed to run a longer simulation than the typical 6 h run intended for WRFDA. GH15/17 implemented a second-order checkpoint mechanism to overcome the memory limit. This approach breaks the whole simulation period into sections and saves restart files at the end of each section by the forward model. This approach requires extra calls of the forward model to recalculate the trajectory for each section during backward integration (see Fig. 3 of GH15). In WRF-CO2 4D-Var, a different approach was implemented to overcome this memory limit: the forward model saves the trajectory at each time step in memory, as WRFPLUS does. After a number of integration steps, the memory on each task processor core is dumped to an external file, and the memory is then reused. Each external file is marked with its starting time stamp and the processor core it belongs to. For instance, a 24 h simulation with 120 s time step will have a total of 720 steps. If the system saves its trajectory to external files each 30 time steps, memory allocation on each task processor core is only needed for 30 steps instead of 720 steps. This will result in 24 (720/30) trajectory files on each task processor core, and the total number of trajectory files depends on the number of processor cores used. These trajectory files are read by both tangent linear and adjoint models in a similar way to standard WRF auxiliary files. In the above example, they are read in at each 30 time steps, substantially reducing I/O time compared with reading in at each step. These trajectory files are different from standard WRF auxiliary files in that each file belongs to an individual processor core, rather than being shared among all processor cores. This means all model runs in an inverse experiment must use the same domain patch configuration, which is the most common practice.

In future development, we plan to implement observation operators for real observations, including those from towers, satellites, and airborne instruments. This is required for applying WRF-CO2 4D-Var with real observations. As a regional inverse system, the correct treatment of tracer lateral boundary conditions is important. We plan to test the lateral boundary condition adjoint code (Sect. 4) in a follow-up study. In addition, future applications of WRF-CO2 4D-Var with real observations must use proper treatment of observation and background error covariance, which was not tackled in the pseudo-observation tests in the present paper.

In addition, we also plan to periodically update the WRF-CO2 4D-Var system to keep up with WRF system updates. Such updates will mainly consist of replacing the forward model with the updated WRF code and developing the tangent linear and adjoint code for the relevant updated procedures. As the variable dependence analysis (Sect. 2.4.1) indicates that the tangent linear and adjoint code is only needed for a portion of WRF procedures, the amount of work required for updating WRF-CO2 4D-Var is manageable. In addition, future development of WRF-CO2 4D-Var will also be dependent on updates to WRFPLUS, which has always been updated along with WRF.

WRF-CO2 4D-Var source code can be retrieved via

The supplement related to this article is available online at:

The authors declare that they have no conflict of interest.

The authors express their appreciation for the WRF/WRF-Chem/WRFDA/WRFPLUS development teams for making their code available in the public domain. Discussion with Joel LeBlanc of Michigan Tech Research Institute (MTRI) improved the optimization schemes' implementation and presentation in this paper. The insightful and detailed comments from the three reviewers greatly improved both the model and the paper. This work was partially supported by a Central Michigan University CST research incentive fund. This is contribution 98 of the Central Michigan University Institute for Great Lakes Research. Edited by: Tomomichi Kato Reviewed by: three anonymous referees