3.1 DATeS architecture

The DATeS architecture abstracts, and provides a set of modules of, the four generic components of any DA system. These components are the linear algebra routines, a forecast computer model that includes the discretization of the physical processes, error models, and analysis methodologies. In what follows, we discuss each of these building blocks in more detail, in the context of DATeS. We start with an abstract discussion of each of these components, followed by technical descriptions.

3.1.6 DATeS layout

The design of DATeS takes into account the distinction between these components and separates them in design following an object-oriented programming (OOP) approach. A general description of DATeS architecture is given in Fig. 1.

The enumeration in Fig. 1 (numbers from 1 to 4 in circles) indicates the order in which essential DATeS objects should be created. Specifically, one starts with an instance of a model. Once a model object is created, an assimilation object is instantiated, and the model object is passed to it. An assimilation process object is then instantiated, with a reference to the assimilation object passed to it. The assimilation process object iterates the consecutive assimilation cycles and saves and/or outputs the results which can be optionally analyzed later using visualization modules.

All DATeS components are independent so as to maximize the flexibility in experimental design. However, each newly added component must comply with DATeS rules in order to guarantee interoperability with the other pieces in the package. DATeS provides base classes with definitions of the necessary methods. A new class added to DATeS, for example, to implement a specific new model, has to inherit the appropriate model base class and provide implementations of the inherited methods from that base class.

Figure 1Diagram of the DATeS architecture.

Download

In order to maximize both flexibility and generalizability, we opted to handle configurations, inputs, and output of DATeS object using “configuration dictionaries”. Parameters passed to instantiate an object are passed to the class constructor in the form of key-value pairs in the dictionaries. See Sect. 4 for examples on how to properly configure and instantiate DATeS objects.

3.2 Linear algebra classes

The main linear algebra data structures essential for almost all DA aspects are (a) model state-size and observation-size vectors (also named state and observation vectors, respectively), and (b) state-size and observation-size matrices (also named state and observation matrices, respectively). A state matrix is a square matrix of order equal to the model state-space dimension. Similarly, an observation matrix is a square matrix of order equal to the model observation space dimension. DATeS makes a distinction between state and observation linear algebra data structures. It is important to recall here that, in large-scale applications, full state covariance matrices cannot be explicitly constructed in memory. Full state matrices should only be considered for relatively small problems and for experimental purposes. In large-scale settings, where building state matrices is infeasible, low-rank approximations or sparse representation of the covariance matrices could be incorporated. DATeS provides simple classes to construct sparse state and observation matrices for guidance.

Table 1DA filtering routines provided by the initial version of DATeS (v1.0).

Download Print Version | Download XLSX

Third-party linear algebra routines can have widely different interfaces and underlying data structures. For reusability, DATeS provides unified interfaces for accessing and manipulating these data structures using Python classes. The linear algebra classes are implemented in Python. The functionalities of the associated methods can be written either in Python or in lower-level languages using proper wrappers. A class for a linear algebra data structure enables updating, slicing, and manipulating an instance of the corresponding data structures. For example, a model state vector class provides methods that enable accessing/slicing and updating entries of the state vector, a method for adding two state vector instances, and methods for applying specific scalar operations on all entries of the state vector such as evaluating the square root or the logarithm. Once an instance of a linear algebra data structure is created, all its associated methods are accessible via the standard Python dot operator. The linear algebra base classes provided in DATeS are summarized in Table 1.

Figure 2Python implementation of state vector, observation vector, state matrix, and observation matrix data structures. Both dense and sparse state and observation matrices are provided.

Download

Python special methods are provided in a linear algebra class to enable iterating a linear algebra data structure entries. Examples of these special methods include__getitem__( ), __setitem__( ), __getslice__( ), __setslice__( ), etc. These operators make it feasible to standardize working with linear algebra data structures implemented in different languages or saved in memory in different forms.

DATeS provides linear algebra data structures represented as NumPy ndarrays, and a set of NumPy-based classes to manipulate them. Moreover, SciPy-based implementation of sparse matrices is provided and can be used efficiently in conjunction with both sparse and non-sparse data structures. These classes, shown in Fig. 2, provide templates for designing more sophisticated extensions of the linear algebra classes.

Table 2DA filtering routines provided by the initial version of DATeS (v1.0).

Download Print Version | Download XLSX

3.3 Forecast model classes

Each numerical model needs an associated class providing methods to access its functionality. The unified forecast model class design in DATeS provides the essential tasks that can be carried out by the model implementation. Each model class in DATeS has to inherit the model base class: models_base.ModelBase or a class derived from it. A numerical model class is required to provide access to the underlying linear algebra data structures and time integration routines. For example, each model class has to provide the method state_vector( ) that creates an instance of a state vector class, and the method integrate_state( ) that takes a state vector instance and time integration settings, and returns a trajectory (list of states) evaluated at specified future times. The base class provided in DATeS contains definitions of all the methods that need to be supported by a numerical model class. The package DATeS v1.0 includes implementations of several popular test models summarized in Table 2.

While some linear algebra and the time integration routines are model-specific, DATeS also implements general-purpose linear algebra classes and time integration routines that can be reused by newly created models. For example, the general integration class FatODE_ERK_FWD is based on FATODE (Zhang and Sandu, 2014) explicit Runge–Kutta (ERK) forward propagation schemes.

3.4 Error model classes

In many DA applications, the errors are additive and are modeled by random variables normally distributed with zero mean and a given or an unknown covariance matrix. DATeS implements NumPy-based functionality for background, observation, and model errors as guidelines for more sophisticated problem-dependent error models. The NumPy-based error models in DATeS are implemented in the module error_models_numpy. These classes are derived from the base class ErrorsModelBase and provide methodologies to sample the underlying probability distribution, evaluate the value of the density function, and generate statistics of the error variables based on model trajectories and the settings of the error model. Note that, while DATeS provides implementations for Gaussian error models, the Gaussian assumption itself is not restrictive. Following the same structure, or by inheritance, one can easily create non-Gaussian error models with minimal efforts. Moreover, the Gaussian error models provided by DATeS support both correlated and uncorrelated errors, and it constructs the covariance matrices accordingly. The covariance matrices are stored in appropriate sparse formats, unless a dense matrix is explicitly requested. Since these covariance matrices are either state or observation matrices, they provide access to all proper linear algebra routines. This means that the code written with access to an observation error model and its components should work for both correlated and uncorrelated observations.

Table 3DA filtering routines provided by the initial version of DATeS v1.0.

Download Print Version | Download XLSX

3.5 Assimilation classes

Assimilation classes are responsible for carrying out a single assimilation cycle (i.e., over one assimilation window) and optionally printing or writing the results to files. For example, an EnKF object should be designed to carry out one cycle consisting of the “forecast” and the “analysis” steps. The basic assimilation objects in DATeS are a filtering object, a smoothing object, and a hybrid object. DATeS provides the common functionalities for filtering objects in the base class filters_base.FiltersBase; all derived filtering classes should have it as a super class. Similarly, smoothing objects are to be derived from the base class smoothers_base.SmoothersBase. A hybrid object can inherit methods from both filtering and smoothing base classes.

A model object is passed to the assimilation object constructor via configuration dictionaries to give the assimilation object access to the model-based data structures and functionalities. The settings of the assimilation object, such as the observation time, the assimilation time, the observation vector, and the forecast state or ensemble, are also passed to the constructor upon instantiation and can be updated during runtime.

Table 3 summarizes the filters implemented in the initial version of the package, which is DATeS v1.0. Each of these filtering classes can be instantiated and run with any of the DATeS model objects. Moreover, DATeS provides simplified implementations of both 3D-Var and 4D-Var assimilation schemes. The objective function, e.g., the negative log posterior, and the associated gradient are implemented inside the smoother class and require the tangent linear model to be implemented in the passed forecast model class. The adjoint is evaluated using FATODE following a checkpointing approach, and the optimization step is carried out using SciPy optimization functions. The settings of the optimizer can be fine-tuned via the configuration dictionaries. The 3D- and 4D-Var implementations provided by DATeS are experimental and are provided as a proof of concept. The variational aspects of DATeS are being continuously developed and will be made available in future releases of the package.

Covariance inflation and localization are ubiquitously used in all ensemble-based assimilation systems. These two methods are used to counteract the effect of using ensembles of finite size. Specifically, covariance inflation counteracts the loss of variance incurred in the analysis step and works by inflating the ensemble members around their mean. This is carried out by magnifying the spread of ensemble members around their mean by a predefined inflation factor. The inflation factor could be a scalar, i.e., space–time independent, or even varied over space and/or time. Localization, on the other hand, mitigates the accumulation of long-range spurious correlations. Distance-based covariance localization is widely used in geoscientific sciences, and applications, where correlations are damped out with increasing distance between grid points. The performance of the assimilation algorithm is critically dependent on tuning the parameters of these techniques. DATeS provide basic utility functions (see Sect. 3.7) for carrying out inflation and localization which can be used in different forms based on the specific implementation of the assimilation algorithms. The work in Attia and Constantinescu (2018) reviews inflation and localization and presents a framework for adaptive tuning of the parameters of these techniques, with all implementations and numerical experiments carried out entirely in DATeS.

Figure 3The assimilation process in DATeS.

Download

3.6 Assimilation process classes

A common practice in sequential DA experimental settings is to repeat an assimilation cycle over a given time span, with similar or different settings at each assimilation window. For example, one may repeat a DA cycle on several time intervals with different output settings, e.g., to save and print results only every fixed number of iterations. Alternatively, the DA process can be repeated over the same time interval with different assimilation settings to test and compare results. We refer to this procedure as an “assimilation process”. Examples of numerical comparisons, carried out using DATeS, can be found in Attia et al. (2018) and Attia and Constantinescu (2018), and in Sect. 4.6.

assimilation_process_base.Assimilation Process is the base class from which all assimilation process objects are derived. When instantiating an assimilation process object, the assimilation object, the observations, and the assimilation time instances are passed to the constructor through configuration dictionaries. As a result, the assimilation process object has access to the model and its associated data structures and functionalities through the assimilation object.

The assimilation process object either retrieves real observations or creates synthetic observations at the specified time instances of the experiment. Figure 3 summarizes DATeS assimilation process functionality.

3.7 Utility modules

Utility modules provide additional functionality, such as the _utility_configs module which provides functions for reading, writing, validating, and aggregating configuration dictionaries. In DA, an ensemble is a collection of state or observation vectors. Ensembles are represented in DATeS as lists of either state or observation vector objects. The utility modules include functions responsible for iterating over ensembles to evaluate ensemble-related quantities of interest, such as ensemble mean, ensemble variance/covariance, and covariance trace. Covariance inflation and localization are critically important for nearly all ensemble-based assimilation algorithms. DATeS abstracts tools and functions common to assimilation methods, such as inflation and localization, where they can be easily imported and reused by newly developed assimilation routines. The utility module in DATeS provides methods to carry out these procedures in various modes, including state-space and observation space localization. Moreover, DATeS supports space-dependent covariance localization; i.e., it allows varying the localization radii and inflation factors over both space and time.

Table 4A sample of the modules wrapped by the main utility module dates_utility.

Download Print Version | Download XLSX

Ensemble-based assimilation algorithms often require matrix representation of ensembles of model states. In DATeS, ensembles are represented as lists of states, rather than full matrices of size N_state×N_ens. However, it provides utility functions capable of efficiently calculating ensemble statistics, including ensemble variances, and covariance trace. Moreover, DATeS provides matrix-free implementations of the operations that require ensembles of states, such as a matrix–vector product, where the matrix is involved in a representation of an ensemble of states.

Figure 4The sequence of essential steps required in order to run a DA experiment in DATeS.

Download

The module dates_utility provides access to all utility functions in DATeS. In fact, this module wraps the functionality provided by several other specialized utility routines, including the sample given in Table 4. The utility module provides other general functions such as handling file downloading, and functions for file I/O. For a list of all functions in the utility module, see the user manual (Attia et al., 2016).

Figure 5Initializing the DATeS run.

Download

4.1 Step1: initialize DATeS

Initializing a DATeS run involves defining the root directory of DATeS as an environment variable and adding the paths of DATeS source modules to the system path. This can be done by executing the code snippet in Fig. 5 in the DATeS root directory.

4.2 Step2: create a model object

QG-1.5 is a nonlinear 1.5-layer reduced-gravity QG model with double-gyre wind forcing and biharmonic friction (Sakov and Oke, 2008).

4.2.1 Quasi-geostrophic model

This model is a numerical approximation of the equations

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{q}_t={\mathit{\psi }}_x-\mathit{\epsilon }J(\mathit{\psi },q)-A{\mathrm{\Delta }}^{\mathrm{3}}\mathit{\psi }+\mathrm{2}\mathit{\pi }\sin \left(\mathrm{2}\mathit{\pi }y\right),\\ {\displaystyle }& {\displaystyle }q=\mathrm{\Delta }\mathit{\psi }-F\mathit{\psi },\\ \text{(9)}& {\displaystyle }& {\displaystyle }J(\mathit{\psi },q)\equiv {\mathit{\psi }}_x{q}_x-{\mathit{\psi }}_y{q}_y,\end{array}$$

where $\mathrm{\Delta }:={\partial }^{\mathrm{2}}/\partial x^{\mathrm{2}}+{\partial }^{\mathrm{2}}/\partial y^{\mathrm{2}}$ and ψ is the surface elevation. The values of the model coefficients in Eq. (9) are obtained from Sakov and Oke (2008) and are described as follows: F=1600, $\mathit{\epsilon }={\mathrm{10}}^{-\mathrm{5}}$, and $A=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{12}}$. The domain of the model is a 1×1 (space units) square, with $\mathrm{0}\le x\le \mathrm{1}$, $\mathrm{0}\le y\le \mathrm{1}$, and is discretized by a grid of size 129×129 (including boundaries). The boundary conditions used are $\mathit{\psi }=\mathrm{\Delta }\mathit{\psi }={\mathrm{\Delta }}^{\mathrm{2}}\mathit{\psi }=\mathrm{0}$. The dimension of the model state vector is N_state=16 641. This is a synthetic model where the scales are not relevant, and we use generic space, time, and solution amplitude units. The time integration scheme used is the fourth-order Runge–Kutta scheme with a time step of 1.25 (time units). The model forward propagation core is implemented in Fortran. The QG-1.5 model is run over 1000 model time steps, with observations made available every 10 time steps.

Figure 6Creating the QG model object.

Download

4.2.2 Observations and observation operators

We use a standard linear operator to observe 300 components of ψ with observation error variance set to 4.0 (units squared). The observed components are uniformly distributed over the state vector length, with an offset that is randomized at each filtering cycle. Synthetic observations are obtained by adding white noise to measurements of the sea surface height level (SSH) extracted from a reference model run with lower viscosity. To create a QG model object with these specifications, one executes the code snippet in Fig. 6.

Figure 7Creating a DEnKF filtering object.

Download

4.3 Step3: create an assimilation object

One now proceeds to create an assimilation object. We consider a deterministic implementation of EnKF (DEnKF) with ensemble size equal to 20, and parameters tuned optimally as suggested in Sakov and Oke (2008). Covariance localization is applied via a Hadamard product (Houtekamer and Mitchell, 2001). The localization function is Gaspari–Cohn (Gaspari and Cohn, 1999) with a localization radius of 12 grid cells. The localization is carried out in the observation space by decorrelating both HB and HBH^T, where B is the ensemble covariance matrix, and H is the linearized observation operator. In the present setup, the observation operator ℋ is linear, and thus H=ℋ.

Ensemble inflation is applied to the analysis ensemble of anomalies at the end of each assimilation cycle of DEnKF with an inflation factor of 1.06. The code snippet in Fig. 7 creates a DEnKF filtering object with these settings.

Figure 8Creating a filter process object to carry out DEnKF filtering using the QG model.

Download

Most of the methods associated with the DEnKF object will raise exceptions if immediately invoked at this point. This is because several keys in the filter configuration dictionary, such as the observation, the forecast time, the analysis time, and the assimilation time, are not yet appropriately assigned. DATeS allows creating assimilation objects without these options to maximize flexibility. A convenient approach is to create an assimilation process object that, among other tasks, can properly update the filter configurations between consecutive assimilation cycles.

Figure 9Running the filtering experiment.

Download

4.4 Step4: create an assimilation process

We now test DEnKF with the QG model by repeating the assimilation cycle over a time span from 0 to 1250 with offsets of 12.5 time units between each two consecutive observation/assimilation times. An initial ensemble is created by the numerical model object. An experimental time span is set for observations and assimilation. Here, the assimilation time instances da_checkpoints are the same as the observation time instances obs_checkpoints, but they can in general be different, leading to either synchronous or asynchronous assimilation settings. This is implemented in the code snippet in Fig. 8.

Here, experiment is responsible for creating synthetic observations at all time instances defined by obs_checkpoints (except the initial time). To create synthetic observations, the truth at the initial time (0 in this case) is obtained from the model and is passed to the filtering process object experiment, which in turn propagates it forward in time to assimilation time points.

Finally, the assimilation experiment is executed by running the code snippet in Fig. 9.

4.5 Experiment results

The filtering results are printed to screen and are saved to files at the end of each assimilation cycle as instructed by the output_configs dictionary of the object experiment. The output directory structure is controlled via the options in the output configuration dictionary output_configs of the FilteringProcess object, i.e., experiment. All results are saved in appropriate subdirectories under a main folder named Results in the root directory of DATeS. We will refer to this directory henceforth as DATeS results directory. The default behavior of a FilteringProcess object is to create a folder named Filtering_Results in the DATeS results directory and to instruct the filter object to save/output the results every file_output_iter whenever the flag file_output is turned on. Specifically, the DEnKF object creates three directories named Filter_Statistics, Model_States_Repository, and Observations_Repository, respectively. The root mean square (rms) forecast and analysis errors are evaluated at each assimilation cycle and are written to a file under the Filter_Statistics directory. The output configurations of the filter object of the DEnKF class, i.e., denkf_filter, instruct the filter to save all ensemble members (both forecast and analysis) to files by setting the value of the option file_output_moment_only to False. The true solution (reference state), the analysis ensemble, and the forecast ensembles are all saved under the directory Model_States_Repository, while the observations are saved under the directory Observations_Repository. We note that, while here we illustrate the default behavior, the output directories are fully configurable.

Figure 10 shows the reference initial state of the QG model, an example of the observational grid used, and an initial forecast state. The initial forecast state in Fig. 10 is the average of an initial ensemble collected from a long run of the QG model.

The true field, the forecast errors, and the DEnKF analyses errors at different time instances are shown in Fig. 11.

Typical solution quality metrics in the ensemble-based DA literature include RMSE plots and rank (Talagrand) histograms (Anderson, 1996; Candille and Talagrand, 2005).

Upon termination of a DATeS run, executable files can be cleaned up by calling the function clean_executable_files( ) available in the utility module (see the code snippet in Fig. 12).

Figure 10The QG-1.5 model. The truth (reference state) at the initial time (t=0) of the assimilation experiment is shown in panel (a). The red dots indicate the locations of observations for one of the test cases employed. The initial forecast state, taken as the average of the initial ensemble at time t=0, is shown in panel (b).

Download

4.6 DATeS for benchmarking

4.6.1 Performance metrics

In the linear settings, the performance of an ensemble-based DA filter could be judged based on two factors. Firstly, convergence is explained by its ability to track the truth and secondly by the quality of the flow-dependent covariance matrix generated given the analysis ensemble.

Figure 11Data assimilation results. The reference field ψ, the forecast errors, and the analysis errors at t=300, t=600, t=900, and t=1200 (time units) are shown. Here, the forecast error is defined as the reference field minus the average of the forecast ensemble, and the analysis error is the reference field minus the average of the analysis ensemble.

Download

Figure 12Cleanup of DATeS executable files.

Download

Figure 13Data assimilation results. In panel (a), “no assimilation” refers to the RMSE of the initial forecast (the average of the initial forecast ensemble) propagated forward in time over the 100 cycles without assimilating observations into it. The rank histogram of where the truth ranks among analysis ensemble members is shown in panel (b). The ranks are evaluated for every 13th variable in the state vector (past the correlation bound) after 100 assimilation cycles.

Download

The convergence of the filter is monitored by inspecting the RMSE, which represents an ensemble-based standard deviation of the difference between reality, or truth, and the model-based prediction. In synthetic experiments, where the model representation of the truth is known, the RMSE reads

$$\begin{array}{}\text{(10)}& \mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{state}}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{state}}}(x_i-x_i^{\mathrm{True}}{)}^{\mathrm{2}}},\end{array}$$

where $\mathbf{x}=(x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{\dots },x_{N_{\mathrm{state}}}{)}^T\in {\mathbb{R}}^{N_{\mathrm{state}}}$ is the prediction at a given time instant, e.g., the forecast ensemble mean, and ${\mathbf{x}}^{\mathrm{True}}\in {\mathbb{R}}^{N_{\mathrm{state}}}$ is the verification, e.g., the true model state at the same time instant. For real applications, the states are generally replaced with observations. The rank (Talagrand) histogram (Anderson, 1996; Candille and Talagrand, 2005) could be used to assess the spread of the ensemble and its coverage to the truth. Generally speaking, the rank histogram plots the rank of the truth (or observations) compared to the ensemble members (or equivalent observations), ordered increasingly in magnitude. A nearly uniform rank histogram is desirable and suggests that the truth is indistinguishable from the ensemble members. A mound rank histogram indicates an overdispersed ensemble, a while U-shaped histogram indicates underdispersion. However, mound rank histograms are rarely seen in practice, especially for large-scale problems. See, e.g., Hamill (2001) for a mathematical description and a detailed discussion on the usefulness and interpretation of rank histograms.

Figure 13a shows an RMSE plot of the results of the experiment presented in Sect. 4.5. The histogram of the rank statistics of the truth, compared to the analysis ensemble, is shown in Fig. 13b.

For benchmarking, one needs to generate scalar representations of the RMSE and the uniformity of a rank histogram of a numerical experiment. The average RMSE can be used to compare the accuracy of a group of filters. To generate a scalar representation of the uniformity of a rank histogram, we fit a beta distribution to the rank histogram, scaled to the interval [0,1], and evaluate the Kullback–Leibler (KL) divergence (Kullback and Leibler, 1951) between the fitted distribution and a uniform distribution. The KL divergence between two beta distributions Beta(α,β), and $\mathrm{Beta}({\mathit{\alpha }}^{\prime },{\mathit{\beta }}^{\prime })$ is $D_{\mathrm{KL}}\left(\mathrm{Beta}(\mathit{\alpha },\mathit{\beta })\left|\mathrm{Beta}\right({\mathit{\alpha }}^{\prime }{\mathit{\beta }}^{\prime })\right)=\ln \mathrm{\Gamma }(\mathit{\alpha }+\mathit{\beta })-\ln \left(\mathit{\alpha }\mathit{\beta }\right)-\ln \mathrm{\Gamma }({\mathit{\alpha }}^{\prime }+{\mathit{\beta }}^{\prime })+\ln \left({\mathit{\alpha }}^{\prime }{\mathit{\beta }}^{\prime }\right)+(\mathit{\alpha }-{\mathit{\alpha }}^{\prime })\left(\mathit{\psi }\left(\mathit{\alpha }\right)-\mathit{\psi }\left({\mathit{\alpha }}^{\prime }\right)\right)+(\mathit{\beta }-{\mathit{\beta }}^{\prime })\left(\mathit{\psi }\left(\mathit{\beta }\right)-\mathit{\psi }\left({\mathit{\beta }}^{\prime }\right)\right)$, where $\mathit{\psi }(\cdot )={\mathrm{\Gamma }}^{\prime }(\cdot )/\mathrm{\Gamma }(\cdot )$ is the digamma function, i.e., the logarithmic derivative of the gamma function. Here, we set $\mathrm{Beta}({\mathit{\alpha }}^{\prime },{\mathit{\beta }}^{\prime })$ to a uniform distribution by setting ${\mathit{\alpha }}^{\prime }={\mathit{\beta }}^{\prime }=\mathrm{1}$. We consider a small, e.g., closer to 0, KL distance to be an indication of a nearly uniform rank histogram and consequently an indication of a well-dispersed ensemble. An alternative measure of rank histogram uniformity is to average the absolute distances of bins' heights from a uniformly distributed rank histogram (Bessac et al., 2018). DATeS provides several utility functions to calculate such metrics for a numerical experiment.

Figure 14Rank histograms with fitted beta distributions. The KL-divergence measure is indicated under each panel.

Download

Table 5Measures of uniformity of the rank histograms shown in Fig. 14.

Download Print Version | Download XLSX

Figure 14 shows several rank histograms, along with uniform distribution, and fitted beta distributions. The KL-divergence measure is indicated under each panel. Results in Fig. 14 suggest that the fitted beta distribution parameters give, in most cases, a good scalar description of the shape of the histogram. Moreover, one can infer the shape of the rank histogram from the parameters of the fitted beta distribution. For example, if α>1 and β>1, the histogram has a mound shape, and is U-shaped if α<1 and β<1. The histogram is closer to uniformity as the parameters α, β both approach 1. Table 5 shows both KL distance between fitted beta distribution with respect to a uniform one, and the average distances between histogram bins and a uniform one.

Figure 15Data assimilation results with DEnKF applied to Lorenz-96 system. RMSE results on a log scale are shown in panel (a). The KL distances between the analysis rank histogram and a uniform rank histogram are shown in panel (b). The localization radius is fixed to 4.

Download

4.6.2 Benchmarking

The architecture of DATeS makes it easy to generate benchmarks for a new experiment. For example, one can write short scripts to iterate over a combination of settings of a filter to find the best possible results. As an example, consider the standard 40-variable Lorenz-96 model (Lorenz, 1996) described by the equations

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}{x}_i}{\mathrm{d}t}}=x_{i-\mathrm{1}}\left(x_{i+\mathrm{1}}-x_{i-\mathrm{2}}\right)-x_i+F;i=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{40},\end{array}$$

where $\mathbf{x}=(x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{\dots },x_{\mathrm{40}}{)}^T\in {\mathbb{R}}^{\mathrm{40}}$ is the state vector, with periodic boundaries, i.e., x₀≡x₄₀, and the forcing parameter is set to F=8. These settings make the system chaotic (Lorenz and Emanuel, 1998) and are widely used in synthetic settings for geoscientific applications. Adjusting the inflation factor and the localization radius for EnKF filter is crucial. Consider the case where one is testing an adaptive inflation scheme and would like to decide on the ensemble size and the benchmark inflation factor to be used. As an example of benchmarking, we run the following experiment over a time interval [0,30] (units), where Eq. (11) is integrated forward in time using a fourth-order Runge–Kutta scheme with model step size 0.005 (units). Assume that synthetic observations are generated every 20 model steps, where every other entry of the model state is observed. We test the DEnKF algorithm, with the fifth-order piecewise-rational function of Gaspari and Cohn (Gaspari and Cohn, 1999) for covariance localization. The localization radius is held constant and is set to l=4, while the inflation factor is varied for each experiment. The experiments are repeated for ensemble sizes $N_{\mathrm{ens}}=\mathrm{5},\mathrm{10},\mathrm{\dots },\mathrm{40}$. We report the results over the last two-thirds of the experiments' time span, i.e., over the interval [10,30], to avoid spinup artifacts. This interval, consisting of the last 200 assimilation cycles out of 300, will be referred to as the “testing time span”. Any experiment that results in an average RMSE of more than 0.65 over the testing time span is discarded, and the filter used is seen to diverge. The numerical results are summarized in Figs. 15 and 16.

Figure 15 shows the average RMSE results and the KL distances between a beta distribution fitted to the analysis rank histogram of each experiment, and a uniform distribution. These plots give a preliminary idea of the plausible regimes of both ensemble size and inflation factor that should be used to achieve the best performance of the filter used under the current experimental settings. For example, for an ensemble size N_ens=20, the inflation factor should be set approximately to 1.01–1.07 to give both a small RMSE and an analysis rank histogram close to uniform.

Figure 16Data assimilation results with DEnKF applied to Lorenz-96 system. The minimum average RMSE over the interval [10,30] is indicated by red triangles. The minimum average KL distance between the analysis rank histogram and a uniformly distributed rank histogram [10,30] is indicated by blue tripods. We show the results for every choice of the ensemble size. The localization radius is fixed to 4.

Download

Concluding the best inflation factor for a given ensemble size, based on Fig. 15, however, could be tricky. Figure 16 shows the inflation factors resulting in minimum average RMSE and minimum KL distance to uniformity. Specifically, for each ensemble size, a red triangle refers to the experiment that resulted in minimum average RMSE over the testing time span, out of all benchmarking experiments carried out with this ensemble size. Similarly, the experiment that yielded minimum KL divergence to a uniform rank histogram is indicated by a blue tripod.

To answer the question about the ensemble size, we pick the ensemble size N_ens=25, given the current experimental setup. The reason is that N_ens=25 is the smallest ensemble size that yields small RMSE and is a well-dispersed ensemble as explained by Fig. 16. As for the benchmark inflation factor, the results in Fig. 16 show that for an ensemble size N_ens=25, the best choice of an inflation factor is approximately 1.03–1.05 for Gaspari–Cohn localization with a fixed radius of 4.

Despite being a relatively easy process, unfortunately, generating a set of benchmarks for all possible combinations of numerical experiments is a time-consuming process and is better carried out by the DA community. Some example scripts for generating and plotting benchmarking results are included in the package for guidance.

Note that, when the Gaussian assumption is severely violated, standard benchmarking tools, such as RMSE and rank histograms, should be replaced with, or at least supported by, tools capable of assessing ensemble coverage of the posterior distribution. In such cases, MCMC methods, including those implemented in DATeS (Attia and Sandu, 2015; Attia et al., 2018; Attia, 2016), could be used as a benchmarking tool (Law and Stuart, 2012).

5.1 Adding a numerical model class

A new model class has to be created as a subclass of ModelsBase or a class derived from it. The base class ModelsBase, similar to all base classes in DATeS, contains headers of all the functions that need to be provided by a model class to guarantee that it interacts properly with other components in DATeS.

Figure 18The leading lines of an implementation of a class for the model MyModel derived from the models' base class ModelsBase. Linear algebra objects are derived from NumPy-based (or SciPy-based) objects.

Download

The first step is to grant the model object access to linear algebra data structures and to error models. Appropriate classes should be imported in a numerical model class:

• Linear algebra includes state vector, state matrix, observation vector, and observation matrix.

• Error models include background, model, and observation error models.

This gives the model object access to model-based data structures and error entities necessary for DA applications. Figure 17 illustrates a class of a numerical model named MyModel, along with all the essential classes imported by it.

The next step is to create Python-based implementations for the model functionalities. As shown in Fig. 17, the corresponding methods have descriptive names in order to ease the use of DATeS functionality. For example, the method state_vector( ) creates (or initializes) a state vector data structure. Details of each of the methods in Fig. 17 are given in the DATeS user manual (Attia et al., 2016).

As an example, suppose we want to create a model class name MyModel using NumPy and SciPy (for sparse matrices) linear algebra data structures. The code snippet in Fig. 18 shows the implementation of such a class.

Note that in order to guarantee extensibility of the package we have to fix the naming of the methods associated with linear algebra classes, and even if only binary files are provided, the Python-based linear algebra methods must be implemented. If the model functionality is fully written in Python, the implementation of the methods associated with a model class is straightforward, as illustrated in Attia et al. (2016). On the other hand, if a low-level implementation of a numerical model is given, these methods wrap the corresponding low-level implementation.

5.2 Adding an assimilation class

The process of adding a new class for an assimilation methodology is similar to creating a class for a numerical model; however, it is expected to require less effort. For example, a class implementation of a filtering algorithm uses components and tools provided by the passed model and by the encapsulated linear algebra data structures and methods. Moreover, filtering algorithms belonging to the same family, such as different flavors of the well-known EnKF, are expected to share a considerable amount of infrastructure. Python inheritance enables the reuse of methods and variables from parent classes.

Figure 19Illustration of a DA filtering class MyFilter and its relation to the filtering base class. A solid arrow represents an “inherit” relation.

Download

To create a new class for DA filtering, one derives it from the base class FiltersBase, imports appropriate routines, and defines the necessary functionalities. Note that each assimilation object has access to a model object and consequently to the proper linear algebra data structures and associated functionalities through that model.

Figure 20The leading lines of an implementation of a DA filter; the MyFilter class is derived from the filters' base class FiltersBase.

Download

Unlike the base class for numerical models (ModelsBase), the filtering base class FiltersBase includes actual implementations of several widely used solvers. For example, an implementation of the method FiltersBase.filtering_cycle( ) is provided to carry out a single filtering cycle by applying a forecast phase followed by an analysis phase (or vice versa, depending on stated configurations).

Figure 19 illustrates a filtering class named MyFilter that works by carrying out analysis and forecast steps in the ensemble-based statistical framework.

The code snippet in Fig. 20 shows the leading lines of an implementation of the MyFilter class.