% ================================================================
% Chapter stochastic parametrization of EOS (STO)
% ================================================================
\chapter{Stochastic parametrization of EOS (STO)}
\label{STO}

Authors: P.-A. Bouttier

\minitoc

\newpage
$\ $\newline    % force a new line

The stochastic parametrization module aims to explicitly simulate uncertainties in the model. More particularly, \cite{Brankart_OM2013} has shown that, because of the nonlinearity of the seawater equation of state, unresolved scales represent a major source of uncertainties in the computation of the large scale horizontal density gradient (from T/S large scale fields), and that the impact of these uncertainties can be simulated by random processes representing unresolved T/S fluctuations.

The stochastic formulation of the equation of state can be written as:
\begin{equation}
 \label{eq:eos_sto}
  \rho = \frac{1}{2} \sum_{i=1}^m\{ \rho[T+\Delta T_i,S+\Delta S_i,p_o(z)] + \rho[T-\Delta T_i,S-\Delta S_i,p_o(z)] \}
\end{equation}
where $p_o(z)$ is the reference pressure depending on the depth and $\Delta T_i$ and $\Delta S_i$ are a set of T/S perturbations defined as the scalar product of the respective local T/S gradients with random walks $\mathbf{\xi}$:
\begin{equation}
 \label{eq:sto_pert}
 \Delta T_i = \mathbf{\xi}_i \cdot \nabla T \qquad \hbox{and} \qquad \Delta S_i = \mathbf{\xi}_i \cdot \nabla S
\end{equation}
$\mathbf{\xi}_i$ are produced by a first-order autoregressive processes (AR-1) with a parametrized decorrelation time scale, and horizontal and vertical standard deviations $\sigma_s$. $\mathbf{\xi}$ are uncorrelated over the horizontal and fully correlated along the vertical.


\section{Stochastic processes}
\label{STO_the_details}

The starting point of our implementation of stochastic parameterizations
in NEMO is to observe that many existing parameterizations are based
on autoregressive processes, which are used as a basic source of randomness
to transform a deterministic model into a probabilistic model.
A generic approach is thus to add one single new module in NEMO,
generating processes with appropriate statistics
to simulate each kind of uncertainty in the model
(see \cite{Brankart_al_GMD2015} for more details).

In practice, at every model grid point, independent Gaussian autoregressive
processes~$\xi^{(i)},\,i=1,\ldots,m$ are first generated
using the same basic equation:

\begin{equation}
\label{eq:autoreg}
\xi^{(i)}_{k+1} = a^{(i)} \xi^{(i)}_k + b^{(i)} w^{(i)} + c^{(i)}
\end{equation}

\noindent
where $k$ is the index of the model timestep; and
$a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are parameters defining
the mean ($\mu^{(i)}$) standard deviation ($\sigma^{(i)}$)
and correlation timescale ($\tau^{(i)}$) of each process:

\begin{itemize}
\item for order~1 processes, $w^{(i)}$ is a Gaussian white noise,
with zero mean and standard deviation equal to~1, and the parameters
$a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are given by:

\begin{equation}
\label{eq:ord1}
\left\{
\begin{array}{l}
a^{(i)} = \varphi \\
b^{(i)} = \sigma^{(i)} \sqrt{ 1 - \varphi^2 } 
 \qquad\qquad\mbox{with}\qquad\qquad
\varphi = \exp \left( - 1 / \tau^{(i)} \right) \\
c^{(i)} = \mu^{(i)} \left( 1 - \varphi \right) \\
\end{array}
\right.
\end{equation}

\item for order~$n>1$ processes, $w^{(i)}$ is an order~$n-1$ autoregressive process,
with zero mean, standard deviation equal to~$\sigma^{(i)}$; correlation timescale
equal to~$\tau^{(i)}$; and the parameters
$a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are given by:

\begin{equation}
\label{eq:ord2}
\left\{
\begin{array}{l}
a^{(i)} = \varphi \\
b^{(i)} = \frac{n-1}{2(4n-3)} \sqrt{ 1 - \varphi^2 } 
 \qquad\qquad\mbox{with}\qquad\qquad
\varphi = \exp \left( - 1 / \tau^{(i)} \right) \\
c^{(i)} = \mu^{(i)} \left( 1 - \varphi \right) \\
\end{array}
\right.
\end{equation}

\end{itemize}

\noindent
In this way, higher order processes can be easily generated recursively using the same piece of code implementing Eq.~(\ref{eq:autoreg}), and using succesively processes from order $0$ to~$n-1$ as~$w^{(i)}$.
The parameters in Eq.~(\ref{eq:ord2}) are computed so that this recursive application
of Eq.~(\ref{eq:autoreg}) leads to processes with the required standard deviation
and correlation timescale, with the additional condition that
the $n-1$ first derivatives of the autocorrelation function
are equal to zero at~$t=0$, so that the resulting processes
become smoother and smoother as $n$ is increased.

Overall, this method provides quite a simple and generic way of generating a wide class of stochastic processes. However, this also means that new model parameters are needed to specify each of these stochastic processes. As in any parameterization of lacking physics, a very important issues then to tune these new parameters using either first principles, model simulations, or real-world observations.

\section{Implementation details}
\label{STO_thech_details}
The computer code implementing stochastic parametrisations is made of one single FORTRAN module,
with 3 public routines to be called by the model (in our case, NEMO):

The first routine ({sto\_par}) is a direct implementation of Eq.~(\ref{eq:autoreg}),
applied at each model grid point (in 2D or 3D),
and called at each model time step ($k$) to update
every autoregressive process ($i=1,\ldots,m$).
This routine also includes a filtering operator, applied to $w^{(i)}$,
to introduce a spatial correlation between the stochastic processes.

The second routine ({sto\_par\_init})
is an initialization routine mainly dedicated
to the computation of parameters $a^{(i)}, b^{(i)}, c^{(i)}$
for each autoregressive process, as a function of the statistical properties
required by the model user (mean, standard deviation, time correlation,
order of the process,\ldots). Parameters for the processes can be specified through the following namelist parameters:
\begin{alltt}
\tiny
\begin{verbatim}
   nn_sto_eos = 1                ! number of independent random walks 
   rn_eos_stdxy = 1.4            ! random walk horz. standard deviation (in grid points)
   rn_eos_stdz  = 0.7            ! random walk vert. standard deviation (in grid points)
   rn_eos_tcor  = 1440.0         ! random walk time correlation (in timesteps)
   nn_eos_ord  = 1               ! order of autoregressive processes
   nn_eos_flt  = 0               ! passes of Laplacian filter
   rn_eos_lim  = 2.0             ! limitation factor (default = 3.0)
\end{verbatim}
\end{alltt}
This routine also includes the initialization (seeding)
of the random number generator.

The third routine ({sto\_rst\_write}) writes a ``restart file''
with the current value of all autoregressive processes
to allow restarting a simulation from where it has been interrupted.
This file also contains the current state of the random number generator.
In case of a restart, this file is then read by the initialization routine
({sto\_par\_init}), so that the simulation can continue exactly
as if it was not interrupted.
Restart capabilities of the module are driven by the following namelist parameters:
\begin{alltt}
\tiny
\begin{verbatim}
   ln_rststo = .false.           ! start from mean parameter (F) or from restart file (T)
   ln_rstseed = .true.           ! read seed of RNG from restart file
   cn_storst_in  = "restart_sto" !  suffix of stochastic parameter restart file (input)
   cn_storst_out = "restart_sto" !  suffix of stochastic parameter restart file (output)
\end{verbatim}
\end{alltt}

In the particular case of the stochastic equation of state, there is also an additional module ({sto\_pts}) implementing Eq~\ref{eq:sto_pert} and specific piece of code in the equation of state implementing Eq~\ref{eq:eos_sto}.