\svnidlong
{$HeadURL: $}
{$LastChangedDate: 2017-08-28 18:10:04 +0200 (Tue, 08 Jul 2014) $}
{$LastChangedRevision: $}
{$LastChangedBy:$}

%opening
\odachapter{OpenDA: How to setup and test your application step by step}


\begin{tabular}{p{4cm}l}
\textbf{Contributed by:} & Nils van Velzen, \vortech, Martin Verlaan, Deltares \\
\textbf{Last update:}    & \svnfilemonth-\svnfileyear\\
\end{tabular}



\section{Introduction}
Setting up a data assimilation framework for a model is a difficult task. Several things contribute to the complexity. The dynamical models are often complex software packages with many options. In addition, we add (real) observed data and blow up the amount of computations and data by an order of 100. 
Another problem is related to our way of using the models in a data assimilation framework. For sequential data assimilation algorithms, such as the Ensemble Kalman Filter or 3D-VAR, We often perform (short) model runs and update the parameters or state of the model between these short runs. This is a way of using the simulation models, for which they are often not developed.

We have seen so often that students/users are struggling to get things to work, because they want to do too much too soon. A recipe for failure is to attempt to setup your data assimilation system for a real (big/huge) model  with real data in one go. So it is clear that you should take some intermediate steps, but which? There is no one-size-fits-all approach but in this document we try to present a recipe, which you can use or adapt to your own needs. In between there are tips/ideas, which are hopefully useful to you.

For simplicity we assume that the user uses a black-box coupling and wants to setup a data assimilation system using a sequential data assimilation algorithm e.g. a flavor of the ensemble Kalman filter.

\section{Preparation}
\subsection{Start small}
Setting up a data assimilation system involves many steps and challenges. It is not advised to directly focus on your final setup, which will involve real observations and probably a large simulation model. The last thing you will introduce are your real observations. You need their specifications e.g. location, quantity, quality and sampling rate in an early stage since it is an important aspect of your system but the measured values will not be used the first 80\% of the time. In order to setup and test your framework it is best to use generate/synthetic data which you understand. This will be explained in more detail, when we explain about setting up twin experiments in Section \ref{Sec:Twin}.

If possible, make various variations of your model. Start with a very, very simplified model that runs blisteringly fast and only incorporates the most basic features. When you have everything working for this small model, you move towards a more complex model. The amount of steps you have to take towards your full model depends on many aspects. Making these "extra" small steps is time well spent and in our experience you will save a lot of time in the end. 

Create experiments with one group of observations at a time when you want to assimilate observations of various quantities and/or sources . You will learn a lot about the behavior of your model when assimilating these different types of observations and it is much easier to identify which kind of observations might cause problems, like model instabilities.

\subsection{Check the restart of a model}
In order to use a model\footnote{Here we refer to the model as the program, not the mathematical description or the set-up for a particular application/area} as part of a sequential data assimilation algorithm, it needs to have a proper restart functionality. This makes it possible to split up a long simulation run into several shorter ones. The model will write the internal state to a restart file or files, at the end of each run. This will contain the model state $x$, but often some other information as well e.g. the information on the integration step size, computed forcing, etc. The restart information will be read from disk at the start of the next run. There should be no differences in the result between the restarted simulations and the original simulation when the restart is implemented correctly. 

The first step is to check whether the restart functionality of your model is working properly. Run a simulation in one go and perform the same simulation with a number of restarts. It is always best to choose the same interval between the restarts as your assimilation interval you are planning to use in your data assimilation framework.

Note: When your model is already available in OpenDA, you can use OpenDA to do this experiment (will be explained later), but do not skip this step, because your model configuration might contain features that have not been used before and for which the restart might be faulty!

%TODO add a figure here on restart checking

Unfortunately, the restart functionality of models is often not perfect. When that is the case you have to look how bad it is. Here is a list of issues we have seen in the past that might cause differences:
\begin{itemize}
\item loss of precision: It can happen that some precision is lost in reading and writing values from the restart files (e.g. computations are in double precision but restart is in single precision). When we expect that the model updates of the data assimilation algorithm are much larger than this loss of precision, it is only annoying (it makes testing/comparing/debugging more difficult), but no show stopper.
\item Incomplete restart information: At some point in the history of the model some functionality has been added but the developers forgot to incorporate the relevant new (state) information in the restart file.
\item Imperfect by design: Sometimes the developers never intended to have a perfect restart functionality, which means the results are not exactly same as without the restart. (Writing a correct restart functionality is in many cases far from easy)
\end{itemize}

Some tips when you notice the restart is imperfect:
\begin{itemize}
\item Experiment with a simplified model. Switch features on and off to figure out where the differences are originating from.
\item Does your model have automatic integration steps? Check the initial integration time steps for your restarted model. Can you run your model with constant integration time steps?
\item How is the model forcing defined? Does the model interpolate your forcing input data? Changing the model time steps might fix your problems.
\item Contact the developers of the code. With some luck, they are willing to help you.
\end{itemize}

In the end you have to figure out whether the errors in the restart are acceptably small. When the deviation between the original run and a run with restarts is much smaller than the expected impact of your data assimilation you might be OK.

\subsection{Uncertainty of your model}
For the ensemble based algorithms, you need to have an ensemble that statistically represents the uncertainty in your model prediction. There are various ways to setup your ensemble. 

When your model is dominated by chaotic behavior, e.g. for most ocean models and atmospheric models, you can generate an initial ensemble by running the model for some time and taking various snapshots of the state. Another approach is to setup an ensemble with some initial perturbation. Then run the ensemble long enough for the chaotic behavior to do its work and use that as the initial ensemble of your experiment.

When the uncertainty is dominated by the forcing, e.g. coastal sea-, rivers-, air pollution-, run-off- and sewage-models, you have to work on describing the uncertainty, including time and spatial correlations of these forcings.

When the uncertainty is in the parameters of the model, e.g.  groundwater and run-off models, (and we are not planning to estimate them), you can carefully generate an ensemble of these parameters that represents their uncertainty. Then you set up your ensemble in such a way that each member has a different set of parameters.  Be aware that this setup is not suited for all flavors of EnKF, since the model state after the update must in some sense correspond to the perturbed set of model parameters!

Combinations of the above are possible as well. It is a good investment of time to generate and explore your (initial) ensemble. Note that the filter can only improve your model based on the uncertainty (sub-space) of your ensemble. When important sources are not captured by your ensemble, the filter will not be able to perform well. 

Finally, your model may have time-dependent systematic errors. We often found it useful to add an artificial forcing to the model to describe these model errors.

We will explain in Section \ref{Sec:SequentialEnsembleSimulation} how these experiments can be carried out using OpenDA.


\section{Twin experiments} \label{Sec:Twin}
In real life application we use data assimilation to estimate the true state of the system. Unfortunately we do not know the true state and that makes it difficult to test your data assimilation system. You can set up a so called twin experiment to overcome this problem and test your system in a controlled way. The observations in a twin experiment are generated by a model run with known internal perturbed state or added noise. The perturbation should correspond to the specified uncertainty of your ensemble. Note: Do not use the mean (or deterministic run), because that realization is special. The true state is known in the twin experiment and has the dynamics of your model. This makes it easy to investigate the performance of your data assimilation framework. The SequentialSimulation algorithm in OpenDA is a useful tool for creating your twin experiment.


\section{Workflow}
OpenDA implements a number of algorithms that can be used to gradually grow from a simulation model to a data assimilation system.
\subsection{org.openda.algorithms.Simulation}
Running this algorithm is equivalent to running the model standalone. The only difference it that is that it runs from within OpenDA. It allows you to test whether the configuration is handled correctly and the output of the model can be processed by OpenDA. 

\subsection{SequentialSimulation}
The SequentialSimulation algorithm \footnote{org.openda.algorithms.kalmanFilter.SequentialSimulation}) is again equivalent to running the model by itself. However this time the model is stopped at each moment in which we have observations (or at predefined intervals).  The interpolated model state to the observations are written to the output. 

This algorithm is used to check whether the restart functionality of the model within the OpenDA framework is working correctly (by comparing the results to a normal simulation). Another usage for this algorithm is to create synthetic observations for a twin experiment. You setup observations with arbitrary values but with the location and time you are interested in. After you have run the SequentialSimulation you can find the model predictions that you can use for your synthetic observations. Note: Do not forget to perturb your observation according to the measurement error and perturb the initial state and/or have the model generate noise on the forcing.

\subsection{SequentialEnsembleSimulation}\label{Sec:SequentialEnsembleSimulation}
The SequentialEnsembleSimulation algorithm \footnote{class name org.openda.algorithms.kalmanFilter.SequentialEnsembleSimulation} will propagate your model ensemble without any data assimilation. This algorithm helps you study the behavior of your ensemble. How is explicit noise propagated in to the model? Or how is the initial ensemble propagated? At the same time it is interesting to study the difference between the mean ensemble and your model run. Due to nonlinearities, your mean ensemble can behave significantly differently from your deterministic run.

\subsection{EnKF}
The title suggest to use EnKF  \footnote{class name org.openda.algorithms.kalmanFilter.EnkF}, but other algorithms, e.g. DEnKF or EnSR, are possible as well. However this is the time to start filtering. Start with a twin experiment so that you know that there are no artifacts in the observation data. Start small! First assimilate a small number of observations and take those of which you think that they have a lot of impact. Then start adding observations and see what happens. When you want to assimilate observations from various quantity or quality, first investigate their impact as group and only mix observations in the final steps.

\subsection{Localization, Kalman smoothing, parallel computing, steady state Kalman etc}
To improve performance you can add additional techniques like localization to cope with spurious correlations and steady state filtering or parallel computing filtering to computational performance. OpenDA can output many of the variables involved, such as the Kalman gain. Please, consult the section on ResultWriters for more info.


