BEATBOX v1.0: Background Error Analysis Testbed with Box Models

The Background Error Analysis Testbed (BEATBOX) is a new data assimilation framework for box models. Based on the BOX Model eXtension (BOXMOX) to the Kinetic Pre-Processor (KPP), this framework allows users to 10 conduct performance evaluations of data assimilation experiments, sensitivity analyses and detailed chemical scheme diagnostics from an Observation Simulation System Experiment (OSSE) point of view. The BEATBOX framework incorporates an observation simulator and a data assimilation system with the possibility of choosing ensemble, adjoint or combined sensitivities. A user-friendly, python-based interface allows tuning of many parameters for atmospheric chemistry and data assimilation research as well as for educational purposes, e.g. observations error, model covariances, ensemble size, 15 perturbation distribution on initial conditions, and so on. In this work, the testbed is described and two case studies are presented to illustrate: the design of a typical OSSE experiment, data assimilation experiments, a sensitivity analysis and a method for diagnosing model errors. BEATBOX is released as an open source tool for the atmospheric chemistry and data assimilation communities.

1 Introduction 20 Current regional and global models of the composition of Earth's atmosphere exhibit a high level of complexity due to the combination of chemical and meteorological processes such as transport, thermodynamics, radiation or precipitation. But 'just' gas phase chemistry itself already has a considerable amount of complexity. The number of variables (chemical compounds) and equations (chemical reactions) in current model representations of tropospheric chemistry ('mechanisms') can vary by two orders of magnitude (from 10 2 to 10 4 ). Compare, for example, two well-known mechanisms: The Model for 25 Ozone And Related Chemical Tracers, version T1 (MOZART-T1; Emmons et al., 2010;Knote et al., 2014) uses 134 species and 250 reactions, whereas the Master Chemical Mechanism, version 3.3 (MCMv3.3; Jenkin et al., 2015) employs over 5000 species and over 15 000 reactions. MCMv3.3 is called a 'near-explicit' chemical mechanism, which describes in detail the gas-phase chemical processes involved in the tropospheric degradation of volatile organic compounds (VOC). On the other hand, MOZART-T1 uses a simplified ('lumped') representation of VOCs. This reduction in complexity is advantageous as it 30 reduces the computational demand, but can lead to significant errors in the prediction of atmospheric composition. The choice of a chemical mechanism is therefore a trade-off: while MCMv3.3 would be desirable due to its fidelity in representing atmospheric chemistry, it cannot currently be used in large scale 3-D atmospheric simulation due to its computational demand. MOZART-T1 is less accurate, but also much more economic.
Investigating the performance of chemical mechanisms is often done using zero-dimensional (0-D) box models. Numerous 5 studies have used box models to study reactive gas phase chemistry and provide intercomparisons and validation of mechanisms. Emmerson et al. (2009) provided a comprehensive comparison of the MCM mechanism with six tropospheric chemistry lumped schemes that could be used within chemistry transport models. Archibald et al. (2010) performed an intercomparison of the gas phase mechanism for isoprene degradation in a box model with various mechanisms widely used in 3-D models. More recently, Coates et al. (2015) compared MCM to simplified VOC mechanisms to look at ozone 10 production for a selection of VOCs representative of urban air masses. Knote et al. (2015) conducted box model simulations using different chemical mechanisms and compared them to each other to understand mechanism-specific biases during a 3-D model intercomparison. Mazzuca et al. (2016) used an observation-constrained box model with a lumped carbon-bond mechanism to study photochemical oxidation and ozone production processes along a research aircraft campaign. Wolfe et al. (2016) presented a tool for 0-D atmospheric modeling that can use different chemical mechanisms and methods of 15 photolysis frequencies calculations.
This non-exhaustive list of recent studies using box models to study tropospheric chemistry shows the importance and need of such tools. Box models are attractive due to their simplicity and low computational cost, but cannot provide a realistic representation of the entire atmosphere because they lack vertical and horizontal diffusion, boundary conditions and numerous other processes that 3-D models take into account. In that regard, box-model studies should not focus on 20 replicating the most accurate predictions but rather aim at gaining significant fundamental and conceptual understanding of a given system of ordinary differential equations (ODE), in the present case the one that governs tropospheric chemistry.
Another aspect of box-models or reduced complexity models is the applicability to data assimilation research. Low dimension and/or box-models are often used to design new data assimilation algorithms and to prove conceptually the advantage of a given method. In data assimilation theory, use of the Lorenz model and other simple systems is common 25 practice: Sandu et al. (2005) used a box model approach to design an adjoint-based sensitivity analysis for reactive gas phase tropospheric chemistry. Ott et al. (2004) used a Lorenz-96 model to introduce new formulation of the ensemble Kalman filter approach. Van Leeuwen (2010) also used the Lorenz model to investigate non-linear advanced data assimilation technique such as the particle filter.
Atmospheric composition data assimilation, and more generally inversion, is complex and is computationally demanding. 30 Reactive gas-phase photochemistry is highly non-linear and has to deal with hundreds to thousands of variables. There is a need for a tool that allows exploring suitable novel data assimilation approaches for atmospheric chemistry, assess uncertainty and errors of a given chemical mechanism and perform chemical sensitivity analysis from a parameter to another. All this should be possible at minimal computational expenses and coding skills requirement. In this paper, we present a new framework based on BOXMOX (BOX MOdel eXtension to KPP; Knote et al., 2015) called BEATBOX (Background Error Analysis Testbed with Box Models) that is able to cycle data assimilation windows, calculate adjoint, ensemble or even more advanced sensitivity analysis and assess box-model errors and uncertainties. In section 2 we present in detail the structure of BEATBOX and its algorithms, exemplified through case studies which we discuss in section 3.

Design of BEATBOX 5
The Background Error Analysis Testbed with Box Models (BEATBOX) is a suite of tools that allows for a simple and fast investigation, comparison and evaluation of any system of Ordinary Differential Equations (ODE) evolving over time. By "Background Error" we designate a general term for model error characterization. In this study, BEATBOX is used within the scope of atmospheric chemical mechanisms. Currently, the BEATBOX framework consists of a forecasting tool, the BOX MOdel eXtension (BOXMOX; Knote et al., 2015) to the Kinetic PreProcessor (KPP; Sandu and Sander, 2006), 10 presented in Sect. 2.1, and a data assimilation tool, which will be introduced in Sect. 2.2.
The BEATBOX structure is built upon Observing System Simulation Experiments (OSSE; Arnold and Dey, 1986). OSSEs are generally used in the field of numerical weather prediction (e.g. Kuo et al., 1998;Wang et al., 2008;Liu et al., 2009) as well as atmospheric composition and air quality predictions (e.g. Edwards et al., 2009;Claeyman et al., 2011;Barré et al., 2015Barré et al., , 2016. OSSEs allow assessing the benefit of a potential new type of instrument for environmental predictions using a 15 data assimilation system and are of crucial importance to define requirements of a given instrument. Space agencies such as the National Aeronautics and Space Agency (NASA) or European Space Agency (ESA) hence support OSSEs as tools to proof scientific readiness levels for proposed space missions. Also, the model and data assimilation requirements should be assessed to meet a required predictive capability. The current BEATBOX framework employs a box modeling approach to avoid the space dimension problem -the dimensionality of the problem is reduced to time and variables only. Hence the 20 geometry, radiative specifications, and spatial resolution of a new instrument type are irrelevant. What can be easily explored within BEATBOX is the type of data assimilation method used (e.g. background error covariance calculations, localizations, inflation methods), the revisit time (or temporal sampling) of a measurement, the variable(s) observed, and straightforward comparisons between sets of ODEs (chemical schemes in our case).
A number of issues on the OSSE technique should be mentioned as well. Performing an OSSE could be costly in terms of 25 setup and design as well as computationally. Numerical integration of the most state of the art representation of the earth system for sampling observations and benchmarking with could be intensively costly and requires highly skilled staff and extensive collaboration between research entities. Approximations are often required to make experiments possible (e.g. the "identical twin" problem ) necessitating careful diagnosis of the results and could limit scientific conclusions.
Ultimately an OSSE should be used to highlight model deficiencies and inaccuracies, and provide direct guidance for model 30 improvement. In that context, BEATBOX could be considered as a derived OSSE framework focused on data assimilation technique and model improvement rather than the benefit of new or future type of observations. Starting from a scientific question or hypothesis to be validated (or rejected) that fits the topics mentioned above, several components are required ( Figure 1): • A nature run (NR) considered as the "true state". The NR supposes to use the best model representation possible considering the state of the art. In this study, we used the Master Chemical Mechanism (MCMv3.3.1) as the NR.
• A control run (CR) as the prior estimated state of the atmosphere. Compared to the NR, a simplified or degraded 5 model should be used for example a set of ODEs that can be implemented in large scale 3-D models. In this study, we use MOZART-T1.
• An observation simulator that generates synthetic observation by sampling the NR. Observation errors also need to be simulated.
• An assimilation run (AR) that is produced using the data assimilation tool merging the synthetic observation with 10 the CR to produce the best estimate possible of the state.
• A suite of diagnostic tools that use NR, CR and AR designed to be able point out model and data assimilation technique limitations, ultimately providing a direct feedback for model improvement. The BEATBOX framework has the capability to loop over several assimilation cycles, also called 'cycling'. Cycling with BEATBOX is schematically displayed in Fig. 2. Every cycle starts with the forecast step. By applying their respective model, the NR, the forecast (F) and the CR are processed from cycle − 1 to cycle . Then, is used to generate synthetic observations through the observation operator (see Sect. 2.2.2). These observations are assimilated into the forecast to produce the analysis with, the data assimilation method of choice using a gain (see . 20 Then, will be used to generate new initial conditions to start a new forecast. serves as a reference to determine the performance of the assimilation method after forecast cycles. The forecast can be taken to quantify the skill of the assimilation method at the current cycle . Afterwards, the cycle + 1 starts with its forecast step. BOXMOX relies on KPP, a code generator which 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 10 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 the 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 15 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, de-/entrainment, turbulent mixing, emission, and deposition are possible input parameters. BOXMOX is a standalone C / 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. 20 BEATBOX uses BOXMOX to rapidly generate a large number of simulations with perturbed input parameters. 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 log-normaldistributed perturbation factors with 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. 5 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/degraded model (CR, AR).

Input data generation
BOXMOX comes with a tool to generate input data out of field campaign observations (the genbox Python package). 10 Translation to mechanism-specific species naming and lumping is achieved using a translation tool (the  (Madronich and Flocke, 1997) in BOXMOX using the tuv Python package. 20 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.

Data assimilation tool
The beatboxtestbed Python package provides a data assimilation tool that sample observations from the NR, with possibility of tuning the observation error parameters. With assimilating observations, different sensitivity analyses can be 25 used, such as adjoint, ensemble or combined sensitivity (in this paper called hybrid) are included in this version of BEATBOX (see Sect. 2.2.3-5). For notation purposes, we represent x and y as variables in model or state space and observation space, respectively.

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 ∈ ℕ and ∈ ℕ 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: 5 where is the analysis state, the background state, the observation, the observation operator (also known as forward operator in the variational formalism) and the gain matrix. The gain matrix handles the transformation from the observation space to the state (or model) space. Conversely handles the transformation from model space to observation space. The above equation can also be expressed in the incremental form such as: 10 with ∆ called the increment and ∆ called the innovation (or departure). The innovation in the observation space is then translated to the model space by applying the gain matrix . Different methods exist to estimate the gain, see Sect. 2.2.3-

5.
Commonly the gain matrix is given by: where and 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 ( = 1) and only the dimension along the chemical variables is considered. Hence, simplifies to a 1 × 1 matrix, a scalar 2 the observation error variance, also simplifies to a scalar 2 (see Sect. 2.2.2 below) the background error variance in the observation space. Then can be seen as 2 where can be called here the sensitivity vector. The gain matrix in Eq. (3) then becomes a vector and can 20 be then reformulated as follows:

Observation generator -synthetic observations
Generating observations consists of the following steps: • Sampling values from the nature run 25 • Perturbing those values to simulate an observation inaccuracy 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 log-normal perturbation is also implemented in this version). Then, with (Μ, Σ) a normal distribution of mean Μ and standard deviation Σ. The latter two quantities can be viewed as the 5 observation bias or accuracy (Μ) and precision (Σ). Finally, an associated observation error is associated with for data assimilation. can account for bias and/or precision, overestimating or underestimating those parameters, depending if 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 = Σ.

Adjoint sensitivity 10
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 (that includes the observed variable ) at the time step relative to the change of the observed variable at a previous time step − 1 (i.e. / −1 or / −1 ) is linear. In that current study, we assume − 1 and at the beginning and the end of the assimilation window. The adjoint sensitivity vector of the state to a given observed variable at time can be computed using the Jacobian vectors as follows: 15 The adjoint assimilation method runs a single forecast with the forward model and also runs the adjoint model and computes the analysis combining Eq. (1), (4) and (6). The innovation in observation space ∆ is calculated and the state is inferred using the gain that includes the adjoint sensitivity . In this method, should be determined with ad-hoc information and will not change during the cycling process. However different methods exist to scale appropriately between each cycle. In 20 the BEATBOX context this can be further explored using the provided OSSE framework. . Then statistical assumptions are made to derive the sensitivity between and . One of the commonly used methods as described by Anderson (2003) is to draw a linear fit in the least-squares sense between and using the ensemble members, such as

Ensemble sensitivity
where and are respectively the correlation coefficient and covariance vectors between and each chemical variable of and ∘ the Schur product operator. Then, the innovation in the observation space ∆ is calculated for each ensemble member. The ensemble state vector is inferred using the same gain that includes the same ensemble sensitivity for each ensemble member.

Hybrid sensitivity and further approaches 5
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 on each ensemble member, such as In that sense, each ensemble member has its own adjoint sensitivity calculation and its own innovation in observation space ∆ . This results to independent or different inference on the state vector for each ensemble member . This method is 10 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 non-linear systems, can be explored with ease with the proposed framework, as well.

Inflation 15
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 by how much the ensemble members are spread out from the If is found smaller than one, 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 25 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 such as: 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 assures the inflation is kept to physical values. In the BEATBOX framework, users can easily implement and explore new inflation methods for non-linear, definite positive and perturbation of highly-sensitive systems, as in the case for reactive gas-phase chemistry. 5

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 should be used to define where to apply the sensitivity , such as: 10 In the current BEATBOX framework, it is possible to specify which species should be inferred, e.g. = [0, 1, … ,1, … ,0].

Flux method for model analysis
A flux-tool has been included in boxmox that calculates the production and loss fluxes for a given chemical component.
Consider the generalized chemical reactions such as, 15 with called the rate constant. The chemical flux can be expressed as For a given chemical component, a detailed budget of chemical production and loss can be made identifying different chemical reactions. An example of the application of the flux tool is shown in Sect. 3.4. 20

Application to an urban pollution case study
In this section, we provide an example case study to showcase the capabilities of the BEATBOX framework.

Control Run and Nature Run
Temperature, concentrations and photolysis rates from the FRAPPE data (see Sect. 2.1.2) are used for initial conditions. In the present simulations, all environmental parameters such as temperature and photolysis rates are kept constant. Interactions 25 between the simulated box and the surrounding are neglected, the box model simulation can be seen as a 'chemistry in a jar', similar to a chamber experiment without consideration of wall losses, where only the temporal evolutions due to chemical reactions is allowed. Initial conditions for primary VOCs and inorganic compounds are provided using the FRAPPE observations.
The present case study focuses on an air mass originating from industrial area Commerce City, near Denver, Colorado. Figure 3 shows the estimated temporal evolution over the first 60 hours of simulation for the VOC/NO x (NO x = NO + NO 2 ) 5 and toluene/benzene-ratio using the MOZART-T1 scheme. The vertical lines suggest the VOC-limited (VOC/NO x ratio <= 4), NO x -limited (VOC/NO x ratio > 15) and the transition region (4 < VOC/NO x ratio <= 15), to show that the simulation transitions through different O 3 production regimes with possibly very different relevant chemical pathways. The measurements show an initially strongly VOC-limited air indicative of an urban air mass. The VOC/NO x -ratio of the aging air increases in time. During the first 15 h the simulation shows a strong VOC-limited regime. A transition regime spans 10 from 15 h to approximately 30 h. After the transition period the chemical regime becomes NO x -limited, representative of more rural / background conditions. The toluene/benzene-ratio is used as qualitative measure of photochemical age. Toluene and benzene are considered to have the same sources but toluene is more quickly oxidized than benzene which leads to a decline in the ratio over time.

Figure 3: Temporal evolution of the VOC/NO x (blue line) and toluene/benzene (red line) ratios over 60h, derived from MOZART-T1 (CR).
The temporal evolution of some key species in the NR and the CR are displayed in Fig. 4. Nitrogen dioxide (NO 2 ) shows a decrease over time, with a stronger decay in the NR than in the CR. Ozone (O 3 ) production is observed over time with hence a stronger production in the NR than in the CR. The hydroxyl-radical (OH) availability is increasing over time especially 20 between 20 and 35 hours which can be identified as the end of the transition from a VOC-limited to a NO x -limited regime. a decrease followed by an increase and finally a decrease. Those variations illustrate the change of chemical regimes from VOC-limited to transition to NO x -limited. If we compare the NR to the CR the change of chemical regimes is lagged. The change of regimes seems to occur faster in the NR than in the CR suggesting a different reactivity and pathways of oxidation in the NR than in the CR. 5

Assimilation runs
We focus on assimilating NO 2 and CH 2 O concentrations in two separate experiments to show the capabilities of the 10 BEATBOX framework. This is motivated by the fact that NO 2 and CH 2 O play a key role in the atmospheric chemistry especially for short-lived chemical compounds. Also, NO 2 and CH 2 O are, and will be, observed from space from both low earth and geostationary orbiters and from in-situ observations. We define a reduced state vector for simplicity of the demonstration with the species described above: NO 2 , O 3 , CH 2 O, OH.
We use an assimilation window of 3 hours. Unbiased observations have been defined and assumed to be at the end of the 15 assimilation window with an observation error of 0.75ppbv for NO 2 and 0.07ppbv for CH 2 O (corresponding to approximately 2.5% and 5% of the concentrations at initial time). All the generated observations have well estimated observation error such as, = Σ (see Sect. 2.2.2). The model error in the observation space has to be specified for the adjoint technique and is set to 2.5% for NO 2 and 5% for CH 2 O of the current forecast concentration (see Sect. 2.2.3). For the ensemble method, the ensemble spread implicitly defines using 50 ensemble members which are generated by perturbing the initial NO 2 concentration by 2.5% for NO 2 assimilation and by perturbing the initial CH 2 O concentration by 5% for 5 CH 2 O assimilation (see Sect. 2.2.4). Adaptive inflation is applied on the ensembles at every assimilation window to maintain a realistic ensemble spread to weight observation and model appropriately (see Sect. 2.2.6).

NO 2 assimilation
We assimilate NO 2 observations using two different localizations. The first experiment infers only NO 2 concentration and the second infers the entire state vector (i.e. NO 2 , O 3 , CH 2 O and OH concentrations). After looking at the evolution of NO 2 10 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 between the three AR can be found.
Recalling Sect. 2.2.1 and Sect. 2.2.2, because of the dimensionality of the problem the sensitivity from observation space NO 2 to the model variable NO 2 is equal to identity, = 1 for every assimilation method. Hence, difference between assimilation runs will only result from differences in and after a subsequent forecast. The ensemble and hybrid 15 methods show a quicker improvement than the single member adjoint method due to the adaptive nature of with the inflation (see Sect. 2.2.6). The single adjoint method keeps = 2.5% while ensemble and hybrid methods can tune this value with the inflation.
When only NO 2 is inferred the other species are modified and this could be called the model response to assimilated changes: NO 2 is decreased that is increasing the OH availability for other oxidation pathways. A slight increase of O 3 is 20 noted and more significantly for CH 2 O during the transition region. The ensemble methods have a stronger effect due to the adaptive nature of . The adjustment of NO2 concentration towards NO 2 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 on O 3 , CH 2 O and OH in the transition region. In general, in the first 10 to 15 hours, when VOC-limited conditions dominate (see Fig. 3) the impact of 25 NO 2 assimilation is low. By definition, VOC-limited air is insensitive to changes in NO 2 , so if the model (CR) predicts VOC-limited conditions either adjoint sensitivities or ensemble sensitivity will remain small. After 15 hours of simulation, the chemical regime transitions significantly to NO x -limited with higher sensitivity of the state to changes in NO 2 . After 40 hours, NO 2 concentrations drop to very low values, NO 2 assimilation increments are very small, hence almost no inference on the other state variables is observed. 30 In that case study, inference from NO 2 observation on the rest of vector is not the major reason for improvement. Model response from NO 2 changes is mostly responsible for improving the state. NO x concentrations drive the chemistry in the

CH 2 O assimilation
We repeat the experiments presented in 3.2.1, but assimilate CH 2 O observations instead of NO 2 . In Fig. 6, the left-hand side plots show the concentration evolution when only CH 2 O is inferred. In the first 40 hours CH 2 O is underestimated by the CR mechanism. All three AR show very similar results. During the analysis phases the CH 2 O concentrations are pushed up towards the NR. Those abrupt increases are systematically compensated by the chemical mechanism nudging the system 5 back into chemical equilibrium, i.e. towards the CR concentrations. This makes the CH 2 O concentration evolution a saw tooth shape. CH 2 O is a short lived chemical compound (shorter than NO 2 in this case study) and it is mainly driven by other chemical species concentrations. Ultimately CH 2 O concentration in the ARs are slightly risen up after each 3-hour forecast.
These changes will very mildly affect the state vector concentrations. The model response of CH 2 O changes on NO 2 , O 3 and OH is slight and definitely smaller than the model response on NO 2 changes (see Sect. 3.2.1). 10 When the entire state vector is inferred, i.e. no localization (Fig. 6, right-hand side plots) we observe very different results.
The two ARs that are using adjoint sensitivities (hybrid and adjoint) show similar results than the previous experiment, the saw tooth shape is still observed. The analysis shows strong increases on CH 2 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, NO 2 , O 3 and OH concentrations are improved compared to the experiment when only NO 2 is inferred. 15 In the ensemble method where no localization is applied, the improvement on CH 2 O is of a different nature. The saw-tooth behavior of the AR is not observed anymore and the chemical mechanism now seems to have changed from systematic CH 2 O loss to CH 2 O production in the forecast. This will drive the AR to fit better the observations and hence the NR. The inference on other species show significant changes that will in general change the sign of the error and sometimes increases it; NO 2 becomes underestimated in the transition region, O 3 is overestimated but then underestimated after 35 hours, and OH 20 is overestimated in the transition region. In the ensemble method for that case study, the computed sensitivities seem significantly different from the adjoint sensitivities. This will allow to change the chemical production/loss rate of CH 2 O to adjust better the CH 2 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 between sensitivities from the two experiments presented above in Sect. 3.3. We also diagnose the chemical behavior change from CH 2 O assimilation using fluxes in Sect. 3.4. 25  Figure 7 shows the sensitivities of the unobserved species to a change in NO 2 and CH 2 O respectively at the end of each 3 h assimilation window with the single member adjoint and the ensemble method (see Sect. 2.2.3 and 2.2.4). To make sensitivities comparable we have normalized (divided) them by the state concentrations. For example, ( 2 , 3 ) the sensitivity of NO 2 changes over O 3 will be most likely orders of magnitude larger than ( 2 , ) since O 3 concentration 5 are orders of magnitude larger than OH concentrations.
For NO 2 assimilation both methods deliver similar results. The sensitivities are in general small, not above 60%. Sensitivities are small or negligible in VOC-limited regime but become more significant during the transition and the NO x -limited regimes, mostly after 30 hours. The changes on NO 2 from assimilation after 30 hours are rather small and hence the inference on other species will be small. This then confirms that the most important part of the changes from NO 2 10 assimilation is due to the model response from NO 2 changes and secondly due to the data assimilation inference on the other species of the state.
For CH 2 O assimilation we saw significant differences in behavior between adjoint and ensemble methods to compute the sensitivities. The adjoint displays rather small sensitivities, peaking at 20% but mostly below 10%. Those sensitives are observed a during the VOC-limited regime when the chemistry is sensitive to changes from VOC concentrations, and 15 disappear during the transition region and become insignificant in the NO x -limited regime. This explains the small but reasonable impacts from CH 2 O changes on the rest of the state. Concerning the ensemble method, the computed sensitivities are much larger and do not decrease after the transition regime. Values switch abruptly from negative to positive and are sometimes above 200%. To understand this unphysical behavior, we display tracer-tracer relationships during the assimilation phase of the 9 th cycle (27 hours), in Fig. 8. 20 In Sect. 2.2.4, we defined the ensemble method sensitivity computation as a linear-fit to the ensemble distribution between two species. Most state-of-the-art EnKF methods use this approximation. One can see the limitation to such linear assumption after looking carefully at Fig. 8. For example, on the right-top corner are displayed the prior and posterior states of the NO 2 -O 3 distributions. The ensemble method represents the observation well; the ensemble distribution is at the NO 2 level of the observation and the ensemble spread fits the observation error (green vertical bars). However, the relationship 25 between NO 2 -O 3 that has formed during the ensemble AR after 9 cycles is strongly curved and non-linear, which is difficult to represent through the linear fit of the ensemble sensitivity In this example, moving the ensemble members along the linear-fit is moving the O 3 distribution slightly away from the NR. At the same time the adjoint sensitivities are in comparison very weak (the slope is along the y-axis), slightly improving the state distributions but not very strongly. Finally, the hybrid distribution show lager spread and different distribution shape than the ensemble distribution. The adjoint 30 inference does not strongly change the rest of the state which maintains the chemical production/loss rates of CH 2 O that fall into a temporary attractor; the model wants to go back to the CR concentrations and creates a sawtooth-shaped concentration evolution over time (see Sect. 3.2.2). After 3 hours, the forecast will be significantly far from the observation and the distribution will be inflated. The ensemble inference strongly changes the rest of the state, which will change the chemical production/loss rates and drive the system out of the attractor. To understand this more clearly, the flux tool is presented in the following section.

Model diagnostics using fluxes
In this last section, we will focus on the CH 2 O fluxes (production and loss) over the 9 th cycle forecast (25 to 27 hours) of the CH 2 O assimilation without localization. In Fig. 9 are displayed the individual contributions of production and loss of CH 2 O for the NR and one same member of the CR, AR ensemble and AR hybrid. The flux tool isolates a few different reactions with NR because the chemical scheme is different and more detailed than the CR. In this case study, CH 3 O → CH 2 O + HO 2 5 in the NR is simply an intermediate reaction of NO + CH 3 O 2 → CH 2 O + NO 2 + HO 2 . Other than this reaction, the isolated reactions are similar.
The NR fluxes to the CR fluxes show the same order of production/loss importance, however the NR fluxes are stronger.
Overall the loss terms in the NR are stronger than the CR leading to a net chemical flux that is negative. The CR rates are significantly faster than the AR ensemble rates, and the slopes of the rates (i.e. the second derivative of the concentration 10 evolution) also differ. . The production terms are stronger leading to almost no loss net flux. The order of importance of the reactions in the budget is also different; CH 2 O + OH → CO + HO 2 + H 2 O is much more important, likely due to the increase of OH during the AR ensemble (see Sect. 3.2.2). Finally, the AR hybrid fluxes are similar to the CR fluxes with the same order of importance between reactions. Except for the carbon monoxide (CO) formation from photolysis, CH 2 O → CO + H 2 that is now stronger and that is responsible for the systematic decrease of CH 2 O in the AR adjoint and AR hybrid (see Sect. 15 3.2.2). This increased flux from this reaction explains the saw tooth behavior of CH 2 O, when no other pathway of loss is possible, i.e. OH is not as strongly changed in AR adjoint and hybrid as in AR ensemble, this way of destroying CH 2 O is then increased.
We demonstrate here the usefulness of the flux tool to diagnose the effect of data assimilation on atmospheric chemistry in a detailed manner. The diagnostic of the flux tool can be used on any species that a chemical scheme contains for any kind of 20 BOXMOX simulation. In this paper, we have presented a new suite of tools for box-models, BEATBOX. The design of BEATBOX is based on an OSSE approach that can integrate simultaneously various chemical schemes to simulate a nature run, control runs and assimilation runs. This framework includes the capability of running assimilation windows of different chemical schemes 5 using a forecasting tool (BOXMOX) and an assimilation tool allowing sensitivity analysis. BEATBOX provide ensemble and adjoint sensitivity analysis that can be combined or modified to explore new inverse or data assimilation methodologies.
Additionally, a flux tool is also integrated into the framework to diagnose and assess in detail models runs differences and ultimately use data assimilation to improve the model (chemical scheme) itself and not only the model outputs. The systematic and detailed assessment of the multivariate data assimilation problem that BEATBOX can tackle important 10 scientific hypotheses at a limited computational cost for future data assimilation configuration on large scale 3D models for the atmospheric chemistry but not only limited to it. Any system of equations can be integrated over time in the current framework.
A typical case study of ozone photochemical production from NO x is presented to showcase the capability of BEATBOX.
Differences between the nature run and the control run are presented followed by a data assimilation experiment of synthetic 15 NO 2 and CH 2 O observations using adjoint and ensemble sensitivity analysis with different localization parameters. The case studies shown in this paper illustrate the need to understand in detail the effect of data assimilation in as complicated and non-linear a model that the atmospheric chemistry requires. In these case studies, we showcased that BEATBOX is a powerful and user-friendly tool for: • Understanding chemical mechanism differences and deficiencies 20 • Performing chemical sensitivity analysis using ensemble or adjoint methods • Envisioning and designing new inverse and data assimilation methods for atmospheric chemistry to optimally constrain as much of the chemical state as possible • Defining requirements for new chemical and data assimilation schemes and ultimately improving them.
• Educational purposes for data assimilation and atmospheric chemistry 25 BEATBOX will continue to evolve depending on user requirements. For example, emission inversion capability is currently being implemented. Setting any given observation time into the assimilation window will also be possible. Assimilating multiple observations in a given assimilation window will be also implemented using sequential and variational minimization techniques, as well. Using real observations, i.e. from field campaigns is also possible with minimal code modifications. Because of the user friendliness, flexibility and the open source nature of most the BOXMOX/BEATBOX, 30 users could also contribute to model development making it a broad atmospheric chemistry community tool.

Code and/or data availability
All code developments presented here are open source tools released under the GNU General Public License v3. BEATBOX consists of a number of Python packages: • beatboxtestbed: the BEATBOX Background Error Analysis Testbed 5 • boxmox: Python interface for the chemical box model BOXMOX • genbox: input data generator for boxmox • frappedata: FRAPPE campaign dataset for genbox • tuv: TUV data connector for genbox • chemspectranslator: translator to translate species between chemical mechanisms and observations 10 • icartt: reader/writer for ICARTT files The underlying chemical box model BOXMOX is a standalone C/Fortran executable and has to be installed before BEATBOX can be used.
Source code version 1.0 used in this manuscript is archived online at digital object identifier (DOI) 10.5281/zenodo.1118244 (Knote et al., 2017). Full documentation for BOXMOX and all Python packages, including examples on how to reproduce 15 the case studies shown in this manuscript can be found at https://boxmodeling.meteo.physik.unimuenchen.de/documentation . 20 6 Author contributions 25 CK and JB contributed equally to the manuscript. CK designed the BOXMOX system, wrote the framework around BEATBOX, supervised ME, and contributed to writing the manuscript. JB came up with the idea for BEATBOX and coded most of the beatbox.py routine, co-supervised ME and wrote most of the manuscript. ME developed parts of the BEATBOX framework during his Masters thesis, prepared and analysed the case studies, and contributed the flux tool.
We thank the many scientists who contributed to the FRAPPE field campaign for providing useful data to initialize the model simulations. Frank Flocke and Rebecca Hornbrook (NCAR) are thanked for creating the 'Commerce City' sample used in the case studies. We acknowledge Bill Carter (UC Riverside) for his work on the emission database. We acknowledge support from the NASA KORUS-AQ grant NNX16AD96G. 5

Competing interests
The authors declare that they have no conflict of interest.