2.1.1 The box model BOXMOX

BOXMOX performs box model simulations with different chemical mechanisms using varying sets of input parameters. BOXMOX relies on KPP, a code generator that simplifies the numerical integration of systems of ODEs. The temporal evolution of concentrations of chemical compounds due to photochemistry is a prime example of such a system. KPP takes a predefined set of chemical equations (in our case a chemical mechanism) written in a symbolic, human-readable language and generates a computer code (FORTRAN, Matlab, C) containing a numerical solver to integrate the system over time. A number of integration methods are available (e.g., Rosenbrock or Runge–Kutta methods). In addition to predicting the evolution of concentrations over time, the resulting solver also delivers the Jacobian and Hessian matrices of the system. Adjoint models can be generated and tuned (Sandu et al., 2003) and the Jacobians of the adjoint of the model can be obtained (see Sect. 2.2.3). Building upon KPP, BOXMOX provides additional processes typically used in box model studies (emissions, photolysis, deposition, mixing) and allows for convenient data input. BOXMOX makes simulations of chamber experiments, Lagrangian-type air parcel studies, and a description of the chemistry in the atmospheric boundary layer feasible without effort. Input is done via simple text files: initial conditions, photolysis rates, temperature, boundary layer height, detrainment and entrainment, turbulent mixing, emission, and deposition are possible input parameters. BOXMOX is a stand-alone C and Fortran program running on Linux or Mac OS X. In this work we have extended BOXMOX with an interface written in the Python language (boxmox package) to interface more easily with BEATBOX.

BEATBOX uses BOXMOX to rapidly generate a large number of simulations with perturbed input parameters. The temperature, (time-varying) photolysis rates, and initial concentrations of each species included in the investigated chemical mechanism can be perturbed independently. Ensemble members can be generated by producing normal- or lognormal-distributed perturbation factors with the possibility to adjust the mean and the standard deviation for each perturbed variable. These perturbation factors are then multiplied by the relevant initial values to produce an ensemble of initial conditions.

In this work, we demonstrate the BEATBOX capabilities using MCMv3.3.1 as NR. Because of its near-explicit representation of atmospheric chemistry, MCMv3.3.1 is a good choice to be seen as the assumed “truth” within the context of an OSSE. The MOZART-T1 chemical mechanism is employed as the simplified and/or degraded model (CR, AR).

2.1.2 Input data generation

BOXMOX comes with a tool to generate input data from field campaign observations (the genbox Python package). Translation to mechanism-specific species naming and lumping is achieved using a translation tool (the chemspectranslator package) originally based on the emission database created by Bill Carter (UC Riverside; http://www.cert.ucr.edu/~carter/emitdb). The current system uses data collected in the FRAPPE (Front Range Air Pollution and Photochemistry Experiment) field campaign (frappedata package).

During FRAPPE, a number of flight measurements with remote sensing and in situ devices of numerous quantities were performed, including concentrations of chemical species, photolysis rates, and temperature. In the examples shown we use measurements taken during FRAPPE with the NCAR C130 research aircraft to initialize the box model. Photolysis rates measured during FRAPPE are used in the examples shown here. For other cases in which photolysis rates are missing, we provide the ability to use photolysis rates calculated by the Tropospheric Ultraviolet Visible (TUV) radiation model version 5.1 (Madronich and Flocke, 1997) in BOXMOX using the tuv Python package.

In this current version of the code only measurements taken during the FRAPPE campaign have been used, but data from other field campaigns can be easily adapted as well with minimal code development.

2.2.1 The data assimilation problem

Data assimilation combines observations and model information (also called forecast or background) to derive an optimal state (analysis) with a reduced error to provide the best initial condition for a subsequent forecast (see, e.g., Lahoz et al., 2010). Consider n∈ℕ and p∈ℕ the dimensions of the state (or model) space and the observation space, respectively. Following Nichols (2010) the solution to the data assimilation problem is commonly expressed as follows:

where x_a is the analysis state, x_b the background state, y_o the observation, H the observation operator (also known as forward operator in the variational formalism), and K the gain matrix. The gain matrix K handles the transformation from the observation space to the state (or model) space. Conversely H handles the transformation from model space to observation space. The above equation can also be expressed in the incremental form, such as

with Δx called the increment and Δy called the innovation (or departure). The innovation in the observation space is then translated to the model space by applying the gain matrix K. Different methods exist to estimate the K gain; see Sects. 2.2.3–2.2.5. Commonly the gain matrix is given by

where B and R are the error covariance matrices of the background (or model) and observation, respectively. In the BEATBOX framework a single observation in a given assimilation window (p=1) and only the dimension along the chemical variables is considered. Hence, R simplifies to a 1×1 matrix, a scalar ${\mathit{\sigma }}_{\text{o}}^{\mathrm{2}}$ the observation error variance, and HBH^T also simplifies to a scalar ${\mathit{\sigma }}_{\text{b}}^{\mathrm{2}}$ (see Sect. 2.2.2 below) the background error variance in the observation space. Then BH^T can be seen as ${\mathit{\sigma }}_{\text{b}}^{\mathrm{2}}\mathit{s}$, where s can be called the sensitivity vector. The gain matrix in Eq. (3) then becomes a vector κ and can be reformulated as follows:

2.2.2 Observation generator – synthetic observations

Generating observations consists of the following steps:

• sampling values from the nature run,

• perturbing those values to simulate an observation inaccuracy, and

• specifying an observation error value.

BEATBOX forecasts have no dimensions in the 3-D atmospheric space (latitude, longitude, altitude). Sampling the NR to simulate observations is straightforward. Conventionally, H is defined as the observation operator and handles the transformation of information defined in the model space to the observation space to compare model and observation quantities. Then, the H operator can be expressed as $H^{\text{T}}=[\mathrm{0},\mathrm{0},\mathrm{\dots },\mathrm{1},\mathrm{\dots },\mathrm{0}]$. If no observation error is simulated the observation is defined as a perfect observation and the observation value is y_o=H(x_NR).

If the observation error needs to be simulated, then y_o is generated by adding a perturbation. In that version of BEATBOX the perturbation is assumed to be a normal distribution. But other probability density functions can be implemented easily (a lognormal perturbation is also implemented in this version). Then,

with N(M,Σ) as a normal distribution of mean M and standard deviation Σ. The latter two quantities can be viewed as the observation bias or accuracy (M) and precision (Σ). Finally, an associated observation error σ_o is associated with y_o for data assimilation; σ_o can account for bias and/or precision, overestimating, or underestimating those parameters depending on whether the effect of observation error needs to be tested in BEATBOX. In the following case study (see Sect. 3) we assume non-biased observation and non-biased observation error with a Gaussian distribution, leading to σ_o=Σ.

2.2.3 Adjoint sensitivity

Adjoint sensitivities can be calculated using the KPP Jacobian matrix of output from the adjoint model at a given time step. In our case, we make the approximation that the change of the state x (that includes the observed variable y) at the time step t relative to the change of the observed variable y at a previous time step t−1 (i.e., $\text{d}{\mathit{x}}_t/\text{d}{y}_{t-\mathrm{1}}$ or $\text{d}{y}_t/\text{d}{y}_{t-\mathrm{1}})$ is linear. In the current study, we assume t−1 and t at the beginning and the end of the assimilation window. The adjoint sensitivity vector of the state x to a given observed variable y at time t can be computed using the Jacobian vectors as follows:

The adjoint assimilation method runs a single forecast with the forward model and also runs the adjoint model and computes the analysis by combining Eqs. (1), (4), and (6). The innovation in observation space Δy is calculated and the state x is inferred using the κ gain that includes the adjoint sensitivity s. In this method, σ_b should be determined with ad hoc information and will not change during the cycling process. However, different methods exist to scale σ_b appropriately between each cycle. In the BEATBOX context this can be further explored using the provided OSSE framework.

2.2.4 Ensemble sensitivity

Ensemble methods run a perturbed set (ensemble) of model realizations in parallel and derive model error and sensitivity using the created ensemble. This gives multiple realizations i of the model in the observation space y_i and model space x_i. The standard deviation of the ensemble (or the ensemble spread) is used to estimate σ_b in the observation space, such as ${\mathit{\sigma }}_y=\sqrt{E\left[{\left(y_i-E\left[y_i\right]\right)}^{\mathrm{2}}\right]}={\mathit{\sigma }}_{\text{b}}$ with E as the averaging operator. Similarly, the ensemble spread in model space can be estimated as ${\mathit{\sigma }}_x=\sqrt{E\left[{\left({\mathit{x}}_i-E\left[{\mathit{x}}_i\right]\right)}^{\mathrm{2}}\right]}$. Then statistical assumptions are made to derive the sensitivity between y and x. One of the commonly used methods as described by Anderson (2003) is to draw a linear fit in the least-squares sense between y and x using the ensemble members, such as

where r_xy and σ_xy are respectively the correlation coefficient and covariance vectors between y and each chemical variable of x and ∘ is the Schur product operator. Then, the innovation in the observation space Δy_i is calculated for each ensemble member. The ensemble state vector x_i is inferred using the same k gain that includes the same ensemble sensitivity s for each ensemble member.

2.2.5 Hybrid sensitivity and further approaches

In the current version of BEATBOX the hybrid sensitivity is defined as a combination between the two approaches mentioned above: it uses an adjoint sensitivity for each ensemble member, such as

In that sense, each ensemble member has its own adjoint sensitivity calculation s_i and its own innovation in observation space Δy_i. This results in an independent or different inference on the state vector for each ensemble member x_i. This method is just an example of what can be explored with the BEATBOX system. Because of the simplicity of the code and the low dimensionality of the problem, advanced techniques can be easily implemented and analyzed. Other hybrid methods and filters, such as a polynomial filter or particle filters and their benefits for highly nonlinear systems, can be explored with ease with the proposed framework as well.

2.2.6 Inflation

For ensemble methods, to avoid the filter divergence problem (Fitzgerald, 1971), inflation algorithms are needed. In the BEATBOX framework we included the inflation method as proposed by Anderson (2007), such as

with λ called the covariance inflation factor that determines how much the ensemble members x_i are spread out from the mean E[x]. Many different methods exist to estimate λ; it can be chosen as constant over time or adaptive given the ensemble spread in observation space σ_b, observation error σ_o, and the innovation norm from the ensemble mean $\mathit{\theta }=\left|E\right[y_i]-y_{\text{o}}|$. In the version of BEATBOX presented here we calculate λ as

If λ is found to be smaller than 1, the ensemble spread is assumed to be large enough and no inflation is calculated. Straightforward computation of λ as above assumes linearity between observation space and model space, which is true in the current BEATBOX framework. More advanced ways to compute λ as described in Anderson (2007, Appendix A) have been tested and implemented in BEATBOX and can be used as well. Finally, to conserve the positive definite nature of the ensembles and also prevent forcing ensemble members to zero that would otherwise be inflated to negative values, we reduce the inflation factor iteratively on every value of the state:

This is one of several methods to conserve the physical aspects of the ensemble: it might under-disperse the ensemble in some cases of low concentrations but ensures that the inflation is kept to physical values. In the BEATBOX framework, users can easily implement and explore new inflation methods for nonlinear, definite positive, and perturbation of highly sensitive systems, as in the case of reactive gas-phase chemistry.

2.2.7 Localization

One should consider to what extent the sensitivities should be relied on. Localization algorithms try to limit the impact on the analysis of errors in the sample covariance between observations and model state variables (Mitchell and Houtekamer, 2000). Depending on the ensemble size or the mathematical assumptions to compute a sensitivity, a localization function C should be used to define where to apply the sensitivity s, such as

In the current BEATBOX framework, it is possible to specify which species should be inferred, e.g., ${\mathbf{C}}^{\text{T}}=[\mathrm{0},\mathrm{1},\mathrm{\dots },\mathrm{1},\mathrm{\dots },\mathrm{0}]$.

3.2.1 NO₂ assimilation

We assimilate NO₂ observations using two different localizations. The first experiment infers only NO₂ concentration and the second infers the entire state vector (i.e., NO₂, O₃, CH₂O, and OH concentrations). After looking at the evolution of NO₂ concentrations (see Fig. 5), all three assimilation methods in both experiments (adjoint, ensemble, and hybrid) tend to move the analysis closer toward the observations and hence the NR. Some differences among the three ARs can be found. Recalling Sects. 2.2.1 and 2.2.2, because of the dimensionality of the problem the sensitivity from observation space NO₂ to the model variable NO₂ is equal to identity, s=1, for every assimilation method. Hence, the difference among assimilation runs will only result from differences in x_b and σ_b after a subsequent forecast. The ensemble and hybrid methods show a quicker improvement than the single member adjoint method due to the adaptive nature of σ_b with the inflation (see Sect. 2.2.6). The single adjoint method keeps σ_b=2.5 %, while the ensemble and hybrid methods can tune this value with the inflation.

When only NO₂ is inferred the other species are modified and this could be called the model response to assimilated changes: NO₂ is decreased, which increases the OH availability for other oxidation pathways. A slight increase in O₃ is noted and more significantly for CH₂O during the transition region. The ensemble methods have a stronger effect due to the adaptive nature of σ_b. The adjustment of NO₂ concentration towards NO₂ observation can be more effective and this can have a stronger effect on the model response.

When the entire state vector is inferred, i.e., no localization, slight additional increases occur in O₃, CH₂O, and OH in the transition region. In general, in the first 10 to 15 h when VOC-limited conditions dominate (see Fig. 3), the impact of NO₂ assimilation is low. By definition, VOC-limited air is insensitive to changes in NO₂, so if the model (CR) predicts VOC-limited conditions either adjoint sensitivities or ensemble sensitivity will remain small. After 15 h of simulation, the chemical regime transitions significantly to NO_x limited with higher sensitivity of the state to changes in NO₂. After 40 h, NO₂ concentrations drop to very low values, NO₂ assimilation increments are very small, and hence almost no inference on the other state variables is observed.

In the case study, inference from NO₂ observation on the rest of vector is not the major reason for improvement. Model response from NO₂ changes is mostly responsible for improving the state. NO_x concentrations drive the chemistry in the transition region and NO_x-limited regimes. NO_x chemistry is well known and similarly represented between the NR and the CR. Hence the model response is likely to improve the state and not likely to produce additional errors.

Figure 5Evolution of the state vector over 60 h, including the nature run (NR, black line), control run (CR, blue), and assimilation runs using adjoint (red), ensemble (pink), and hybrid (turquoise) methods. Two experiments with different localizations are displayed: only NO₂ inferred (a–d) and whole state vector inferred (e–h). Shaded areas show corresponding ensemble spread for ensemble and hybrid methods.

Download

3.2.2 CH₂O assimilation

We repeat the experiments presented in Sect. 3.2.1, but assimilate CH₂O observations instead of NO₂. Figure 6a–d show the concentration evolution when only CH₂O is inferred. In the first 40 h CH₂O is underestimated by the CR mechanism. All three ARs show very similar results. During the analysis phases the CH₂O concentrations are pushed up towards the NR. Those abrupt increases are systematically compensated for by the chemical mechanism nudging the system back into chemical equilibrium, i.e., towards the CR concentrations. This makes the CH₂O concentration evolution a sawtooth shape. CH₂O is a short-lived chemical compound (shorter than NO₂ in this case study) and it is mainly driven by other chemical species concentrations. Ultimately CH₂O concentrations in the ARs are slightly elevated after each 3 h forecast. These changes will very mildly affect the state vector concentrations. The model response to CH₂O changes in NO₂, O₃, and OH is slight and definitely smaller than the model response to NO₂ changes (see Sect. 3.2.1).

When the entire state vector is inferred, i.e., no localization (Fig. 6e–h), we observe very different results. The two ARs that are using adjoint sensitivities (hybrid and adjoint) show similar results as the previous experiment, and the sawtooth shape is still observed. The analysis shows strong increases in CH₂O and the chemical mechanism tries to come back to the CR state. The inference on the other species of the state vector is noticeable, and NO₂, O₃, and OH concentrations are improved compared to the experiment when only NO₂ is inferred.

In the ensemble method in which no localization is applied, the improvement in CH₂O is of a different nature. The sawtooth behavior of the AR is not observed anymore and the chemical mechanism now seems to have changed from systematic CH₂O loss to CH₂O production in the forecast. This will drive the AR to better fit the observations and hence the NR. The inference on other species shows significant changes that will in general change the sign of the error and sometimes increase it; NO₂ becomes underestimated in the transition region, O₃ is overestimated but then underestimated after 35 h, and OH is overestimated in the transition region. In the ensemble method for the case study, the computed sensitivities seem significantly different from the adjoint sensitivities. This will allow us to change the chemical production and loss rate of CH₂O to better adjust the CH₂O concentration, but at the risk of disturbing the rest of the state vector and increasing errors that might lead to unphysical results. We diagnose the difference among sensitivities from the two experiments presented above in Sect. 3.3. We also diagnose the chemical behavior change from CH₂O assimilation using fluxes in Sect. 3.4.

Figure 6Same as Fig. 5 but with CH₂O observations assimilated.

Download