<figure><img src="image_1.png" /></figure>

HAL
open science

# Introducing time parallelisation within data assimilation using parareal method

Rishabh Bhatt, Laurent Debreu, Arthur Vidard

► To cite this version:

Rishabh Bhatt, Laurent Debreu, Arthur Vidard. Introducing time parallelisation within data assimilation using parareal method. 2022. hal-03540480v2

HAL Id: hal-03540480

https://inria.hal.science/hal-03540480v2

Preprint submitted on 6 Apr 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Introducing time parallelisation within data assimilation

Rishabh Bhatt, Laurent Debreu, Arthur Vidard

► To cite this version:

Rishabh Bhatt, Laurent Debreu, Arthur Vidard. Introducing time parallelisation within data assimilation. SIAM Journal on Scientific Computing, 2025, 47 (2), pp.B533-B557. 10.1137/24M1651903 hal-03540480v3

HAL Id: hal-03540480

https://inria.hal.science/hal-03540480v3

Submitted on 22 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# INTRODUCING TIME PARALLELISATION WITHIN 4D-VAR USING PARAREAL METHOD*

RISHABH BHATT†, LAURENT DEBREU†, AND ARTHUR VIDARD†

**Abstract.** Four dimensional variational data assimilation (4DVAR) methods, in their incremental formulation, are based on optimisation algorithms which require the integration of the forward and adjoint versions of the original model in order to compute its gradient. For their use on parallel computers, these models are classically parallelised only on the spatial dimension and this is a limiting factor on the maximum number of computing cores that can be used. We here present a novel approach of introducing additional time parallelisation using the parareal method. This approach is used here for the integration of the forward model. We use a modified version of the inexact conjugate gradient method where the matrix-vector multiplication is supplied through the parareal algorithm. The use of this inexact conjugate gradient and the associated convergence conditions allows to determine precisely the stopping criterion of the iterations of the parareal method. The results are demonstrated by considering a one dimensional shallow water model. They are presented in terms of the accuracy (in comparison with the original exact conjugate gradient) and in terms of the number of required iterations of the parareal algorithm.

**Key words.** Variational data assimilation, Parareal, Parallel-in-time methods, Optimisation, Shallow water equations

**AMS subject classifications.** 86A22, 68W10, 65L05, 65K10

**1. Introduction.** Predicting weather forecasts accurately for long periods remains a challenge for meteorological institutions. The observations available in the form of in-situ and satellite data are often error-prone and scarcely scattered. While the numerical models expressing the state of the atmosphere comes with assumptions and are discretisations of partial differential equations. The data assimilation methods are run with the objective that an optimum initial state could be recovered which corrects the previous forecasts and gives an improved prediction of the atmosphere. Operational variational data assimilation algorithms [16] are highly sequential in the sense that one needs to perform multiple forward and a backward model integrations to do gradient based minimisation. With billions of observations available along with an extremely high dimensional state variable, this process requires immense computing power.

This is one of the reasons why modern computation emphasises on reducing the clock time. With the advancement of massively fast computers millions of processors can be utilised at once. Speeding up data assimilation with space parallelisation by the means of domain-decomposition is a natural option [29]. Rantakokko [23] provided strategies to calculate the cost function and gradient in parallel by distributing the data (observations, grid point variables) using partitioning algorithms.

Space parallelisation however has its own limits and for large time-dependent problems it alone is insufficient. After one point it gets saturated even before effectively using up all the processors. Again for a climate simulation problem with $10^5$ unknowns per model levels, using $10 \times 10$ unknowns for the horizontal domain would end up using $10^3$ cores. Assigning lesser number of unknowns to each processor would make the computational time so less that the communication time would dominate the overall computing time. This means that after using certain number of cores, there will be a performance loss if more cores are utilised. This would make no sense in weather simulation if it gives us the result on the day for which we actually wanted to predict the weather.

One of the options is to look for the much needed parallelisation in the time domain. The time domain is notorious for its property of causality which makes parallelisation a difficult task. It simply means that the solution at a given time can only be computed if we have some prior information about the solution at previous times. This dependency makes all the algorithms serial or sequential and this is where the concept of the Parallel-in-Time algorithms comes to minimise the influence of causality. It is not surprising that the first ideas of time parallelisation were already introduced in the early 60s by Nivergelt [21] who proposed a multiple shooting method where the solutions at the subintervals are joined together by some interpolation. Much work has been done since then and the state of the art review by Gander [9] provides a comprehensive knowledge of most of the important parallel-in-time algorithms. For a more

<footnote>**Funding:** This work was supported by ANR-ADOM (ANR-18-CE46-0008).</footnote>

<footnote>†Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France (rishabh.bhatt@inria.fr, laurent.debreu@inria.fr, arthur.vidard@inria.fr).</footnote>
# INTRODUCING TIME PARALLELISATION WITHIN DATA ASSIMILATION*

RISHABH BHATT†, LAURENT DEBREU†, AND ARTHUR VIDARD†

**Abstract.** Four dimensional variational data assimilation (4DVAR), in its incremental formulation, is based on optimisation algorithms which require the integration of the forward and adjoint versions of the original model in order to compute the gradient. For their use on parallel computers, these models are classically parallelised only in spatial dimension and this is a limiting factor on the maximum number of cores that can be utilised. We present here a novel approach of introducing additional time parallelisation using the Parareal algorithm. This approach is used here for integration of the forward model. We use a modified version of the inexact conjugate gradient method where the matrix-vector multiplication is supplied through Parareal. The use of this inexact conjugate gradient and the associated convergence conditions allows to precisely determine the stopping criterion of the Parareal iterations. The results are demonstrated by considering a one dimensional shallow water model. They are presented in terms of the accuracy (in comparison with the original exact conjugate gradient) and in terms of the number of required iterations of the Parareal algorithm.

**Key words.** Data assimilation, High performance computing, Optimisation, Numerical analysis, Parareal algorithm

**MSC codes.** 86A22, 68W10, 65L05, 65K10

# 1 Introduction
Modern computation emphasises reducing the clock time. With the advancement of massively parallel computers, millions of cores can be utilised at once. However, for solving very large time-dependent problems space parallelisation alone becomes insufficient. This is particularly relevant for applications such as weather prediction where even small gains in computational speed are significant.

Leading meteorological institutes like Météo France, ECMWF, Met Office use variational data assimilation algorithms like 3D or 4D-Var [27, 5, 39, 40, 35] to obtain optimal initial states for short to medium-range weather forecasts. Data assimilation [8, 24, 1] is essentially an initial condition control problem which combines all the information from the model trajectory, the observations, and some error statistics to give an optimal initial state of the atmosphere/ocean. The evolving nature of the numerical models and the observations in terms of quality and size poses a challenge to the data assimilation systems to meet the demands for faster and accurate predictions. In ocean data assimilation, the state vector and the observation vector are usually of order $O(10^6 - 10^9)$ and $O(10^6)$ respectively. This significantly increases the overall computational cost and a way to deal with it is to improve the use of massively parallel supercomputers.

For instance, ocean models can typically utilise minimum domain sizes of 20 × 20 grid points, beyond which the simulation time increases due to the communication time becoming predominant over the computation time within each subdomain. A global simulation of the NEMO model [29] using ORCA-R025 grid, which has a horizontal domain of 1442 × 1021 grid points and 75 vertical levels, would end up using approximately a maximum of 3672 cores. This means that even if the simulation can have access to more cores, after some point the spatial parallelisation would reach saturation and would not allow to effectively use more cores.

<footnote>*</footnote>

<footnote>**Funding:** This work was supported by ANR-ADOM (ANR-18-CE46-0008).</footnote>

<footnote>†Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France (r.bhatt@reading.ac.uk, laurent.debreu@inria.fr, arthur.vidard@inria.fr).</footnote>
applications-based review article talking about the efficiency and speedup see [22].

There have been a few attempts to introduce time parallelism in data assimilation. Rao and Sandhu [24] presented the strong-constrained 4D-Var as a problem of minimising the Lagrangian with the model trajectories at the sub-intervals acting as constraints. This way the cost function and the gradient evaluations can be performed in parallel. Fischer and Gürol developed a parallel weak-constraint 4D-Var algorithm which works on finding the saddle-point of the Lagrangian [7]. To the authors' knowledge no work has been done to couple proper parallel-in-time method with 4D-Variational data assimilation. Our aim is to introduce time parallelisation with one such method (Parareal) in 4D-Var. We then see if we can exploit some of the available computational power to produce any possible speedup. We present our results for a 1-D linear shallow water model and show a theoretical measure of the speedup when an inexact version of the conjugate gradient method is used for the minimisation.

The paper is structured as follows. In section 2 we briefly talk about Variational data assimilation and how its incremental approach is related to our implementation for a linear model. We will talk about how the time parallelisation could be introduced in its framework. Next we discuss about a specific kind of parallel-time-method called the parareal algorithm which has been central to most of the recent developments in section 3. After the set-up we will talk about the optimisation techniques for minimising the 4D-Var cost function in Section 4. To maintain efficiency we will use the inexact Krylov subspace methods and later provide a modified version of an existing method to make it work with the parareal algorithm. Section 5 is left for the numerical experiments and results. Section 6 concludes the paper with some remarks and talking points for further work.

Throughout this article, we will adopt the following notation

* Upper case calligraphic fonts for nonlinear functions/operators: $\mathcal{F}, \mathcal{J}$

* Upper case italics for linear functions/operators and matrices: $G, H$

* Bold lower case for vectors: $\mathbf{x}, \mathbf{y}$

* Greek letters and lower case for constants and variables

A list of the notation for all the important quantities used in the article is provided in the table below.

<table><thead><tr><th colspan="2">Table of quantities used</th></tr><tr><th>Symbol</th><th>Meaning</th></tr></thead><tbody><tr><td><b>x<sub>i</sub></b></td><td>state variable at time t<sub>i</sub></td></tr><tr><td><b>Δx</b></td><td>state variable increment</td></tr><tr><td><b>Δx</b></td><td>spatial grid stepping</td></tr><tr><td><b>y<sub>i</sub></b></td><td>observation at time t<sub>i</sub></td></tr><tr><td><b>F</b></td><td>fine propagator/model</td></tr><tr><td><b>G</b></td><td>coarse propagator</td></tr><tr><td><b>E</b></td><td>perturbation matrix</td></tr><tr><td><b>r<sub>j</sub></b></td><td>inexact residual at inexact CG iteration j</td></tr><tr><td><b>P(k)</b></td><td>parareal operator at kth iteration</td></tr><tr><td><b>p<sub>j</sub></b></td><td>conjugate direction at inexact CG iteration j</td></tr><tr><td><b>ω<sub>j</sub></b></td><td>tolerance for ||<b>E<sub>j</sub></b>||<sub>A<sup>-1</sup>,A</sub></sub></td></tr><tr><td><b>ξ<sub>j</sub></b></td><td>tolerance for ||<b>E<sub>j</sub>p<sub>j</sub></b>||<sub>A<sup>-1</sup></sub></sub></td></tr><tr><td><b>κ(A)</b></td><td>condition number of a matrix A</td></tr></tbody></table>

2. **Data assimilation.** Let us describe our dynamical model depending on a state variable vector $\mathbf{x} \in \mathbb{R}^n$ by

$$ (2.1) \qquad \begin{aligned} \mathbf{x}(t_0) &= \mathbf{x}_0 \\ \mathbf{x}_i &= \mathcal{F}_{t_{i-1} \to t_i}(\mathbf{x}_{i-1}) \quad i = 1, \dots, N \end{aligned} $$

where $\mathcal{F}$ is the discrete non-linear model operator, $\mathbf{x}_0$ is the initial condition and $N$ is the number of time windows.
In large time-dependent problems such as the one described in the above example, it is tempting to use the time dimension as an additional source of parallelism. This is one of the main reasons behind turning our attention to the parallel-in-time (PinT) methods [34, 14] in addition to spatial domain decomposition methods. This observation applies to direct modelling, but also, and this is what will be of interest in this paper, to data assimilation algorithms.

Previous works have explored a “time-parallel 4D-Var”, with attempts to parallelise weak-constraint 4D-Var in [12] by solving a saddle point problem. A PDE constrained optimisation problem to solve the coupled system of the forward and backward evolution problems with the help of Parareal (a parallel-in-time algorithm) is proposed in [16]. However, a drawback of this approach is that it is not easy to control the initial condition of the model, which is central to data assimilation. Other works have also investigated Parareal in combination with minimisation iterations to solve optimal control problems [32, 9, 28]. Our main objective here is to study the exploitation of time parallelism by coupling Parareal with incremental 4D-Var. We are going to focus specifically on the inner loop of incremental 4D-Var and so we consider the case where our model is linear.

The paper is structured as follows. Section 2 talks about the basics of variational data assimilation and Parareal algorithm. Section 3 delves into the coupling of Parareal with incremental 4D-Var, focusing on a minimisation problem involving inexact matrix-vector products. The development of a modified inexact conjugate gradient method is outlined in Section 4. Numerical experiments based on a 1D shallow water model are presented in Section 5. The conclusion follows in Section 6, offering some perspectives on our work to envision more realistic applications. Throughout this article, we will adopt the notations from Ide et al.[22].

## 2 Background

**2.1 Data Assimilation** Four-dimensional variational data assimilation, or 4D-Var [39, 40, 35] minimises a cost function based on discrepancies between the model trajectory and the observations, and between the obtained analysis and the previous known background state. For the model state vector $\mathbf{x}$ we consider the following discrete non-linear dynamical model

$$ (2.1) \qquad \begin{aligned} \mathbf{x}_0 &= \mathbf{x}(t_0) \\ \mathbf{x}_i &= F_i(\mathbf{x}_{i-1}) \quad i = 1, 2, \dots, N \end{aligned} $$

where $F_i : \mathbb{R}^n \to \mathbb{R}^n$ is the model operator describing the state evolution from time $t_{i-1}$ to $t_i$, and $N$ is the total number of time windows. Given a prior estimate/background state $\mathbf{x}^b \in \mathbb{R}^n$ and a set of observations $\mathbf{y}_i^o \in \mathbb{R}^m$ at time $t_i$, the 4D-Var cost function is written as

$$ (2.2) \quad J(\mathbf{x}_0) = \frac{1}{2} (\mathbf{x}_0 - \mathbf{x}^b)^T \mathbf{B}^{-1} (\mathbf{x}_0 - \mathbf{x}^b) + \frac{1}{2} \sum_{i=0}^{N} (H_i(\mathbf{x}_i) - \mathbf{y}_i^o)^T \mathbf{R}_i^{-1} (H_i(\mathbf{x}_i) - \mathbf{y}_i^o) $$

Here $H_i : \mathbb{R}^n \to \mathbb{R}^m$ is the observation operator mapping the evaluated model state $\mathbf{x}_i$ to the observation $\mathbf{y}_i^o$ at the correct time $t_i$. The errors in the background and the observations are characterised by the corresponding covariance matrices $\mathbf{B}$ and $\mathbf{R}_i$.

**2.2 Incremental 4D-Var formulation** Incremental 4D-Var introduced in [6] is an algorithm obtained by linearising the model and observation operators around
The problem of 4D-Variational data assimilation or 4D-Var amounts to correcting the model trajectory by striking a right balance between the a priori background information $\mathbf{x}_0^b \in \mathbb{R}^n$ and a set of observations $\mathbf{y}_i^o \in \mathbb{R}^m$ at time $t_i$ [3]. This can be expressed as a least-squares fit problem whose cost function is

$$ \mathcal{J}(\mathbf{x}_0) = \frac{1}{2}(\mathbf{x}_0 - \mathbf{x}_0^b)^T B^{-1}(\mathbf{x}_0 - \mathbf{x}_0^b) + \frac{1}{2} \sum_{i=0}^{N} (\mathbf{y}_i^o - \mathcal{H}_i(\mathbf{x}_i))^T R_i^{-1} (\mathbf{y}_i^o - \mathcal{H}_i(\mathbf{x}_i)) $$

Here $B$ and $R$ are the background and observation covariance matrices respectively. $\mathcal{H}: \mathbb{R}^n \to \mathbb{R}^m$ is the observation operator which maps the state vectors to the observation space. Since the entire model trajectory can be described by the initial condition $\mathbf{x}_0$, the cost function can be minimised only with respect to $\mathbf{x}_0$. Thus the minimisation problem can be reformulated as

$$ (2.2) \qquad \begin{aligned} \min_{\mathbf{x}_0 \in \mathbb{R}^n} \mathcal{J}(\mathbf{x}_0) &= \frac{1}{2} \|\mathbf{x}_0 - \mathbf{x}_0^b\|_b^2 + \frac{1}{2} \sum_{i=0}^{N} \|\mathbf{y}_i^o - \mathcal{H}_i(\mathcal{F}_{t_0 \to t_i}(\mathbf{x}_0))\|_o^2 \\ &= \underbrace{\mathcal{J}_b(\mathbf{x}_0)}_{\text{background}} + \underbrace{\mathcal{J}_o(\mathbf{x}_0)}_{\text{observation}} \end{aligned} $$

This is known as the strongly constrained 4D-Var since the numerical model is assumed to be perfect. The norms $\| \cdot \|_b$ and $\| \cdot \|_o$ are the energy norms of the matrices $B^{-1}$ and $R^{-1}$ respectively.

Looking at (2.2) it is evident that the cost function $\mathcal{J}$ is non-convex because of the nonlinearites in $\mathcal{F}$ and $\mathcal{H}$. As a result in the absence of a unique minimum, the minimisation becomes very difficult. Instead of solving the large-scale minimisation process at one go, the problem is solved incrementally by taking successive approximations [4]. This is done by taking the following tangent linear approximations around the iterates with sufficiently small increments $\Delta\mathbf{x}^{(k)}$

$$ \begin{aligned} \mathcal{F}(\mathbf{x}^{(k)} + \Delta\mathbf{x}^{(k)}) &= \mathcal{F}(\mathbf{x}^{(k)}) + F|_{\mathbf{x}^{(k)}} \Delta\mathbf{x}^{(k)} \\ \mathcal{H}(\mathbf{x}^{(k)} + \Delta\mathbf{x}^{(k)}) &= \mathcal{H}(\mathbf{x}^{(k)}) + H|_{\mathcal{F}(\mathbf{x}^{(k)})} \Delta\mathbf{x}^{(k)} \\ \mathcal{H} \circ \mathcal{F}(\mathbf{x}^{(k)} + \Delta\mathbf{x}^{(k)}) &= \mathcal{H} \circ \mathcal{F}(\mathbf{x}^{(k)}) + HF\Delta\mathbf{x}^{(k)} \end{aligned} $$

where $F = \frac{\partial \mathcal{F}}{\partial \mathbf{x}}|_{\mathbf{x}^{(k)}}$ and $H = \frac{\partial \mathcal{H}}{\partial \mathbf{x}}|_{\mathcal{F}(\mathbf{x}^{(k)})}$ are the linearised model and observation operator respectively.

The cost function obtained is now quadratic and much easier to solve. This formulation is known as the Incremental 4D-Var. We present the procedure as provided in [15]

1. For $k = 0$ we usually set $\mathbf{x}_0^{(0)} = \mathbf{x}_0^b$.

2. At iteration $k$ for a given increment defined as $\Delta\mathbf{x}_0^{(k)} = \mathbf{x}_0^{(k+1)} - \mathbf{x}_0^{(k)}$ the cost function (2.2) can be written as

$$ (2.3) \qquad \begin{aligned} J(\Delta\mathbf{x}_0^{(k)}) &= \frac{1}{2} \|\Delta\mathbf{x}_0^{(k)} - (\mathbf{x}_0^b - \mathbf{x}_0^{(k)})\|_b^2 + \frac{1}{2} \sum_{i=0}^{N} \|\mathbf{y}_i - \mathcal{H}_i \circ \mathcal{F}(\mathbf{x}_i^{(k)} + \Delta\mathbf{x}_i^{(k)})\|_o^2 \\ &= \frac{1}{2} \|\Delta\mathbf{x}_0^{(k)} - (\mathbf{x}_0^b - \mathbf{x}_0^{(k)})\|_b^2 + \frac{1}{2} \sum_{i=0}^{N} \|\mathbf{d}_i - H_i F \Delta\mathbf{x}_i^{(k)}\|_o^2 \end{aligned} $$

where we define the innovation vector or the departure as $\mathbf{d}_i^{(k)} = \mathbf{y}_i - \mathcal{H}_i \circ \mathcal{F}(\mathbf{x}_i^{(k)})$

3. The constraint is given by

$$ (2.4) \qquad \begin{aligned} \Delta\mathbf{x}_i^{(k)} &= F_{t_0 \to t_i} \Delta\mathbf{x}_0^{(k)} \\ &= F_i F_{i-1} \dots F_2 F_1 \Delta\mathbf{x}_0^{(k)} \end{aligned} $$

4. Note that the control variable now is the correction $\Delta\mathbf{x}$ to the model parameter $\mathbf{x}$. The new cost function now is quadratic with respect to $\Delta\mathbf{x}_0^{(k)}$. The guess is updated by

$$ \mathbf{x}_0^{(k+1)} = \mathbf{x}_0^{(k)} + \Delta\mathbf{x}_0^{(k)} $$
the model state. For a sequence of linearisation states $\mathbf{x}^{(0)}, \mathbf{x}^{(1)}, \dots, \mathbf{x}^{(l)}$ with corresponding increments $\delta\mathbf{x}^{(0)}, \delta\mathbf{x}^{(1)}, \dots, \delta\mathbf{x}^{(l)}$ the following tangent linear hypothesis (a first order Taylor's expansion) is expressed as

$$ (2.3) \qquad \begin{aligned} F(\mathbf{x}^{(l)} + \delta\mathbf{x}^{(l)}) &\approx F(\mathbf{x}^{(l)}) + \mathbf{F} \delta\mathbf{x}^{(l)} \\ H(\mathbf{x}^{(l)} + \delta\mathbf{x}^{(l)}) &\approx H(\mathbf{x}^{(l)}) + \mathbf{H} \delta\mathbf{x}^{(l)} \end{aligned} $$

where $\mathbf{F}$ and $\mathbf{H}$ are the linearisation of the model and observation operators respectively around the linearisation state $\mathbf{x}^{(l)}$; that is,

$$ (2.4) \qquad \mathbf{F} = \frac{\partial F}{\partial \mathbf{x}}|_{\mathbf{x}^{(l)}} \qquad \mathbf{H} = \frac{\partial H}{\partial \mathbf{x}}|_{F(\mathbf{x}^{(l)})} $$

The control variable now changes from the initial model state $\mathbf{x}_0$ to the initial increment $\delta\mathbf{x}_0$ of the model state. Using the above assumption we get our cost function

$$ (2.5) \qquad J(\delta\mathbf{x}_0) = \frac{1}{2}\{\delta\mathbf{x}_0 - (\mathbf{x}^b - \mathbf{x}_0)\}^T \mathbf{B}^{-1} \{\delta\mathbf{x}_0 - (\mathbf{x}^b - \mathbf{x}_0)\} + \frac{1}{2} \sum_{i=0}^{N} (\mathbf{H}_i \delta\mathbf{x}_i - \mathbf{d}_i)^T \mathbf{R}_i^{-1} (\mathbf{H}_i \delta\mathbf{x}_i - \mathbf{d}_i) $$

where $\mathbf{d}_i = \mathbf{H}_i(\mathbf{x}_i) - \mathbf{y}_i$ is called the innovation vector or just departures. Clearly one has to run the non-linear model integration and store the trajectory to compute the departures since the non-linear model $M$ is embedded within $\delta\mathbf{x}_i$. The increment vector evolves in time through the *tangent linear model*, obtained by linearising the discrete non-linear model

$$ (2.6) \qquad \delta\mathbf{x}_{i+1} = \frac{\partial F[\mathbf{x}(t_i)]}{\partial \mathbf{x}} \delta\mathbf{x}_i = \mathbf{F}_i \delta\mathbf{x}_i $$

In a nutshell, the incremental 4D-Var involves running a high resolution model in the outer loop, which retains all the physics, to compare the model state trajectory with observations for cost function evaluation. The inner loop involves running a low resolution model with simplified physics, to minimise the cost function [35]. For an algorithmic implementation of incremental 4D-Var we refer to [25].

**2.3 Parareal Method** Parareal initially introduced in [26] and later reformulated [2] is a parallel-in-time method where the solution at a given time step is computed independently of the solution at previous time steps. For a general class of ODE

$$ (2.7) \qquad \mathbf{x}'(t) = f(t, \mathbf{x}(t)), \quad t \in (0, T], \quad \mathbf{x}(0) = \mathbf{x}_0 $$

the strategy relies on partitioning the time domain $[0, T]$ into $N$ time windows or sub-intervals $[t_i, t_{i+1}]$ of length $\Delta T = T/N$. Within each such sub-interval, a two-level grid system is defined consisting of:

* a cheap coarse propagator $G(t_{i+1}, t_i, \mathbf{x}_i)$ which is fast but gives a rough solution at $t_{i+1}$ from the initial condition $\mathbf{x}_i$.

* a fine propagator $F(t_{i+1}, t_i, \mathbf{x}_i)$ which is computationally expensive but gives very accurate solution (in theory the exact solution).

The Parareal algorithm proceeds iteratively by computing the solution $\mathbf{x}_{i+1}^k \approx \mathbf{x}(t_{i+1})$ at the sub-intervals for given $\mathbf{x}_i^k$ at time $t_i$ and iteration $k$. An initial configuration is
**Note:** It is more efficient for solving quadratic problems but this does not count out the possibility that it can get stuck in a local minimum.

**2.1. Introducing time parallelisation.** Each minimisation iteration of a 4D-Var (or incremental 4D-Var inner) cycle begins by a forward integration of the model (2.1) for calculating the cost function value. The obtained trajectory is stored for running the adjoint model for gradient computation. A potential speedup can be achieved if these trajectories are instead calculated by integrating the model over the time windows parallel in time. To see how this could be done, let us suppose the length of one cycle to be $T$ hours (say) which is made up of $N$ time windows.

<figure><img src="image_2.png" /><figcaption>FIG. 2.1. Using time parallelisation for running the forward model</figcaption></figure>

Let place ourselves in the case of a linear forward model

$$ (2.5) \qquad \begin{aligned} \mathbf{x}(t_0) &= \mathbf{x}_0 \\ \mathbf{x}_{i+1} &= F\mathbf{x}_i \quad i = 1, \dots, N \end{aligned} $$

The model operator $F$ is the same as the one used in (2.4). The only difference is that we have now replaced the state variable $\Delta\mathbf{x}$ with $\mathbf{x}$ to simplify notations.

Let us now set up our data assimilation problem. Again for the sake of simplicity we assume that we have no background information about the model, that is there is no background term in the cost function. Also let there be only one true observation $\mathbf{y}$ at the end of the integration time $t_N = T$. This means that the covariance matrix $R_i$ is equal to the identity matrix while the observation operator $\mathcal{H}_i$ is an identity map for each time $t_i$. Thus our cost function can be written as

$$ (2.6) \qquad J(\mathbf{x}_0) = \frac{1}{2} \|F^N \mathbf{x}_0 - \mathbf{y}\|_2^2 $$

with gradient,

$$ \nabla J(\mathbf{x}_0) = (F^N)^T (F^N \mathbf{x}_0 - \mathbf{y}) $$

where $(F^N)^T$ is the adjoint matrix.

As stated before, if we now carry out the computation of the forward model (involved in the misfit between the model trajectory and the observation) by a parallel-in-time based operator $P$ (can be iterative or
<figure><img src="image_2.png" /><figcaption>FIG. 1. Illustration of the Parareal algorithm for a sub-interval $[t_i, t_{i+1}]$. For a Parareal iterate $\mathbf{x}_i^k$ at iteration $k$, $G(t_{i+1}, t_i, \mathbf{x}_i^k)$ gives solution by coarse propagator $G$ using time-step $\Delta t$. $F(t_{i+1}, t_i, \mathbf{x}_i^k)$ gives solution by fine propagator $F$ with time step $\delta t$.</figcaption></figure>

obtained by a serial run of the coarse propagator $G$ over the whole time domain. So, the initial Parareal iteration $k = 0$ looks like

$$ (2.8) \qquad \mathbf{x}_{i+1}^0 = G(t_{i+1}, t_i, \mathbf{x}_i^0) $$

With the solution now available at all the grid points, the fine propagator $F$ is run to provide more accurate solution. As $F$ is expensive, all the heavy fine computations performed by $F$ are done in parallel by assigning a core to each sub-interval. The iterates are updated by a correction procedure

$$ (2.9) \qquad \begin{aligned} \mathbf{x}_0^{k+1} &= \mathbf{x}_0 \\ \mathbf{x}_{i+1}^{k+1} &= \underbrace{G(t_{i+1}, t_i, \mathbf{x}_i^{k+1})}_{\text{Prediction}} + \underbrace{F(t_{i+1}, t_i, \mathbf{x}_i^k) - G(t_{i+1}, t_i, \mathbf{x}_i^k)}_{\text{Correction}} \end{aligned} $$

Parareal can be seen as a predictor-corrector algorithm where each iteration predicts with $G$ and corrects with the difference between the solutions from $F$ and $G$ at the previous iteration. An in-depth discussion of the convergence properties of Parareal can be found in [11, 3, 17, 15, 36].

### 2.3.1 Speedup
We briefly discuss about the speedup of Parareal algorithm which we will be used as a criteria to check the performance of our experiments. The theoretical speedup $S$ is defined as the ratio of the time taken to perform a task sequentially to the time to do the same task in parallel,

$$ (2.10) \qquad S = \frac{\text{Serial computation time}}{\text{Parallel computation time}} $$

Assuming $N$ cores are available, each assigned to one time window, and negligible communication time between any two cores, the speedup is bounded by [33, 37]

$$ (2.11) \qquad S \le \min \left( \frac{N}{k}, \frac{\text{fine time}}{\text{coarse time}} \right) $$

The bound (2.11) shows that the coarse propagator should be really cheap and fast in order to get a good speed up but at the same time there will be more Parareal iterations which will reduce the speedup according to the first bound. Thus the optimal speedup must balance the two bounds in the best possible way.
direct), we have

$$ (2.7) \qquad \begin{aligned} J(\mathbf{x}_0) &\approx \frac{1}{2} \|P\mathbf{x}_0 - \mathbf{y}\|_2^2 \\ \nabla J(\mathbf{x}_0) &\approx (F^N)^T (P\mathbf{x}_0 - \mathbf{y}) \end{aligned} $$

Minimising this cost function requires setting the gradient $\nabla J(\mathbf{x}_0)$ to 0. This is equivalent to solving the system

$$ A\mathbf{x} = \mathbf{b} \quad \text{with} \quad A = (F^N)^T P, \quad \mathbf{b} = (F^N)^T \mathbf{y} $$

At the end of each minimisation iteration we expect to have an initial condition with a lower cost function value. This in turn will be used by $P$ to provide new trajectories for the next iteration. The process will end depending on the stopping criterion of the minimisation.

When $\kappa(A)$, the conditioning of $A$, is large, a regularisation term $\alpha$ can be used in the cost function which depends on the spatial derivative of the initial state $\mathbf{x}_0$. Thus,

$$ J(\mathbf{x}_0) \approx \frac{1}{2} \|P\mathbf{x}_0 - \mathbf{y}\|_2^2 + \frac{\alpha}{2} \left\| \frac{\mathrm{d}\mathbf{x}_0}{\mathrm{d}\mathbf{x}} \right\|_2^2 $$

The spatial derivative can be written as

$$ \frac{\mathrm{d}\mathbf{x}_0}{\mathrm{d}\mathbf{x}} = \frac{\mathbf{x}_{0,i+1} - \mathbf{x}_{0,i}}{\Delta x} = \frac{1}{\Delta x} Q_1 \mathbf{x}_0 $$

where

$$ Q_1 = \begin{pmatrix} -1 & 1 & & & \\ & -1 & 1 & & \\ & & \ddots & \ddots & \\ & & & -1 & \end{pmatrix} $$

Therefore,

$$ \left\| \frac{\mathrm{d}\mathbf{x}_0}{\mathrm{d}\mathbf{x}} \right\|_2^2 = \frac{1}{\Delta x^2} \langle Q_1 \mathbf{x}_0, Q_1 \mathbf{x}_0 \rangle $$

The simplified cost function can be written as

$$ J(\mathbf{x}_0) \approx \frac{1}{2} \|P\mathbf{x}_0 - \mathbf{y}\|_2^2 + \frac{\alpha}{2\Delta x^2} \|Q_1 \mathbf{x}_0\|_2^2 $$

Written like this, the gradient should be

$$ \nabla J(\mathbf{x}_0) \approx (F^N)^T (P\mathbf{x}_0 - \mathbf{y}) + \alpha Q_2 \mathbf{x}_0 $$

where

$$ Q_2 = \frac{1}{\Delta x^2} Q_1^T Q_1 = \frac{1}{\Delta x^2} \begin{pmatrix} 2 & -1 & & & \\ -1 & 2 & -1 & & \\ & \ddots & \ddots & \ddots & \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \end{pmatrix} $$

Again, this can be put up as a linear system $A\mathbf{x} = \mathbf{b}$ with

$$ A = (F^N)^T P + \alpha Q_2, \quad \mathbf{b} = (F^N)^T \mathbf{y} $$

**REMARK 2.1.** Using an iterative minimisation algorithm for the above linear system is equivalent to solving the incremental 4D-Var problem.
### 3 Introducing time parallelisation
The inner loop of the incremental 4D-Var remains the least scalable component due to the use of lower resolution models [23]. This limitation arises from having fewer grid points available for distributing across cores when producing the analysis increment. In fact, parallelising the inner loop is by no means trivial. The most obvious choice for minimisation within the inner loop is the conjugate gradient method, one of the Krylov-subspace based methods. It is constructed on the principle of finding higher level Krylov subspaces at each iteration, relying on the initial gradient (or residual) of the cost function. These gradients are updated at the end of the iteration and are used to build the next Krylov subspace [13]. However these iterations are inherently sequential, leaving parallelising the computations within the iteration as the only alternative to parallelising the minimisation iterations themselves.

Our approach is as follows. Each minimisation iteration of an incremental 4D-Var (or inner loop) cycle begins with a forward integration of the linear model (2.6) to calculate the quadratic cost function value. An immediate way to introduce time parallelisation is to use Parareal algorithm for the forward model integration. To see how this could be done, let us suppose the length of one data assimilation cycle to be T which is made up of N time windows. Consider the case of a linear forward model for the state vector $\delta\mathbf{x} \in \mathbb{R}^n$

$$ (3.1) \qquad \begin{aligned} \delta\mathbf{x}_0 &= \delta\mathbf{x}(t_0) \\ \delta\mathbf{x}_{i+1} &= \mathbf{F}_i \delta\mathbf{x}_i, \quad i = 1, \dots, N \end{aligned} $$

The model operator $\mathbf{F}_i : \mathbb{R}^n \to \mathbb{R}^n$ is the same as the one used in (2.6) and so

$$ (3.2) \qquad \mathbf{F}_i = \frac{\partial F}{\partial \mathbf{x}} \Big|_{\mathbf{x}=\mathbf{x}_i} $$

For the sake of simplicity we remove the background term $\mathbf{x}^b$ from the data assimilation cost function since it is not affected by time integration. We also assume that there is only one observation $\mathbf{y} \in \mathbb{R}^n$ of the whole state vector at the end of integration time $t_N = T$. Doing so the observation operator $H$ is an identity map for each time $t_i$ and the minimisation problem is well posed. As a further simplification, the covariance matrix $\mathbf{R}_i$ is set to be equal to the identity matrix. Our cost function can be written as

$$ (3.3) \qquad J(\delta\mathbf{x}_0) = \frac{1}{2} \left\| \left( \prod_{i=1}^{N} \mathbf{F}_i \right) \delta\mathbf{x}_0 - \mathbf{y} \right\|_2^2 $$

*Remark 3.1.* In the rest of this manuscript instead of writing the composition of the discrete model operators $\mathbf{F}_i$ at a model state $\delta\mathbf{x}_i$ as $\prod_{i=1}^N \mathbf{F}_i$ each time we are going to use the block $\mathbf{F}^N$ for $\mathbf{F}_1\mathbf{F}_2\cdots\mathbf{F}_N$ for convenience and it is just an abuse of the notation.

Thus we have,

$$ (3.4) \qquad J(\delta\mathbf{x}_0) = \frac{1}{2} \|\mathbf{F}^N \delta\mathbf{x}_0 - \mathbf{y}\|_2^2 $$

with gradient,

$$ (3.5) \qquad \nabla J(\delta\mathbf{x}_0) = (\mathbf{F}^N)^T (\mathbf{F}^N \delta\mathbf{x}_0 - \mathbf{y}) $$
Depending on the choice of the operator P the two iterative processes of forward model integration and minimisation can be coupled in different ways. We are going to focus on the case where P is replaced by the parareal algorithm, which we will present in the next section.

**3. Parareal method.** Parareal was introduced in [17] and later reformulated in [1] as a kind of parallel-in-time method where the solution at a given time step is computed independently of the solution at previous time steps. For a general discrete model

$$ (3.1) \qquad \begin{aligned} \mathbf{x}(t_0) &= \mathbf{x}_0 \\ \mathbf{x}_{i+1} &= \mathcal{F}\mathbf{x}_i \quad i = 1, \dots, N \end{aligned} $$

the strategy relies on partitioning the time domain $[0, T]$ into $N$ time windows or sub-intervals $[t_n, t_{n+1}]$ of length $\Delta T = T/N$. Within each such sub-interval, a two-level grid system is defined using a coarse propagator $\mathcal{G}$ and a fine propagator $\mathcal{F}$. $\mathcal{G}$ is chosen so that it is a very cheap approximation of $\mathcal{F}$, with large space/time steps making it less accurate. On the other hand, $\mathcal{F}$ has small space/time steps making it a highly accurate and computationally expensive solver. In the examples of this paper, only the time discretisation is different between $\mathcal{F}$ and $\mathcal{G}$, respectively denoted $\delta t$ and $\Delta t$.

<figure><img src="image_3.png" /></figure>

The parareal algorithm proceeds iteratively by computing the solution $\mathbf{x}_{n+1}^k \approx \mathbf{x}(t_{n+1})$ at the sub-intervals for given $\mathbf{x}_n^k$ at time $t_n$ and iteration $k$. An initial configuration ($k=0$) is obtained by a serial run of the coarse propagator $\mathcal{G}$ over the whole time domain. Thus,

$$ \begin{aligned} \mathbf{x}_0^0 &= \mathbf{x}_0 \\ \mathbf{x}_{n+1}^0 &= \mathcal{G}(t_{n+1}, t_n, \mathbf{x}_n^0) \quad n = 0, 1, \dots, N-1 \end{aligned} $$

With the solution values now available at all the grid points, the fine propagator $\mathcal{F}$ can be run to provide more accurate estimates. As $\mathcal{F}$ is expensive, all the computations performed by $\mathcal{F}$ are done in parallel by assigning a processor to each sub-interval. Finally the iterates are updated by a correction procedure

$$ (3.2) \qquad \begin{aligned} \mathbf{x}_0^{k+1} &= \mathbf{x}_0 \\ \mathbf{x}_{n+1}^{k+1} &= \underbrace{\mathcal{G}(t_{n+1}, t_n, \mathbf{x}_n^{k+1})}_{\text{Prediction}} + \underbrace{\mathcal{F}(t_{n+1}, t_n, \mathbf{x}_n^k) - \mathcal{G}(t_{n+1}, t_n, \mathbf{x}_n^k)}_{\text{Correction}} \end{aligned} $$

Thus parareal can be seen as a predictor-corrector algorithm where at each iteration the predicted term $\mathcal{G}(t_{n+1}, t_n, \mathbf{x}_n^{k+1})$ is corrected by the solutions from $\mathcal{F}$ and $\mathcal{G}$ at the previous iteration. The iterate update step (3.2) is sequential as one has to do a serial coarse integration first to get the prediction term and only then the pre-computed correction terms can be added.

Theoretically, the maximum accuracy that the parareal algorithm can achieve is that of the fine solver. Taking the limit $k \to \infty$, the solutions from the coarse solver in (3.3) cancel out each other and what is left is just the solution from $\mathcal{F}$. Another thing to note is that after each parareal iteration (p-iteration) the trajectory is corrected to its exact solution for exactly one time window. After $N$ iterations, the parareal algorithm will produce exactly the same solution as the solution obtained from the serial integration by the fine solver. For the algorithm to be practical it is important that $k \ll N$, otherwise there will be no speedup.
where $(\mathbf{F}^N)^T$ is the adjoint of the operator $\mathbf{F}^N$. To carry out computation of the forward model by a Parareal based operator, we accordingly define a Parareal based integrator $\mathbf{P}(k)$ over $[0, T]$ such that

$$ (3.6) \qquad \mathbf{P}(k) \delta\mathbf{x}_0 = \delta\mathbf{x}_N^k $$

which acts on an initial condition $\delta\mathbf{x}_0$ and gives Parareal solution at the last time window $N$ after $k$ iterations. The gradient of the cost function is then approximated by

$$ (3.7) \qquad \begin{aligned} \nabla J(\delta\mathbf{x}_0) &\approx (\mathbf{F}^N)^T (\mathbf{P}(k) \delta\mathbf{x}_0 - \mathbf{y}) \\ &\stackrel{\text{def}}{=} \overline{\nabla J(\delta\mathbf{x}_0)} \end{aligned} $$

Note that we kept the same adjoint for the gradient i.e. $(\mathbf{F}^N)^T$ and not the adjoint of $\mathbf{P}(k)$. Now minimising the original cost function in (3.4) requires finding $\delta\mathbf{x}_0$ such that gradient $\nabla J(\delta\mathbf{x}_0)$ is 0. This equates to solving the linear system

$$ (3.8) \qquad \mathbf{A}\delta\mathbf{x} = \mathbf{b} $$

where $\mathbf{A} = (\mathbf{F}^N)^T \mathbf{F}^N$, $\mathbf{b} = (\mathbf{F}^N)^T \mathbf{y}$. Similarly, the corresponding approximate linear system writes

$$ (3.9) \qquad \mathbf{A}_{\text{para}} \delta \mathbf{x} = \mathbf{b} $$

where $\mathbf{A}_{\text{para}} = (\mathbf{F}^N)^T \mathbf{P}(k)$.

Remark 3.2. In realistic applications, matrix $\mathbf{H}$ is typically low-rank rather than the identity, meaning we observe only a subset of the full state variable. This affects the conditioning of the minimisation which can be improved by adding more observations at different times. Additionally, we need to introduce observation-error covariance matrix $\mathbf{R}$ to account for uncertainties. If observations $\mathbf{y}_i$ are available at times $t_i$, our approach can be extended to this general case by solving the linear system $\mathbf{A}\mathbf{x} = \mathbf{b}$ where

$$ (3.10) \qquad \mathbf{A} = \frac{1}{p} \sum_{i=1}^{p} (\mathbf{H}_i \mathbf{F}^{N_i})^T \mathbf{R}_i^{-1} \mathbf{H}_i \mathbf{P}^i(k), \qquad \mathbf{b} = \frac{1}{p} \sum_{i=1}^{p} (\mathbf{H}_i \mathbf{F}^{N_i})^T \mathbf{R}_i^{-1} \mathbf{y}_i $$

and $\mathbf{F}^{N_i}, \mathbf{P}^i(k), \mathbf{R}_i$, and $\mathbf{H}_i$ are evaluated at time $t_i$.

The exact linear system $\mathbf{A}\mathbf{x} = \mathbf{b}$ is usually solved using the conjugate gradient (CG) method since the matrix $\mathbf{A}$ is symmetric. But as the corresponding approximate matrix $\mathbf{A}_{\text{para}} = (\mathbf{F}^N)^T \mathbf{P}(k)$ is non-symmetric, there is a big possibility that CG will not give the desired solution and thus the correct gradient value of our cost function. We address this challenge by using an *inexact Krylov subspace method* [19, 4] which is based on a counter-intuitive principle that the error in the matrix-vector multiplication is allowed to grow as iterations progress [38, 42]. By doing so, same accuracy levels for the minimisation process can be maintained as expected with the usual (exact) Krylov subspace methods. In the next section we discuss about the *inexact conjugate gradient method*, a kind of *inexact Krylov subspace method*.

**4 Inexact conjugate gradient method** In this section we talk about the inexact conjugate gradient method proposed by Gratton et al. [20] and then introduce our own modified version of this method with the Parareal operator $\mathbf{P}(k)$ embedded in it.
<caption>**Algorithm 3.1 Parareal Algorithm**</caption>

1: **Initialisation**
2: **for** n ← 0 to N **do**
3: $\mathbf{x}_{n+1}^0 = \mathcal{G}(t_{n+1}, t_n, \mathbf{x}_n^0)$
4: **end for**
5: **{Iterations}**
6: $k = 0$
7: **repeat**
8: **{parallel prediction step}**
9: **for** n ← 0 to N **do**
10: $\tilde{\mathbf{x}}_{n+1}^k = \mathcal{F}(t_{n+1}, t_n, \mathbf{x}_n^k)$
11: **end for**
12: **{Serial correction step}**
13: **for** n ← 0 to N **do**
14: $\mathbf{x}_{n+1}^{k+1} = \mathcal{G}(t_{n+1}, t_n, \mathbf{x}_n^{k+1}) + \tilde{\mathbf{x}}_{n+1}^k - \mathcal{G}(t_{n+1}, t_n, \mathbf{x}_n^k)$
15: **end for**
16: $k = k + 1$
17: **until** $k = k_{\max}$

If the model (3.1) is linear, we can re-write the correction step (3.2) as

$$ (3.3) \qquad \begin{aligned} \mathbf{x}_{n+1}^{k+1} &= F\mathbf{x}_n^k + G\mathbf{x}_n^{k+1} - G\mathbf{x}_n^k \\ &= G\mathbf{x}_n^{k+1} + (F-G)\mathbf{x}_n^k \end{aligned} $$

where F and G are the matrix representation of the fine and coarse propagator.

There are different interpretations of parareal where it can be seen as a type of a multiple shooting method or like a time-multigrid method [12]. A thorough convergence study for the algorithm has been done by Gander and Vandewalle [12]. For a linear system with bounded time domain one gets superlinear convergence but it becomes linear when the time domain is unbounded. For a general nonlinear system the convergence remains superlinear [10].

**3.1. Error analysis.** We now focus on obtaining a general expression of the error for a linear system. We write the correction step (3.3)

$$ \mathbf{x}_n^k = G\mathbf{x}_{n-1}^k + (F-G)\mathbf{x}_{n-1}^{k-1} $$

Let us assume that the parareal reference solution is the solution obtained by integrating the model using the fine solver. Let $\mathbf{x}_f$ be the reference solution at the start of the time window n (i.e. $\mathbf{x}_f = F^n\mathbf{x}_0$). Then the error $\mathbf{e}_n^k = \mathbf{x}_f - \mathbf{x}_n^k$ satisfies the evolution equation

$$ \mathbf{e}_n^k = G\mathbf{e}_{n-1}^k + (F-G)\mathbf{e}_{n-1}^{k-1} $$

and thus

$$ \mathbf{e}_n^k = (F - G) \sum_{p=k}^{n-1} G^{n-p-1} \mathbf{e}_p^{k-1} $$

where we have $\mathbf{e}_0^0 = 0$. If the initial states are obtained by a coarse grid integration (i.e. $\mathbf{x}_n^0 = G^n\mathbf{x}_0$), we have $\mathbf{e}_n^0 = (F - G)\mathbf{x}_n^0$ and we get a general expression for $\mathbf{e}_n^k$:

$$ \mathbf{e}_n^k = \mathbf{x}_0 \sum_{p=k+1}^{n} C_p^n (F-G)^p G^{n-p}, \quad C_p^n = \frac{n!}{p!(n-p)!} $$

Since our cost function (2.7) involves an observation at the end of the integration time, we can accordingly define an approximate parareal based integrator over $[0, T]$ by

$$ (3.4) \qquad P(k) \mathbf{x}_0 = \mathbf{x}_N^k $$
**4.1 Formulation** The inexact conjugate gradient method (hereafter referred to as *inexact CG*) provides a way of allowing inaccuracies in matrix-vector multiplications when solving convex optimisation problem. The method monitors the decrease in the quadratic $q(\delta\mathbf{x}) = \frac{1}{2} \delta\mathbf{x}^T \mathbf{A} \delta\mathbf{x} - \mathbf{b}^T \delta\mathbf{x}$. The quadratic change is directly linked to the energy norm of the residual $\mathbf{r}(\delta\mathbf{x}_j) = \mathbf{A}\delta\mathbf{x}_j - \mathbf{b}$ [20, §2], offering an improved stopping minimisation criterion in terms of avoiding under- or over-solving.

Let $\mathbf{E}_j$ denote the error in the matrix $\mathbf{A}$ at inexact CG iteration (hereafter referred to as *icg-iteration*) $j$ such that for corresponding conjugate direction vector $\mathbf{p}_j$ we have access to the matrix-vector product $(\mathbf{A} + \mathbf{E}_j)\mathbf{p}_j$. The presence of error $\mathbf{E}_j$ results in the computed residual $\mathbf{r}_j$ being different from $\mathbf{r}(\delta\mathbf{x}_j)$. Let the size of perturbation matrix size $\mathbf{E}_j$ be measured in the primal-dual norm defined as

$$ (4.1) \qquad \|\mathbf{E}_j\|_{\mathbf{A}^{-1}, \mathbf{A}} = \sup_{\delta\mathbf{x} \neq 0} \frac{\|\mathbf{E}_j \delta\mathbf{x}\|_{\mathbf{A}^{-1}}}{\|\delta\mathbf{x}\|_{\mathbf{A}}} = \|\mathbf{A}^{-1/2} \mathbf{E}_j \mathbf{A}^{-1/2}\|_2 $$

For some permissible bound on the error norm $\|\mathbf{E}_j\|_{\mathbf{A}^{-1}, \mathbf{A}}$ at icg-iteration $j$, if the inexact residual norm $\|\mathbf{r}_j\|_{\mathbf{A}^{-1}}$ and the residual gap norm $\|\mathbf{r}(\delta\mathbf{x}_j) - \mathbf{r}_j\|_{\mathbf{A}^{-1}}$ can be suitably bounded, then the change in the quadratic can be controlled. The following theorem captures the essence of this idea.

**THEOREM 4.1** (From [20]). Let $\epsilon > 0$ and let $\phi = (\phi_0, \phi_1, \dots, \phi_{j-1})^T \in \mathbb{R}^j$ be a positive vector satisfying

$$ (4.2) \qquad \sum_{i=1}^{j} \frac{1}{\phi_i} \le 1 $$

Suppose that

$$ (4.3) \qquad \|\mathbf{E}_i\|_{\mathbf{A}^{-1}, \mathbf{A}} \le \omega_i = \frac{\sqrt{\epsilon} \|\mathbf{b}\|_{\mathbf{A}^{-1}} \|\mathbf{p}_i\|_{\mathbf{A}}}{2 \phi_{i+1} \|\mathbf{r}_i\|_2^2 + \sqrt{\epsilon} \|\mathbf{b}\|_{\mathbf{A}^{-1}} \|\mathbf{p}_i\|_{\mathbf{A}}} $$

for all $i \in \{0, \dots, j-1\}$, with $\mathbf{p}_i$ being the conjugate direction at icg-iteration $i$. Then

$$ (4.4) \qquad \|\mathbf{r}(\delta\mathbf{x}_j) - \mathbf{r}_j\|_{\mathbf{A}^{-1}} \le \frac{\sqrt{\epsilon}}{2} \|\mathbf{b}\|_{\mathbf{A}^{-1}} $$

Additionally if

$$ (4.5) \qquad \|\mathbf{r}_j\|_{\mathbf{A}^{-1}} \le \frac{\sqrt{\epsilon}}{2} \|\mathbf{b}\|_{\mathbf{A}^{-1}} $$

then

$$ (4.6) \qquad |q(\delta\mathbf{x}_j) - q(\delta\mathbf{x}_*)| \le \epsilon |q(\delta\mathbf{x}_*)|. $$

**Remark 4.2.** Due to increasing error in the matrix-vector product, search directions $\mathbf{p}_j$ are no longer conjugate to each other and the (inexact) residuals lose their orthogonality. While the theorem assumes that the inexact residual norm $\|\mathbf{r}_j\|_{\mathbf{A}^{-1}}$ eventually becomes smaller, it is not guaranteed. Adding a reorthogonalisation step could ensure that the residuals converge to zero after at most $n$ steps.
and thus

$$ (3.5) \qquad F^N - P(k) = \sum_{p=k+1}^{N} C_p^N (F - G)^p G^{N-p} $$

**3.2. Speedup.** We assume that the we have N number of processors available so that each processor can be assigned to one time window. We also assume that the communication time between any two processors is so small that it can be ignored.

The theoretical speedup $\psi$ can be defined as the ratio of the time taken to perform a task sequentially to the time to do the same task in parallel.

$$ (3.6) \qquad \psi = \frac{\text{Serial computation time}}{\text{Parallel computation time}} $$

To derive the speedup for the parareal algorithm we follow along the lines of [20]. We begin by defining some quantities. Let $N_f$ and $N_g$ denote the number of fine and coarse time steps per time window respectively. Let $\tau_f$ be the time taken by the processor to perform a single step of the numerical scheme for the fine propagator $F$ and similarly $\tau_g$ for the coarse propagator $G$. Thus for a single time window, $N_f\tau_f$ is the total cost of $F$ and $N_g\tau_g$ for $G$. To simplify notations, let $\gamma_f = N_f\tau_f$ and $\gamma_g = N_g\tau_g$.

We take the solution from the fine propagator as our reference solution. So the serial integration time is nothing but the time taken by the fine solver to solve the problem sequentially. Since the initial configuration for the parareal is fed through a serial coarse integration, it is given by $N\gamma_g$. Now one iteration of parareal involves a prediction by $G$ which is again sequential and the correction step in parallel. The total cost per parareal iteration is given as $N\gamma_g + \gamma_f$. Therefore for $k$ p-iterations,

$$ \begin{aligned} \psi &= \frac{N\gamma_f}{N\gamma_g + k(N\gamma_g + \gamma_f)} \\ &= \frac{\gamma_f}{\gamma_g + k(\gamma_g + \frac{\gamma_f}{N})} \end{aligned} $$

Let us define $\beta = \frac{\gamma_g}{\gamma_f}$, we have

$$ \psi = \frac{1}{\beta + k \left( \beta + \frac{1}{N} \right)} $$

As stated in [26] we get different bounds for the speedup and they compete with each other.

$$ \psi \le \min \left( \frac{N}{k}, \frac{\text{fine time}}{\text{coarse time}} \right) $$

The coarse propagator should be really cheap and fast in order to get a good speedup but at the same time there will be more p-iterations which will reduce the speedup according to the first bound. Thus the optimal speedup must balance the two bounds in the best possible way. For a more rigorous explanation about the speedup and efficiency of the parareal see [2].

**3.3. Coupling.** We now replace the parallel-in-time operator $P$ in (2.7) by the parareal algorithm. In this case the sequential model trajectories are obtained by the reference solution which is nothing but the solution from the fine propagator $F$. We now look at the cost function (2.6)

$$ J(\mathbf{x}_0) = \frac{1}{2} \|F^N \mathbf{x}_0 - \mathbf{y}\|_2^2 $$

and its gradient

$$ \nabla J(\mathbf{x}_0) = (F^N)^T (F^N \mathbf{x}_0 - \mathbf{y}) $$

Again instead of using the reference solution for the forward integration, if we use the parareal based integrator $P(k)$, (2.6) becomes
**4.2 Controlling error with Parareal stopping criteria** In Theorem 4.1, the residual gap is kept under controlled precision provided that $||\mathbf{E}_j||_{\mathbf{A}^{-1},\mathbf{A}}$ is bounded by the quantity $\omega_j$. In our case the $\mathbf{E}$ matrix, by construction, at any Parareal iteration $k$ is given as

$$ (4.7) \qquad \begin{aligned} \mathbf{E}(k) &= \mathbf{A}_{\text{para}} - \mathbf{A} \\ &= (\mathbf{F}^N)^T \mathbf{P}(k) - (\mathbf{F}^N)^T \mathbf{F}^N \\ &= -(\mathbf{F}^N)^T (\mathbf{F}^N - \mathbf{P}(k)) \end{aligned} $$

It can be said that $\omega_j$ acts as a threshold for how large an error can be allowed at a particular icg-iteration $j$ by controlling the number of Parareal iterations $k$. As it will be seen in our numerical experiments (section 5), the use of $||\mathbf{E}_j(k)||_{\mathbf{A}^{-1},\mathbf{A}}$ does not provide a refined enough stopping criterion for our application. We show in the following that it is preferable to bound the product of the error matrix $\mathbf{E}$ with the vector $\mathbf{p}_j$ and that an estimate for this product can be easily obtained in the context of Parareal. In the proof of Theorem 4.1, the bound (4.4) is obtained by using the inequality

$$ (4.8) \qquad ||\mathbf{E}_j \mathbf{p}_j||_{\mathbf{A}^{-1}} \le ||\mathbf{E}_j||_{\mathbf{A}^{-1},\mathbf{A}} ||\mathbf{p}_j||_{\mathbf{A}} \le \omega_j ||\mathbf{p}_j||_{\mathbf{A}} $$

Upon observation, we found that utilising $||\mathbf{E}_j \mathbf{p}_j||_{\mathbf{A}^{-1}}$ instead of $||\mathbf{E}_j||_{\mathbf{A}^{-1},\mathbf{A}}$ for the Parareal stopping criterion gives much better results. To bound our new norm we introduce a new quantity $\xi_j$ which satisfies,

$$ (4.9) \qquad ||\mathbf{E}_j \mathbf{p}_j||_{\mathbf{A}^{-1}} \le \omega_j ||\mathbf{p}_j||_{\mathbf{A}} = \xi_j $$

Thus from (4.3) the value of $\xi_j$ can be obtained as

$$ (4.10) \qquad \xi_j = \frac{\sqrt{\epsilon} ||\mathbf{b}||_{\mathbf{A}^{-1}} ||\mathbf{p}_j||_{\mathbf{A}}^2}{2\phi_{j+1} ||\mathbf{r}_j||_2^2 + \sqrt{\epsilon} ||\mathbf{b}||_{\mathbf{A}^{-1}} ||\mathbf{p}_j||_{\mathbf{A}}} $$

The following theorem encapsulates the modifications discussed above. The proof is almost identical to the proof of theorem 1 in [20], and uses lemma 2 of the same paper.

**THEOREM 4.3.** In Theorem 4.1, condition (4.3) is replaced by

$$ (4.11) \qquad ||\mathbf{E}_j(k)||_{\mathbf{A}^{-1},\mathbf{A}} = \sup_{\delta\mathbf{x}\neq\mathbf{0}} \frac{||\mathbf{E}_j(k)\delta\mathbf{x}||_{\mathbf{A}^{-1}}}{||\delta\mathbf{x}||_{\mathbf{A}}} \le \xi_j/||\mathbf{p}_j||_{\mathbf{A}} $$

However, using the fact that

$$ (4.12) \qquad ||\mathbf{E}_j(k)\mathbf{p}_j||_{\mathbf{A}^{-1}} \le ||\mathbf{E}_j(k)||_{\mathbf{A}^{-1},\mathbf{A}} ||\mathbf{p}_j||_{\mathbf{A}} $$

we can get a tighter bound, while retaining the validity of the theorem.

**4.3 Other feasible approximations** So far the method has been discussed from a theoretical point of view with an assumption that the quantities or norms concerning the matrix $\mathbf{A}$ and its inverse are accessible. But realistically one might not even have access to the matrix $\mathbf{A}$, let alone its inverse. Thus in practice Gratton et al. [20, §3.2-3.4] proposed the following approximations:

$$ \bullet \quad ||\mathbf{p}_j||_{\mathbf{A}} \approx \sqrt{\frac{1}{n} \text{Tr}(\mathbf{A})} ||\mathbf{p}_j||_2 $$

$$ \bullet \quad q(\delta\mathbf{x}_j) \approx q_j \stackrel{\text{def}}{=} -\frac{1}{2}\mathbf{b}^T \delta\mathbf{x}_j $$
$$ (3.7) \qquad J(\mathbf{x}_0) \approx \frac{1}{2} \|P(k)\mathbf{x}_0 - \mathbf{y}\|_2^2 $$

In this respect, the gradient can be written as

$$ \nabla J(\mathbf{x}_0) \approx (F^N)^T (P(k)\mathbf{x}_0 - \mathbf{y}) $$

The corresponding linear system $Ax = b$ to solve has matrix $A = (F^N)^T P(k)$ and the vector $b = (F^N)^T y$

**REMARK 3.1.** Since the fine propagator is replaced by the parareal algorithm in (3.7), we will denote the exact matrix by $A_*$ and the inexact (approximate) matrix by $A$. Thus $A_* = (F^N)^T F^N$ and $A = (F^N)^T P(k)$.

**4. Minimisation.** For solving large linear systems, choosing a suitable Krylov subspace method is a logical option [14]. Krylov methods are iterative in nature with the matrix-vector multiplication as their most dominant operation. In some situations, there might be the case when the matrix-vector product is not exact because the matrix comes from a discretisation scheme of a differential equation. Another reason could be the lack of access to the matrix, and it is only available as a product by a vector.

However, under certain conditions, this undesirable inexactness could be used to better effects in terms of some efficiency gains. In literature, such methods are known as the inexact Krylov subspace methods. They are based on a counterintuitive principle that the error in the matrix-vector multiplication is permitted to grow as we move ahead with the iterations [27, 30]. By doing so, the same level of accuracy can be maintained as expected by using the usual Krylov subspace methods.

In our case we want to minimise a quadratic cost function of the data assimilation problem which leads to solving a symmetric system $Ax = b$. As seen before if parareal is used for the forward integration then $A = (F^N)^T P(k)$ and $b = (F^N)^T y$. $P(k)$ is an operator which acts on the initial condition to give the parareal solution at the end of the integration time. Thus, the matrix-vector multiplication is provided by the parareal method which is itself iterative and bound to some stopping tolerance. We are going to use the inexact conjugate gradient method which is effective for linear symmetric systems.

**4.1. Inexact conjugate gradient.** We can rewrite the 4D-Var cost function in terms of $A_*$ and $b$ leading to the following convex quadratic minimisation problem

$$ \min_{\mathbf{x} \in \mathbb{R}^n} J(\mathbf{x}) = \frac{1}{2} \mathbf{x}^T A_* \mathbf{x} - \mathbf{b}^T \mathbf{x} $$

with symmetric positive definite matrix $A_* \in \mathbb{R}^{n \times n}$ and $b \in \mathbb{R}^n$

As mentioned before, this problem can also be seen as solving the symmetric positive-definite linear system $A_* \mathbf{x} = b$ by a suitable krylov subspace method. We would like to look for any efficiency gains in the sense that relaxing the accuracy of the matrix-vector products (thus allowing inexactness) would still give the same levels of efficiency.

The Inexact Conjugate Gradient (inexact CG) method proposed by Gratton et al. [13] monitors the decrease in the quadratic when the matrix-vector multiplications are inexact. The reason behind focusing on the quadratic is the direct relationship of its change with the energy norm of the residual which leads to a better stopping criterion. If $\mathbf{r}(\mathbf{x}) = \mathbf{b} - A_* \mathbf{x}$ is defined to be the residual of the system then this relationship can be described as [13]

$$ (4.1) \qquad \frac{1}{2} \|\mathbf{r}(\mathbf{x})\|_{A_*^{-1}}^2 = J(\mathbf{x}) - J(\mathbf{x}_*) $$

where $\mathbf{x}_* = A_*^{-1} \mathbf{b}$ is taken as the exact solution of the system.

We will adopt this idea to deal with the inexactness introduced by the use of a parareal method. Let us represent the source of inexactness in the matrix $A$ by the perturbation matrix $E$. In our context the error manifests from the parareal method in equation (3.5). By construction, the $E$ matrix at any
$$ \bullet \quad \|\mathbf{b}\|_{\mathbf{A}^{-1}} = \sqrt{2|q(\delta\mathbf{x}_j)|} \approx \begin{cases} \frac{\|\mathbf{b}\|_2}{\mu_{\max}(\mathbf{A})}, & j = 0 \\ \sqrt{2|q_j|}, & j = 1, \dots, j_{\max} \end{cases} $$

• Termination criterion (4.5) by $q_{j-d} - q_j \le \frac{1}{4} \epsilon |q_j|$ for some integer $d$

Remark 4.4. The quality of approximation for $q$ has been discussed in Sec 3.3 of Gratton et al. In the presence of matrix-vector product errors, the equality

$$ (4.13) \qquad q(\delta\mathbf{x}_j) = -\frac{1}{2}\mathbf{b}^T\delta\mathbf{x}_j $$

may no longer hold. However, $q(\delta\mathbf{x}_j)$ remains close to $q_j$ if the inexactness of the matrix-vector products are controlled in the same way as in Theorem 4.1 (see Sec 3.3, Theorem 2).

Since we have a modified stopping criterion involving $\mathbf{A}^{-1}$ norm, we also require an approximation for $\|\mathbf{E}_j\mathbf{p}_j\|_{\mathbf{A}^{-1}}$

**4.4 Approximation to $\|\mathbf{E}_j\mathbf{p}_j\|_{\mathbf{A}^{-1}}$** We found that we can replace $\|\mathbf{E}_j\mathbf{p}_j\|_{\mathbf{A}^{-1}}$ by another equivalent and computable quantity as proved in the theorem below.

**THEOREM 4.5.** If $\mathbf{F}$ is invertible, $\|\mathbf{E}_j(k)\mathbf{p}_j\|_{\mathbf{A}^{-1}} = \|\mathbf{P}(k)\mathbf{p}_j - \mathbf{F}^N\mathbf{p}_j\|_2$ irrespective of the choice of Parareal approximation. $\mathbf{P}(k)\mathbf{p}_j$ is the Parareal approximation at iteration $k$ with $\mathbf{p}_j$ as the initial condition and $\mathbf{F}^N\mathbf{p}_j$ is the corresponding exact solution.

*Proof.* We have,

$$ \begin{aligned} \|\mathbf{E}_j(k)\mathbf{p}_j\|_{\mathbf{A}^{-1}} &= \langle \mathbf{E}_j(k)\mathbf{p}_j, \mathbf{A}^{-1}\mathbf{E}_j(k)\mathbf{p}_j \rangle \\ &= \langle (\mathbf{F}^N)^T \mathbf{F}^N - (\mathbf{F}^N)^T \mathbf{P}(k)] \mathbf{p}_j, [(\mathbf{F}^N)^T \mathbf{F}^N]^{-1} \\ & \qquad \qquad [(\mathbf{F}^N)^T \mathbf{F}^N - (\mathbf{F}^N)^T \mathbf{P}(k)] \mathbf{p}_j \rangle \\ &= \langle (\mathbf{F}^N)^T [(\mathbf{F}^N - \mathbf{P}(k))\mathbf{p}_j], [(\mathbf{F}^N)^T \mathbf{F}^N]^{-1} (\mathbf{F}^N)^T [(\mathbf{F}^N - \mathbf{P}(k))\mathbf{p}_j] \rangle \\ &= \langle (\mathbf{F}^N - \mathbf{P}(k))\mathbf{p}_j, \mathbf{F}^N[(\mathbf{F}^N)^T \mathbf{F}^N]^{-1} (\mathbf{F}^N)^T [(\mathbf{F}^N - \mathbf{P}(k))\mathbf{p}_j] \rangle \\ &= \langle (\mathbf{F}^N - \mathbf{P}(k))\mathbf{p}_j, [\mathbf{F}^N(\mathbf{F}^N)^\dagger] [(\mathbf{F}^N - \mathbf{P}(k))\mathbf{p}_j] \rangle \end{aligned} $$

where $(\mathbf{F}^N)^\dagger$ is the Moore-Penrose generalised inverse or pseudo-inverse of $\mathbf{F}^N$. Thus when $\mathbf{F}$ is invertible (i.e. $(\mathbf{F}^N)^\dagger = (\mathbf{F}^N)^{-1}$) we have

$$ \begin{aligned} \|\mathbf{E}_j(k)\mathbf{p}_j\|_{\mathbf{A}^{-1}} &= \langle (\mathbf{F}^N - \mathbf{P}(k))\mathbf{p}_j, (\mathbf{F}^N - \mathbf{P}(k))\mathbf{p}_j \rangle \\ &= \langle \mathbf{F}^N\mathbf{p}_j - \mathbf{P}(k)\mathbf{p}_j, \mathbf{F}^N\mathbf{p}_j - \mathbf{P}(k)\mathbf{p}_j \rangle \\ &= \|\mathbf{P}(k)\mathbf{p}_j - \mathbf{F}^N\mathbf{p}_j\|_2 \end{aligned} \quad \square $$

In the general case there is no guarantee that $\mathbf{F}$ is invertible so we can only say that $\|\mathbf{P}(k)\mathbf{p}_j - \mathbf{F}^N\mathbf{p}_j\|_2$ may be a good approximation for $\|\mathbf{E}_j\mathbf{p}_j\|_{\mathbf{A}^{-1}}$.

We managed to replace a quantity which involves the $\mathbf{A}^{-1}$ norm to a quantity which needs the 2-norm. Obviously, what remains is to find an approximation for $\mathbf{F}^N\mathbf{p}_j$ since running the fine solver $\mathbf{F}$ each time is computationally not feasible. The idea is to approximate $\|\mathbf{E}_j\mathbf{p}_j\|_{\mathbf{A}^{-1}}$ as

$$ (4.14) \qquad \|\mathbf{E}_j\mathbf{p}_j\|_{\mathbf{A}^{-1}} = \|\mathbf{P}(k)\mathbf{p}_j - \mathbf{F}^N\mathbf{p}_j\|_2 \approx \|\mathbf{P}(k)\mathbf{p}_j - \mathbf{P}(k+1)\mathbf{p}_j\|_2 $$
parareal iteration k can be written as

$$
\begin{align*}
E(k) &= A - A_* \\
&= (F^N)^T P(k) - (F^N)^T F^N \\
&= -(F^N)^T (F^N - P(k))
\end{align*}
$$

Using (3.5) we get,

$$
(4.2) \qquad E(k) = -(F^N)^T \sum_{p=k+1}^{N} C_p^N (F-G)^p G^{N-p}
$$

From (4.1) it is clear that a sufficiently accurate residual **r**(**x**) and its $A_*^{-1}$ norm is required to monitor
||**r**(**x**)‖<sub>A<sup>-1</sup></sub>. But the presence of errors at each inexact CG iteration (icg-iteration) *j* also results in the
computed residual **r**<sub>*j*</sub> being different from the exact residual **r**(**x**<sub>*j*</sub>) = **b** − *A*<sub>*</sub>**x**<sub>*j*</sub>. In order to still keep
track of the quadratic change, it is necessary to appropriately bound the inexact residual norm ||**r**<sub>*j*</sub>||<sub>A<sup>-1</sup></sub>
and the residual gap norm ||**r**(**x**<sub>*j*</sub>) − **r**<sub>*j*</sub>||<sub>A<sup>-1</sup></sub>. We are also going to need the primal-dual norm of the
perturbation matrix defined as

$$
\|E_j\|_{A^{-1}, A} = \sup_{\mathbf{x} \neq 0} \frac{\|E_j \mathbf{x}\|_{A^{-1}}}{\|\mathbf{x}\|_A} = \|A^{-1/2} E_j A^{-1/2}\|_2
$$

We have dropped the variable *k* while defining the error norm for clarity. So in the later sections it
will be present only when it is explicitly required in the context. Now for some permissible error at
icg-iteration *j*, if the inexact residual norm and the residual gap norm can be suitably bounded, then
the change in the quadratic can be controlled. The following theorem captures the essence of the above
idea [13], where **p** is taken as the conjugate direction vector.

THEOREM 4.1. Let $\epsilon > 0$ and let $\phi \in \mathbb{R}^j$ be a positive vector satisfying

$$
\sum_{i=1}^{j} \frac{1}{\phi_i} \le 1
$$

Suppose that

$$
(4.3) \qquad \|E_i\|_{A^{-1}, A} \le \omega_i = \frac{\sqrt{\epsilon} \|\mathbf{b}\|_{A^{-1}} \|\mathbf{p}_i\|_A}{2\phi_{i+1} \|\mathbf{r}_i\|_2^2 + \sqrt{\epsilon} \|\mathbf{b}\|_{A^{-1}} \|\mathbf{p}_i\|_A}
$$

for all $i \in \{0, \dots, j-1\}$. Then

$$
(4.4) \qquad \| \mathbf{r}(\mathbf{x}_j) - \mathbf{r}_j \|_{A^{-1}} \le \frac{\sqrt{\epsilon}}{2} \| \mathbf{b} \|_{A^{-1}}
$$

Additionally if

$$
(4.5) \qquad \| \mathbf{r}_j \|_{A^{-1}} \le \frac{\sqrt{\epsilon}}{2} \| \mathbf{b} \|_{A^{-1}}
$$

then $|J(\mathbf{x}_j) - J(\mathbf{x}_*)| \le \epsilon |J(\mathbf{x}_*)|$.

REMARK 4.2. The errors are permitted to grow as we proceed with the icg-iterations. As a consequence the search directions are no longer conjugate to each other and the (inexact) residuals lose their orthogonality. In the theorem it is assumed that the inexact residual norm $||\mathbf{r}_j||_{A^{-1}}$ eventually becomes smaller but it is not guaranteed. Adding a reorthogonalisation step could ensure that the residuals converge to zero after at most *n* steps.

REMARK 4.3. Since the inexactness is introduced by another iterative method (parareal in this case), for each outer minimisation iteration *j* there are many inner iterations *k* of the parareal. The number of outer iterations depend on the stopping criteria (4.5) and the number of inner iterations on $\omega_j$ in equation (4.3).
This is achieved by computing the difference between successive Parareal iterates in the 2-norm. As the Parareal iteration $k$ increases the iterate $\mathbf{P}(k+1)\mathbf{p}_j$ gets closer and closer to $\mathbf{F}^N\mathbf{p}_j$. Consequently, we only have to run Parareal to a point where the Parareal stopping criterion is met (bounded by the value of $\xi_j$). That is, at each icg-iteration $j$ we increment the Parareal iterations $k$ by 1 until the condition $\|\mathbf{P}(k)\mathbf{p}_j - \mathbf{P}(k+1)\mathbf{p}_j\|_2 \le \xi_j$ is satisfied. However, as a trade-off for using an approximation, we end up using one more Parareal iteration than needed since the latest Parareal iterate $\mathbf{P}(k+1)\mathbf{p}_j$ replaces $\mathbf{F}^N\mathbf{p}_j$ to compute the difference with the last Parareal iterate $\mathbf{P}(k)\mathbf{p}_j$ (see Algorithm 4.1).

<caption>**Algorithm 4.1 Approximation for $\mathbf{F}^N\mathbf{p}_j$**</caption>

1: **Given:** icg-iteration $j$, direction vector $\mathbf{p}_j$
2: Run 2 Parareal iterations for initial state $\mathbf{p}_j$
3: Set $k = 1$
4: **while** true **do**
5: Compute $e = \|\mathbf{P}(k)\mathbf{p}_j - \mathbf{P}(k+1)\mathbf{p}_j\|_2$
6: **if** $e > \xi_j$ **then**
7: Run one more Parareal iteration and set $k = k+1$
8: **else**
9: **break**
10: **end if**
11: **end while**
12: Store the Parareal iterate $\mathbf{P}(k)\mathbf{p}_j$ for computing the matrix-vector product

**4.5 Inaccuracy budget** Until now we have not talked about the role of $\phi_j$ in Theorem 4.1, which can be used to manage the inaccuracy budget [20, §3.1]. By inaccuracy budget we mean adjusting the error tolerance $\omega_j$ for icg-iteration $j$ by distributing the *unused inaccuracy* between the error bound $\xi_j$ and the actual error norm $\|\mathbf{E}_j\mathbf{p}_j\|_{\mathbf{A}^{-1}}$ to subsequent icg-iterations. This allows for a larger error in the matrix-vector product and thus using fewer number of Parareal iterations.

Suppose for a given $\phi_{j+1}$ we obtain $\|\mathbf{E}_j\mathbf{p}_j\|_{\mathbf{A}^{-1}} = \hat{\xi}_j \le \xi_j$. The corresponding $\xi_j$ can be obtained from $\hat{\xi}_j(\hat{\phi}_{j+1})$ in (4.10) as

$$ (4.15) \qquad \hat{\phi}_{j+1} = \frac{(\|\mathbf{p}_j\|_{\mathbf{A}} - \hat{\xi}_j) \sqrt{\epsilon} \|\mathbf{b}\|_{\mathbf{A}^{-1}} \|\mathbf{p}_j\|_{\mathbf{A}}}{2 \|\mathbf{r}_j\|_2^2} > \phi_{j+1} $$

We distribute the inaccuracy evenly in the remaining $j_{max} - j - 1$ icg-iterations the same way as in the original method [20] with $j_{max}$ being the maximum icg-iterations allowed. Algorithm 4.2 implements inexact CG using all the manipulations done by including the inaccuracy budget, other feasible approximations, and Parareal stopping criterion involving $\|\mathbf{E}_j\mathbf{p}_j\|_{\mathbf{A}^{-1}}$. In the next section we are going to analyse the results of the Parareal and data assimilation coupling using various settings of the minimisation algorithm.

**5 Numerical experiments** We remind the reader that the primary goal of our approach is to improve the theoretical speedup (as introduced in subsection 2.3.1) of the inner loop in incremental 4D-Var using Parareal. This is accomplished by using a modified Parareal stopping criterion within inexact CG.

In this section, we numerically demonstrate how the theoretical speedup is improved by running a series of minimisation experiments with different versions of CG and inexact CG. Their performance is compared in terms of the total number of
**4.2. Stopping criterion.** In Theorem 4.1, the quantity $\omega_j$ acts as a threshold for how large a perturbation can be allowed at a particular icg-iteration $j$. In the beginning one starts with a small value of $\omega_j$ with more and more error allowed (and thus gradual increase in $\omega_j$) in the subsequent iterations. When $\kappa(A)$ is large, the primal-dual norm of the perturbation matrix $\|E_j\|_{A^{-1},A}$ is quite large and it usually goes down below 1 only at the last few p-iterations. In terms of efficiency the number of p-iterations being used in criterion (4.3) does not reflect upon the actual number of p-iterations needed. We make a few modifications in the original method to fix this issue. In the proof of Theorem 4.1, the bound (4.4) is obtained using the inequality

$$ \|E_j \mathbf{p}_j\|_{A^{-1}} \le \|E_j\|_{A^{-1},A} \|\mathbf{p}_j\|_{A^{-1}} \le \omega_j \|\mathbf{p}_j\|_{A^{-1}} $$

We noticed that utilising $\|E_j \mathbf{p}_j\|_{A^{-1}}$ instead of $\|E_j\|_{A^{-1},A}$ for the parareal stopping criterion gives us fairly reasonable results. It can be clearly seen that this in no way affect the original method. To bound our new norm we introduce a quantity $\xi_j$ which satisfies,

$$ \|E_j \mathbf{p}_j\|_{A^{-1}} \le \omega_j \|\mathbf{p}_j\|_{A^{-1}} = \xi_j $$

Thus from (4.3),

$$ (4.6) \qquad \xi_j = \frac{\sqrt{\epsilon} \|\mathbf{b}\|_{A^{-1}} \|\mathbf{p}_j\|_A^2}{2\phi_{j+1} \|\mathbf{r}_j\|_2^2 + \sqrt{\epsilon} \|\mathbf{b}\|_{A^{-1}} \|\mathbf{p}_j\|_A} $$

**Note:** $\|E_j \mathbf{p}_j\|_{A^{-1}}$ is the product of the error matrix and the conjugate direction vector $\mathbf{p}_j$ at a given parareal iteration. This means that $\mathbf{p}_j$ is the initial condition for the parareal algorithm.

**REMARK 4.4.** The modified criterion does not make use of $\|E_j\|_{A^{-1},A}$ anywhere in the algorithm. Thus, we don't need to know the perturbation matrix explicitly. It is sufficient for us if we know it indirectly in the form of a matrix vector product. With a given conjugate direction $\mathbf{p}_j$ the product $P(k)\mathbf{p}_j - A\mathbf{p}_j$ is nothing but $E_j(k)\mathbf{p}_j$.

Figure 4.1 below shows the number of p-iterations for a shallow water model if we use the original criterion and the modified criterion for a large $\kappa(A)$. The red marked symbols indicate the values on the either side of the tolerance. It can be seen that there is a huge difference between the two criteria. On the left image we see that the bound is satisfied only at the second last iteration whereas on the right image it is already satisfied in 5 iterations.

<figure><img src="image_4.png" /><figcaption>FIG. 4.1. Estimates and bounds for icg-iteration 6, $N = 20, \kappa(A) = 1099070.55$</figcaption></figure>

From the figure on the left the curve resembles to the kind of parareal convergence results as reported in [8] for pure advection problems. It says that even if the bound is superlinear the constants in the bounds
<caption>**Algorithm 4.2** Inexact CG with all approximations</caption>

1: **Given:** Symmetric positive definite matrix $\mathbf{A} \in \mathbb{R}^{n \times n}$, right hand side vector $\mathbf{b} \in \mathbb{R}^n$, tolerance $\epsilon_{icg}$
2: Set $\delta\mathbf{x}_0 = 0$, $\mathbf{r}_0 = -\mathbf{b}$, $\mathbf{p}_0 = \mathbf{r}$, $\beta_0 = \|\mathbf{b}\|_2^2$, $\mathbf{u}_1 = \mathbf{b}/\beta_0$, $\phi_0 = j_{\max}$, $\Phi_0 = 1$
3: **for** $j = 0, \dots, j_{\max}$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;Compute the approximation to $\|\mathbf{p}_j\|_{\mathbf{A}}$
5: &nbsp;&nbsp;&nbsp;&nbsp;Determine $\xi_j$ from the equation (4.10)
6: &nbsp;&nbsp;&nbsp;&nbsp;**for** $k = 1, \dots, N-1$ **do**
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Run Parareal with $k+1$ iterations
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Calculate $e = \|\mathbf{P}(k+1)\mathbf{p}_j - \mathbf{P}(k)\mathbf{p}_j\|_2$
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $e < \xi_j$ **then**
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Set $\hat{\xi}_j = e$
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**break**
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end if**
13: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
14: Compute the product $\mathbf{c}_j = (\mathbf{F}^N)^T \mathbf{P}(k)\mathbf{p}_j$
15: Find $\hat{\phi}_j$ from $\hat{\xi}_j$
16: $\Phi_{j+1} = \Phi_j - \hat{\phi}_j^{-1}$
17: **if** $j < j_{\max}$ **then**
18: &nbsp;&nbsp;&nbsp;&nbsp;$\phi_{j+1} = (j_{\max} - j)/\Phi_{j+1}$
19: **else**
20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\phi_{j+1} = \phi_j$
21: **end if**
22: $\alpha_j = \beta_j / \mathbf{p}_j^T \mathbf{c}_j$
23: $\delta\mathbf{x}_{j+1} = \delta\mathbf{x}_j + \alpha_j \mathbf{p}_j$
24: $\mathbf{r}_{j+1} = \mathbf{r}_j + \alpha_j \mathbf{c}_j$
25: $J_{j+1} = \frac{1}{2} \mathbf{b}^T \delta\mathbf{x}_{j+1}$
26: **if** $(J_{j+1-d} - J_{j+1}) \le \frac{1}{4} \epsilon_{icg} |J_{j+1}|$ **then**
27: **break**
28: **end if**
29: **if** (reorth) **then**
30: &nbsp;&nbsp;&nbsp;&nbsp;**for** $i = 1, \dots, j$ **do**
31: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{r}_{j+1} = \mathbf{r}_{j+1} - (\mathbf{u}_i^T \mathbf{r}_{j+1}) \mathbf{u}_i$
32: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
33: $\beta_{j+1} = \mathbf{r}_{j+1}^T \mathbf{r}_{j+1}$
34: $\mathbf{u}_{j+1} = \mathbf{r}_{j+1} / \sqrt{\beta_{j+1}}$
35: **else**
36: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\beta_{j+1} = \mathbf{r}_{j+1}^T \mathbf{r}_{j+1}$
37: **end if**
38: $\mathbf{p}_{j+1} = -\mathbf{r}_{j+1} + (\beta_{j+1}/\beta_j)\mathbf{p}_j$
39: **end for**

Parareal iterations and the average Parareal iterations per CG/inexact CG iterations.
The results are presented for the following versions:

* CG with exact matrix-vector products.

* CG with inexact matrix-vector products (through Parareal).

* Inexact CG with original Parareal stopping criterion.

* Inexact CG with modified Parareal stopping criterion.

* Inexact CG with all approximations, modified Parareal stopping criterion, (subsections 4.3 and 4.4), and inaccuracy budget (subsection 4.5).

To carry out these experiments, we first set up our numerical model, construct the Parareal propagators **F** and **G**, and define the data assimilation cost function.

**5.1 Numerical setup** Before proceeding, we make the reader aware that we are going to use the notation “**x**” for the increment “$\delta\mathbf{x}$” throughout this section for simplicity. For illustrations we use the linearised one-dimensional shallow water equa-
for numerical method make it hard to get the speedup. The error starts increasing and only comes down at the last iteration.

The relative decrease in the curve of $||E_j p_j||_{A^{-1}}$ can be explained by the fact that even if there is an increase in $||E_j||_{A^{-1},A}$, theoretically the conjugate direction norm $||p_j||_A \to 0$ as $j \to j_{max}$.

Till now there has been no talk about the role of $\phi_j$ which is discussed in the original method [13, see Sec 3.1]. By definition $\phi_{j+1}$ puts a restriction on the bound $\xi_j$ and so choosing a constant value of $\phi_j = j_{max}$ limits the possibility of larger error allowance. In fact the smaller the value of $\phi_{j+1}$, the larger the bound for $\xi$ will be. This can actually be done by managing the inaccuracies obtained by the difference between the computed $\xi_j$ and the allowed inexactness $\xi_j(\text{say})$. This difference is unused and could be utilised by distributing to the subsequent icg-iterations. The corresponding $\hat{\phi}_{j+1}$ can be obtained from $\hat{\xi}_j = \xi_j(\hat{\phi}_{j+1})$ in (4.6) as

$$ (4.7) \qquad \hat{\phi}_{j+1} = \frac{(\|\mathbf{p}_j\|_A - \hat{\xi}_j) \sqrt{\epsilon} \|\mathbf{b}\|_{A^{-1}} \|\mathbf{p}_j\|_A}{2 \|\mathbf{r}_j\|_2^2} > \phi_{j+1} $$

The inaccuracy can be distributed the same way as in the original algorithm.

**4.3. Modified inexact conjugate gradient.** So far we have discussed about the method from a theoretical point of view with an assumption that we have access to the quantities or norms concerning the matrix $A$ and its inverse. But realistically one might not even have access to the matrix $A$, let alone its inverse. Thus in practice we will use the following approximations [13] of all the quantities which are unfeasible to compute.

* $\|\mathbf{p}_j\|_A \approx \sqrt{\frac{1}{n} \mathrm{Tr}(A)} \|\mathbf{p}_j\|_2$

* $J(\mathbf{x}_j) \approx J_j \stackrel{\text{def}}{=} -\frac{1}{2} \mathbf{b}^T \mathbf{x}_j$

* $\|\mathbf{b}\|_{A^{-1}} \approx \frac{\|\mathbf{b}\|_2}{\sqrt{\lambda_{\max}}} \quad j=0, \mathbf{x}_0 = 0 \quad \text{and} \quad \|\mathbf{b}\|_{A^{-1}} \approx \sqrt{2|J_j|} \quad j=1, 2, \dots, j_{\max}$

* Termination test (4.5) by $J_{j-d} - J_j \le \frac{1}{4} \epsilon |q_j|$ for some integer $d$

**Note:** Here $\lambda_{\max}$ is the maximum eigenvalue of the matrix $A$ and should not be confused with the state variable $\mathbf{x}$.

While $||E_j p_j||_{A^{-1}}$ is a suitable stopping criterion for the p-iterations, it is not easy to compute as well. It turns out that it can be replaced by another equivalent and computable quantity.

**THEOREM 4.5.** If $F$ is invertible, $||E_j(k)\mathbf{p}_j||_{A^{-1}}$ and $||P(k)\mathbf{p}_j - \mathbf{p}_*||_2$ are rigorously the same whatever the parareal approximation is. $P(k)\mathbf{p}_j$ is the parareal approximation at iteration $k$ with $\mathbf{p}_j$ as the initial condition and $\mathbf{p}_*$ is the corresponding exact solution.

*Proof.* We have,

$$ \begin{aligned} ||E_j(k)\mathbf{p}_j||_{A^{-1}} &= \langle E_j(k)\mathbf{p}_j, A^{-1}E_j(k)\mathbf{p}_j \rangle \\ &= \langle (F^N)^T F^N - (F^N)^T P(k)] \mathbf{p}_j, [(F^N)^T F^N]^{-1} [(F^N)^T F^N - (F^N)^T P(k)] \mathbf{p}_j \rangle \\ &= \langle (F^N)^T [(F^N - P(k))\mathbf{p}_j], [(F^N)^T F^N]^{-1} (F^N)^T [(F^N - P(k))\mathbf{p}_j] \rangle \\ &= \langle (F^N - P(k))\mathbf{p}_j, F^N[(F^N)^T F^N]^{-1} (F^N)^T [(F^N - P(k))\mathbf{p}_j] \rangle \\ &= \langle (F^N - P(k))\mathbf{p}_j, [F^N(F^N)^\dagger] [(F^N - P(k))\mathbf{p}_j] \rangle \end{aligned} $$

where $(F^N)^\dagger$ is the Moore-Penrose generalised inverse or pseudo-inverse of $F^N$. Thus when $F$ is invertible (i.e. $(F^N)^\dagger = (F^N)^{-1}$) we have
tions [21, 18, 41], a simple model describing some important geophysical phenomena. The equations are given by

$$ (5.1) \qquad \begin{aligned} \frac{\partial \eta}{\partial t} &= -H \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial t} &= -g \frac{\partial \eta}{\partial x} + \mu \frac{\partial^2 u}{\partial x^2} \end{aligned} $$

Here $u$ is horizontal velocity, $\eta$ is free surface elevation, $g$ is gravity, $H$ is a constant reference depth, and $\mu$ is the diffusion constant. An explicit diffusion term is added to account for the unresolved small scales in terms of the resolved large scales also known as the closure of the system [7, 43]. The physical domain $\Lambda = [0, L]$ is closed and the physical values are provided in Table 1.

<table><caption>TABLE 1
Domain and physical parameters for the linearised 1D shallow water model</caption><thead><tr><th>Shallow water parameters</th><th>values</th></tr></thead><tbody><tr><td>Length of the basin</td><td>$L = 120$m</td></tr><tr><td>Reference depth</td><td>$H = 0.9$m</td></tr><tr><td>Gravity</td><td>$g = 10m.s^{-2}$</td></tr><tr><td>Diffusion constant</td><td>$\mu = 0.15m^2.s^{-1}$</td></tr></tbody></table>

### 5.1.1 Spatial discretisation
For our model (5.1) we consider a staggered grid discretisation [10] with $\Delta x$ as the spatial stepping size and $N_x$ as the total number of spatial grid points. Using a second-order centered finite difference scheme in space leads to an ODE system

$$ (5.2) \qquad \frac{dx}{dt} = Cx $$

where

$$ (5.3) \qquad C = \begin{pmatrix} 0 & \begin{array}{c} -H/\Delta x \\ H/\Delta x \end{array} & \begin{array}{c} -H/\Delta x \\ \vdots \\ H/\Delta x \end{array} & \begin{array}{c} -H/\Delta x \\ H/\Delta x \end{array} \\ \begin{array}{c} g/\Delta x & -g/\Delta x \\ g/\Delta x & -g/\Delta x \\ \vdots & \vdots \\ g/\Delta x & -g/\Delta x \end{array} & \begin{array}{c} -2\mu/\Delta x^2 \\ \mu/\Delta x^2 \end{array} & \begin{array}{c} \mu/\Delta x^2 \\ -2\mu/\Delta x^2 \\ \vdots \\ \mu/\Delta x^2 \end{array} & \begin{array}{c} -2\mu/\Delta x^2 \\ \mu/\Delta x^2 \\ \vdots \\ -2\mu/\Delta x^2 \end{array} \end{pmatrix} $$

and $x = (\eta_{1/2}, \dots, \eta_{N_x-1/2}, u_1, \dots, u_{N_x-1})^T$ is the state vector. In the numerical experiments, we will use number of grid points $N_x$ equal to 120 ($\Delta x = L/N_x = 1$m). The dimension of the state vector is thus equal to $n = 2N_x - 1 = 239$.

### 5.1.2 Temporal discretisation and Parareal propagators
Now to solve the linear ODE system (5.2) numerically in time, we use an implicit theta scheme. Let $q$ be the index for the temporal grid and $\Delta T$ be the time-stepping size. We have

$$ (5.4) \qquad \mathbf{x}^{q+1} = [\mathbf{I} - \theta \Delta T \mathbf{C}]^{-1} [\mathbf{I} + (1 - \theta) \Delta T \mathbf{C}] \mathbf{x}^q, \quad q = 0, 1, \dots $$

where $0 \le \theta \le 1$. Next, we construct the fine and coarse propagators **F** and **G** respectively for implementing the Parareal algorithm. Let $\Omega = [0, T]$ be the time
$$ \begin{aligned} \|E_j(k)\mathbf{p}_j\|_{A^{-1}} &= \langle(F^N - P(k))\mathbf{p}_j, (F^N - P(k))\mathbf{p}_j\rangle \\ &= \langle F^N \mathbf{p}_j - P(k)\mathbf{p}_j, F^N \mathbf{p}_j - P(k)\mathbf{p}_j\rangle \\ &= \|P(k)\mathbf{p}_j - \mathbf{p}_*\|_2 \end{aligned} $$

One last hurdle is to find a computable estimate for $\mathbf{p}_*$ since we do not want the parareal algorithm to run with a very small tolerance each time. Looking at Figure 4.2 we notice that the curve profile of $\|E_j\mathbf{p}_j\|_{A^{-1}}$ does not change much for given $j$. Then replacing $\|P(k)\mathbf{p}_j - \mathbf{p}_*\|_2$ by $\|P(k)\mathbf{p}_j - \mathbf{p}_*\|_2$ can still lead to convergence where $\mathbf{p}_*$ is the last parareal iterate. All we need is an accurate parareal solution in the beginning of the inexact CG to ensure that the computed last parareal iterates are reliable as well.

<figure><img src="image_5.png" /><figcaption>FIG. 4.2. Profile of $\|E_j\mathbf{p}_j\|_{A^{-1}}$ for different CG iterations</figcaption></figure>

For a clear explanation, suppose for instance that the initial (accurate) parareal solution $P(k)\mathbf{p}_0$ takes $\bar{k}$ iterations. Then for the next icg-iteration, $\|P(k)\mathbf{p}_1 - \mathbf{p}_*\|_2$ is approximated by $\|P(k)\mathbf{p}_0 - \mathbf{p}_*\|_2$ where $\mathbf{p}_* \approx P(\bar{k})\mathbf{p}_0$. We then find the required $k$ which satisfies $\|P(k)\mathbf{p}_0 - \mathbf{p}_*\|_2 \le \xi_1$, say $k = k_1$. Now the parareal solution $P(k_1)\mathbf{p}_1$ will be used to compute the matrix-vector multiplication. Going to the next iteration, $\|P(k)\mathbf{p}_2 - \mathbf{p}_*\|_2$ will be approximated by $\|P(k)\mathbf{p}_1 - \mathbf{p}_*\|_2$ with $\mathbf{p}_* \approx P(k_1)\mathbf{p}_1$. We proceed similarly for the subsequent minimisation iterations till convergence.

Since we are using the solution from the last parareal iterate we would like to make sure we always have reliable estimates when bounding $\|P(k)\mathbf{p}_j - \mathbf{p}_*\|_2$ with $\xi_j$. For this it is important that the value of $\xi_j$ does not lie between the last two values of $k$ (say). If this happens, then we move on to the current parareal iterate where we compute the solution with the current $\mathbf{p}_j$ as the initial condition. We proceed the same way as before but and check again if the bound $\xi_j$ lies above the last two iterates. If not, we keep doing one more parareal iteration until our condition is satisfied.

Again from the above example, we use the current parareal iterate if $\|P(k_1)\mathbf{p}_0 - \mathbf{p}_*\|_2 \le \xi_1$ but $\xi_1 \in [\bar{k} - 1, \bar{k})$. Like before we run the parareal for $k_1$ iterations to get $P(k_1)\mathbf{p}_1$. What changes now is that we look for $k$ such that $\|P(k)\mathbf{p}_1 - P(k_1)\mathbf{p}_1\|_2 \le \xi_1$. If the required $k$ still falls between the last two values, we perform one more parareal iteration and keep doing it until our condition is satisfied. So if it takes $k_*$ more iterations then we update $k_1 := k_1 + k_*$. This discussion has been put in the form of an algorithm for $\mathbf{p}_*$ approximation which can be implemented.

This test is illustrated in Figure 4.3 for icg-iteration 17. We can see that when for the last parareal
domain such that $T$ is the total time period of integration, $N$ is the total number of time windows and $\Delta T$ is the length of one time window. For each time window, let $N_{\text{fine}}$ and $N_{\text{coarse}}$ be the number of fine and coarse time steps respectively. Clearly $N_{\text{coarse}}$ has to divide $N_{\text{fine}}$. If $\delta t$ is the fine step length, then the coarse step length $\Delta t$ can be calculated as

$$ (5.5) \qquad \Delta t = \frac{\Delta T}{N_{\text{coarse}}} = \frac{\delta t \times N_{\text{fine}}}{N_{\text{coarse}}} $$

Also, $T = N_{\text{coarse}}\Delta tN = N_{\text{fine}}\delta tN$. The fine and coarse propagators are defined as

$$ (5.6) \qquad \mathbf{F} = \left\{ \left[ \mathbf{I} - \theta \delta t \mathbf{C} \right]^{-1} \left[ \mathbf{I} + (1 - \theta) \delta t \mathbf{C} \right] \right\}^{N_{\text{fine}}} \\ \mathbf{G} = \left\{ \left[ \mathbf{I} - \theta \Delta t \mathbf{C} \right]^{-1} \left[ \mathbf{I} + (1 - \theta) \Delta t \mathbf{C} \right] \right\}^{N_{\text{coarse}}} $$

In our experiments the fine time step $\delta t$ is defined based on the criterion:

$$ (5.7) \qquad \delta t \max |\lambda_{\mathbf{C}}| \le \omega_{\max} $$

where $\lambda_{\mathbf{C}}$ is an eigenvalue of $\mathbf{C}$ and $\omega_{\max}$ is some constant. The parameter values used for Parareal are put in Table 2.

<table><caption>TABLE 2
Parareal parameters for the linearised 1D shallow water model</caption><thead><tr><td>Parareal parameters</td><td>values</td></tr></thead><tbody><tr><td>Number of time windows, N</td><td>20</td></tr><tr><td>Theta parameter, θ</td><td>0.51</td></tr><tr><td>Number of fine time steps per time window, N<sub>fine</sub></td><td>100</td></tr><tr><td>Number of coarse time steps per time window, N<sub>coarse</sub></td><td>20</td></tr><tr><td>Fine time step, δt</td><td>0.05s</td></tr><tr><td>Coarse time step, Δt</td><td>0.25s</td></tr><tr><td>Courant number, ω<sub>max</sub></td><td>0.29</td></tr><tr><td>Total time period of integration, T</td><td>100s</td></tr></tbody></table>

Since an explicit diffusion term ($\mu \frac{\partial^2 u}{\partial x^2}$) is included in our 1D model, it is important to examine how varying the diffusion constant $\mu$ affects Parareal's convergence. Figure 2 shows the impact of $\mu$ on the total Parareal iterations required to solve system (5.2). Note that $\mu$ must be chosen such that

$$ (5.8) \qquad \mu \frac{\delta t}{\delta x} \le \frac{1}{2} $$

where the quantity on the left is the parabolic Courant number. We find that Parareal requires 15 iterations when no diffusion is applied ($\mu = 0$), compared to 6 Parareal iterations when $\mu = 1$. Figure 3 further shows the impact of diffusion values on the solution at the end of integration time. As expected, increasing $\mu$ reduces the Parareal iterations but also strongly damps the solution. For our experiments, $\mu = 0.15$ strikes a balance between reducing Parareal iterations and introducing only a slight damping (dotted line in Figure 3). This choice corresponds to the parabolic Courant number equal to 0.0075.

5.1.3 **Data assimilation cost function** Having defined $\mathbf{x} = (\eta, \mathbf{u})^T$ and $\mathbf{F}$, we can finally write our data assimilation cost function as

$$ (5.9) \qquad \min J(\mathbf{x}_0) = \frac{1}{2} \|\mathbf{F}^N \mathbf{x}_0 - \mathbf{y}\|_2^2 + \frac{\alpha}{2} \left\| \frac{\mathrm{d}\mathbf{x}_0}{\mathrm{d}x} \right\|_2^2 $$
**Algorithm 4.1** $p_*$ approximation test

1: **Given:** icg-iteration $j$

2: Find $\xi_j$ from equation (4.6) {Last parareal test}

3: **if** $j = 0$ **then**

4: Compute the parareal solution $P(k)\mathbf{p}_0$ with $\epsilon_p = \epsilon_{cg}/10$. Let the total iterations taken be $k_0$

5: Set $k_{sol} = k_0$, $\mathbf{p}_{**} = P(k)\mathbf{p}_0(k_{sol})$

6: **else**

7: Find $0 \le k < k_{sol}$ such that $\|P(k)\mathbf{p}_{j-1}(k) - \mathbf{p}_{**}\|_2 \le \xi_j$; save LHS as $\hat{\xi}_j$

8: $k_{sol} = k$

9: **if** $k = k_{sol} - 1$ **then**

10: cp = True {Use current parareal iterate}

11: **end if**

12: **end if**

13: Compute the current parareal solution $p_j$ for $k_{sol}$ iterations

14: Set $\mathbf{p}_{**} = P(k)\mathbf{p}_j(k_{sol})$

15: **if** cp = True **then**

16: Find $0 \le k < k_{sol}$ such that $\|P(k)\mathbf{p}_j - P(k)\mathbf{p}_j(k_{sol})\|_2 \le \xi_j$; save LHS as $\hat{\xi}_j$

17: **if** $k = k_{sol} - 1$ **then**

18: Perform one more parareal iteration and set $k_{sol} = k_{sol} + 1$

19: **end if**

20: Repeat steps in lines 15-17 till $k \neq k_{sol} - 1$

21: $k_{sol} = k$

22: **end if**

23: Calculate the inexact matrix-vector multiplication $\mathbf{c}_j = (F^N)^T \mathbf{p}_j(k_{sol})$

**Algorithm 4.2** Modified inexact conjugate gradient

1: **Given:** symmetric positive definite matrix $A \in \mathbb{R}^{n \times n}$, $\mathbf{b} \in \mathbb{R}^n$, $\epsilon > 0$, reorth, $\mathbf{x}_0 = 0$, $r_0 = -\mathbf{b}$, $\mathbf{p}_0 = \mathbf{b}$, $\beta_0 = \|\mathbf{b}\|_2^2$, $\mathbf{u}_1 = \mathbf{b}/\beta_0$, $\phi_0 = j_{max}$ and $\Phi_0 = 1$

2: {Iterations}

3: $j = 0$

4: **while** $j \le j_{max}$ **do**

5: Use the $\mathbf{p}_*$ approximation test in Algorithm 3. Also store the value of $\hat{\xi}_j$.

6: Compute $\hat{\phi}_j$ from (4.7)

7: $\Phi_{j+1} = \Phi_j - \hat{\phi}_j^{-1}$

8: **if** $j < j_{max}$ **then**

9: $\phi_{j+1} = (j_{max} - j - 1)/\Phi_{j+1}$

10: **else**

11: $\phi_{j+1} = \phi_j$

12: **end if**

13: $\alpha_j = \beta_j / \mathbf{p}_j^T \mathbf{c}_j$

14: $\mathbf{x}_{j+1} = \mathbf{x}_j + \alpha_j \mathbf{p}_j$

15: $\mathbf{r}_{j+1} = \mathbf{r}_j + \alpha_j \mathbf{c}_j$

16: **if** $(J_{j+1-d} - J_{j+1}) \le \frac{1}{4}\epsilon|J_{j+1}|$ **then**

17: **break**

18: **end if**

19: **if** (reorth) **then**

20: **for** $i = 1, 2, \dots, j$ **do**

21: $\mathbf{r}_{j+1} = \mathbf{r}_{j+1} - (\mathbf{u}_i^T \mathbf{r}_{j+1})\mathbf{u}_i$

22: **end for**

23: $\beta_{j+1} = \mathbf{r}_{j+1}^T \mathbf{r}_{j+1}$

24: $\mathbf{u}_{j+1} = \mathbf{r}_{j+1} / \sqrt{\beta_{j+1}}$

25: **else**

26: $\beta_{j+1} = \mathbf{r}_{j+1}^T \mathbf{r}_{j+1}$

27: **end if**

28: $\mathbf{p}_{j+1} = -\mathbf{r}_{j+1} + (\beta_{j+1}/\beta_j)\mathbf{p}_j$

29: **end while**
<figure><img src="image_3.png" /><figcaption>FIG. 2. Total Parareal iterations when the diffusion constant $\mu$ varies in the 1d model. A sharp decline is seen in the Parareal iterations when $\mu$ increases before reaching a point where the Parareal iterations almost remain unchanged.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3. Effect of varying $\mu$ on the model solution at the end of integration time for free surface $\eta$ (left) and velocity $u$ (right). As expected, a significant damping is observed when a higher $\mu$ is used on the model.</figcaption></figure>

where we use a regularisation term with regularisation constant $\alpha$ to improve the conditioning of $A = (F^N)^T F^N$. Without regularisation and with the computational parameters given in Table 1 and Table 2, $\kappa(A) \approx 2.98 \times 10^{18}$ which is extremely large. The following parameters values are used for data assimilation as shown in Table 3.

<table><caption>TABLE 3
Data assimilation parameters for the linearised 1D shallow water model</caption><thead><tr><th>Data assimilation parameters</th><th>values</th></tr></thead><tbody><tr><td>Regularisation constant</td><td>$\alpha = 10^{-5}$</td></tr><tr><td>Initial free surface value</td><td>$\eta_0(x) = \exp\{[(x - L/2)/(L/15)]^2\}$</td></tr><tr><td>Initial velocity value</td><td>$u_0(x) = 0 \quad \forall x$</td></tr></tbody></table>

Remark 5.1. A standard regularisation term could have been $\alpha \|x_0\|$, but our pri-
iterate, $\xi_{17}$ lies between the second last and the last iteration $k$. So we proceed for the current parareal iterate (on the right) and now $\xi_{17}$ is above the last two $k$.

<figure><img src="image_6.png" /><figcaption>FIG. 4.3. Using last and current parareal iterate</figcaption></figure>

**5. Numerical experiments.** One of the main objectives of coupling two inherently sequential methods is to use it for the applications related to ocean or climate modelling. In general, such problems are solved at a low resolution by integrating the model for a very long period of time for example some decades or centuries. This gives us an opportunity to investigate the parallel-in-time methods for large grid sizes. We start with the simplest case of a linearised one-dimensional shallow water equations.

$$ (5.1) \qquad \begin{aligned} \frac{\partial \eta}{\partial t} &= -h \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial t} &= -g \frac{\partial \eta}{\partial x} \end{aligned} $$

where $\eta$ and $u$ are the free surface and velocity respectively. $h$ is the characteristic height and $g$ is gravity. This is clearly a hyperbolic problem and to numerically solve it we are going to add an explicit diffusion. The reason being to explain the effects of the unresolved scales in terms of the large scales. This so called the closure of the system is performed by the Reynolds averaging [5]. The modified equations are given as

$$ (5.2) \qquad \begin{aligned} \frac{\partial \eta}{\partial t} &= -h \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial t} &= -g \frac{\partial \eta}{\partial x} + \mu \frac{\partial^2 u}{\partial x^2} \end{aligned} $$

where $\mu$ is the diffusion constant. The numerical discretisation is done using Arakawa C-grid or a staggered grid.

$$ (5.3) \qquad \begin{aligned} \frac{\eta_{i+1/2}^{k+1} - \eta_{i+1/2}^{k}}{\Delta t} &= -h \frac{u_{i+1}^{k} - u_{i}^{k}}{\Delta x} \\ \frac{u_{i+1}^{k+1} - u_{i+1}^{k}}{\Delta t} &= g \frac{\eta_{i+1/2}^{k} - \eta_{i-1/2}^{k}}{\Delta x} + \frac{u_{i+1}^{k} - 2u_{i}^{k} + u_{i-1}^{k}}{\Delta x^2} \end{aligned} $$

The above discretisation can be expressed as an ODE system

$$ (5.4) \qquad \frac{dx}{dt} = Cx $$
mary interest is to explore the model dynamics using Parareal and how it impacts the minimisation without constraining the initial state. Instead we use a regularisation term based on the spatial second derivative of the initial state to improve the conditioning of the Hessian. It also ensures the smoothness or stability of the numerical experiments by controlling high-frequency variation or abrupt spatial or temporal changes in the solution.

The conditioning of the regularised cost function after using the regularisation constant given in Table 3 is $\kappa(\mathbf{A}) \simeq 10^6$. The gradient also changes to

$$ (5.10) \qquad \nabla J(\mathbf{x}_0) = (\mathbf{F}^N)^T (\mathbf{F}^N \mathbf{x}_0 - \mathbf{y}) + \alpha \mathbf{Q}_2 \mathbf{x}_0 $$

where $\mathbf{Q}_2 = \frac{1}{\Delta x^2} \mathbf{Q}_1^T \mathbf{Q}_1$ and

$$ (5.11) \qquad \mathbf{Q}_1 = \begin{pmatrix} \begin{array}{cccc} -1 & 1 & & \\ & \ddots & \ddots & \\ & & -1 & 1 \\ & & & -1 \end{array} & \mathbf{0} \\ \hline \mathbf{0} & \begin{array}{cccc} -1 & 1 & & \\ & \ddots & \ddots & \\ & & -1 & 1 \\ & & & -1 \end{array} \end{pmatrix}. $$

Each diagonal block matrix in $\mathbf{Q}_1$ is of size $N_x - 2 \times N_x - 2$ and so $\mathbf{Q}_1, \mathbf{Q}_2 \in \mathbb{R}^{2(N_x-2) \times 2(N_x-2)}$. Therefore, $\mathbf{A} = (\mathbf{F}^N)^T \mathbf{F}^N + \alpha \mathbf{Q}_2$.

In the subsequent subsections, we demonstrate how the performance of minimisation changes when transitioning from CG without Parareal to CG with Parareal, and finally to inexact CG with Parareal. We also show how the incorporation of feasible approximations and the adoption of a modified Parareal stopping criterion within inexact CG significantly enhance Parareal's performance. The results for inexact CG follow from Theorem 4.1. Specifically once the quadratic change $q(\mathbf{x}_j) - q(\mathbf{x}_*)$ is bounded by the required precision level (4.6), the corresponding $\mathbf{x}_j$ represents the retrieved optimal initial state.

## 5.2 Conjugate gradient (CG)
We begin by examining the performance of CG using the 1D shallow water model with exact matrix-vector multiplication is used. Note that this corresponds to using CG without Parareal. For a given CG tolerance $\epsilon_{cg}$, the minimisation ends when

$$ (5.12) \qquad \|\mathbf{r}(\mathbf{x}_j)\|_2 \le \epsilon_{cg}. $$

Figure 4 below shows the evolution of the residual $\|\mathbf{r}(\mathbf{x}_j)\|_2$ (solid line) with respect to the CG iterations for $\epsilon_{cg} = 10^{-4}$ (dashed line). It can be seen that the minimisation terminates after 24 iterations.

## 5.3 Conjugate gradient with Parareal
Our next step is to evaluate the performance of CG when Parareal is introduced. In this case, the matrix-vector products are computed using Parareal, with their accuracy depending on the fixed Parareal stopping tolerance $\epsilon_p$. It is important to note that the exact Parareal precision required to match the performance of CG without Parareal (exact matrix-vector product) is not known beforehand.

To determine the optimal $\epsilon_p$, a hit-and-trial strategy is followed by running CG with arbitrary Parareal tolerances. For simplicity, we choose $\epsilon_p = 10^{-6}$, which corresponds to the minimum of the $\{\xi_j\}_j$ values produced by inexact CG (as discussed in
where

$$ C = \begin{pmatrix} 0 & \begin{array}{c} -h/\Delta x \\ h/\Delta x \end{array} & \begin{array}{c} -h/\Delta x \\ h/\Delta x \end{array} & \cdots & \cdots & h/\Delta x & -h/\Delta x \\ \hline g/\Delta x & -g/\Delta x & \mu/\Delta x^2 & \mu/\Delta x^2 & \cdots & \mu/\Delta x^2 \\ g/\Delta x & -g/\Delta x & -2\mu/\Delta x^2 & -2\mu/\Delta x^2 & \cdots & -2\mu/\Delta x^2 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ g/\Delta x & -g/\Delta x & \mu/\Delta x^2 & \mu/\Delta x^2 & \cdots & \mu/\Delta x^2 \end{array} \end{pmatrix} $$

and $\mathbf{x} = (\eta_0, \eta_1, \dots, \eta_{59}, u_0, u_1, \dots, u_{59})^T$. We take $h = 0.9$, $g = 10$ and the characteristic length $L = 120$. If the number of spatial grid points are 120, then $\Delta x = L/120 = 1$. To solve the system (5.4) we use a weighted $\theta$-scheme

$$ \frac{\mathbf{x}^{k+1} - \mathbf{x}^k}{\Delta t} = C[\theta\mathbf{x}^{k+1} + (1-\theta)\mathbf{x}^k] $$

which we can rewrite as

$$ \mathbf{x}^{k+1} = \frac{[I + C(1 - \theta)\Delta t]}{I - C\theta\Delta t} \mathbf{x}^k $$

We can now define the propagators for implementing the parareal algorithm. For each time window, we associate $N_f$ and $N_g$ as the number of fine time steps and coarse time steps respectively. Clearly $N_g$ has to divide $N_f$. If $\delta t$ is the fine step length, then the coarse step length $\Delta t$ is given as

$$ \Delta t = \frac{\delta t \times N_f}{N_g} $$

while the total integration time is $T = \delta t \times N_f \times N$. Here, $N = 20$, $N_f = 100$ and $N_g = 25$. The fine and coarse propagators are defined as

$$ F = \left( \frac{I + (1 - \theta)\delta t C}{I - \theta \delta t C} \right)^{N_f} , \quad G = \left( \frac{I + (1 - \theta)\Delta t C}{I - \theta \Delta t C} \right)^{N_g} $$

With our model in hand, we solve the following data assimilation problem

$$ \min J(\mathbf{x}_0) = \frac{1}{2} \| F^N \mathbf{x}_0 - \mathbf{y} \|_2^2 + \frac{\alpha}{2} \left\| \frac{d\mathbf{x}_0}{d\mathbf{x}} \right\|_2^2 $$

The following values of the parameters are used:

* $\theta = 0.51$

* regularisation parameter $\alpha = 10^{-5}$

* viscosity coefficient $\mu = 0.15$

* fine time step $\delta t = 0.05$

* initial condition: a Gaussian at the middle of the domain $\mathbf{x}_0 = e^{[(x-L/2)/(L/15)]^2}$

**5.1. Analysis.** First we check how our data assimilation problem performs with our shallow water model when we run the exact CG method with exact matrix-vector multiplication. The stopping criterion depends on the 2-norm of the residual. For a tolerance $\epsilon_{cg}$ of $10^{-4}$, it takes 24 minimisation iterations.
<figure><img src="image_5.png" /><figcaption>FIG. 4. CG when matrix-vector products are exact and Parareal is not used. Minimisation terminates after 24 iterations when the residual norm $||r(x_j)||_2$ (solid) becomes less than the chosen tolerance $\epsilon_{cg} = 10^{-4}$ (dashed).</figcaption></figure>

subsection 5.3.3). It is important to emphasise note that this choice is used solely for comparison purposes and is not practically available, as it requires a prior run of in-exact CG. Figure 5 shows the result of minimisation when $\epsilon_{cg} = 10^{-4}$ and $\epsilon_p = 10^{-6}$. The minimisation ends after 24 CG iterations with a total of 191 Parareal iterations. This corresponds to an average of 7.96 Parareal iterations per CG iteration, which will

<figure><img src="image_6.png" /><figcaption>FIG. 5. CG when matrix-vector products come from Parareal prescribed a fixed Parareal tolerance $\epsilon_p$. Left image shows identical residual norm values as in Figure 4 when CG tolerance $\epsilon_{cg} = 10^{-4}$ (dashed) and $\epsilon_p = 10^{-6}$ are used. Right image shows the evolution of Parareal iterations with CG iterations.</figcaption></figure>

be shown to be slightly higher than the average of 6.25 Parareal iterations (as in subsection 5.3.3) when using inexact CG. This result highlights the need for an *adaptive* Parareal within the inexact CG framework, demonstrating its advantage over fixed $\epsilon_p$ obtained through hit-and-trial.

To further assess CG performance in terms of the relative error in the 2-norm of
<figure><img src="image_7.png" /><figcaption>FIG. 5.1. Exact CG with exact matrix-vector product</figcaption></figure>

From the above exact CG solution with $\epsilon_{cg} = 10^{-4}$ in Figure 5.1, we find the value of $\|\mathbf{r}_k\|_{A^{-1}}$ and use it in (4.5) to calculate a suitable stopping criteria for the inexact CG. We find that $\epsilon_{icg} = 1.12 \times 10^{-7}$. As talked earlier the residuals generated are no longer exact and for that purpose we have added a reorthogonalisation step at each icg-iteration.

We first show the result when we use all theoretical estimates as well as the norm $\|E_j\|_{A^{-1},A}$ and why it is not a good idea to keep the primal-dual norm as the stopping criteria. While we are still able to converge after 24 minimisation iterations, but only at the expense of 464 p-iterations. This makes it an average of 18.5 p-iterations per cg-iteration which is not feasible at all. Replacing $\|E_j\|_{A^{-1},A}$ with $\|E_j \mathbf{p}_j\|_{A^{-1},A}$ as the stopping criterion, it takes 25 minimisation iterations but the total p-iterations goes down to 165 and we have an average of 6.6 p-iterations per cg-iteration.

<figure><img src="image_8.png" /><figcaption>FIG. 5.2. Theoretical inexact CG with $\|E_j\|_{A^{-1},A}$ as the stopping criteria</figcaption></figure>
the residual and the total number of Parareal iterations, Figure 6 presents the results for CG with $\epsilon_p = 10^{-5}, 10^{-6}$ and $10^{-7}$ are used. The relative error in the 2-norm of the residual is defined as

$$ \text{Relative error in the 2-norm of residual} = \frac{\|\mathbf{r}_j\|_2 - \|\mathbf{r}(\mathbf{x}_j)\|_2}{\|\mathbf{r}(\mathbf{x}_j)\|_2} $$

Initially, the relative error is low and nearly identical for all three cases. However, after a few CG iterations it begins to vary, becoming smaller for smaller $\epsilon_p$. The relative error also gradually increases with each CG iteration, suggesting that a lower Parareal precision is needed near the end of the minimisation. As expected, a smaller $\epsilon_p$ implies more Parareal iterations and vice-versa.

<figure><img src="image_7.png" /><figcaption>FIG. 6. Relative error in the 2-norm of the residual (left) and the corresponding Parareal iterations k (right) for fixed values of $\epsilon_p = 10^{-5}$ (dashed), $\epsilon_p = 10^{-6}$ (dotted), and $\epsilon_p = 10^{-7}$ (dashdot) when using CG.</figcaption></figure>

### 5.3.1 Inexact conjugate gradient using Parareal
We now present the results of our experiments using inexact CG, with the CG results serving as a reference for comparison. As noted earlier, CG and inexact CG use different stopping criteria: CG uses the 2-norm of the residual, while inexact CG relies on the $A^{-1}$ norm (see Theorem 4.1). To ensure both methods follow the same stopping criterion we see from (4.5) that inexact CG terminates when, for a given value of $\epsilon_{icg}$,

$$ (5.13) \qquad \|\mathbf{r}_j\|_{A^{-1}} \le \frac{\sqrt{\epsilon_{icg}}}{2} \|\mathbf{b}\|_{A^{-1}}. $$

Rearranging the terms for equality, we find

$$ (5.14) \qquad \epsilon_{icg} = \left( \frac{2 \|\mathbf{r}_j\|_{A^{-1}}}{\|\mathbf{b}\|_{A^{-1}}} \right)^2. $$

To determine the value of $\epsilon_{icg}$ corresponding to the CG stopping tolerance, we use the value of $\|\mathbf{r}_j\|_{A^{-1}}$ at convergence (last CG iteration) (from Figure 5) and deduce $\epsilon_{icg}$ using equation (5.14). In this case, $\epsilon_{icg} = 1.12 \times 10^{-7}$.

We first show the results of inexact CG when the original Parareal stopping criterion ($\|\mathbf{E}_j\|_{A^{-1},A}$ bounded by $\omega_j$ in (4.3)) is used in Figure 7.
<figure><img src="image_9.png" /><figcaption>FIG. 5.3. Theoretical inexact CG with $||E_j p_j||_{A-1}$ as the stopping criteria</figcaption></figure>

With the approximations of the quantities from Section 4.3 and the $||E_j p_j||_{A-1}$, the practical inexact CG performs equally well. The integer value taken for the stopping criterion is $d = 2$. The minimisation ends in 26 iterations with a slightly lower average of 6.2 p-iterations per cg-iteration.

<figure><img src="image_10.png" /><figcaption>FIG. 5.4. Practical inexact CG with $||E_j p_j||_{A-1}$ as the stopping criteria</figcaption></figure>

Finally, we show the main result by combining the practical estimates and the approximation for the modified stopping criterion. We choose $d = 10$. Still, the method gives almost identical results as seen in the last two cases with the minimisation and p-iterations being 25 and 159 respectively. The average p-iterations per CG iterations remains at 6.36. The p-iterations oscillate by $\pm 1$ which can be expected because we assume that the value of $\xi_j$ should not lie between the last two parareal iterates.
<figure><img src="image_8.png" /><figcaption>FIG. 7. Inexact CG following Gratton et al. with original Parareal stopping criterion $||E_j||_{A^{-1},A} \le \omega_j$. The left image shows the residual norm $||r_j||_{A^{-1}}$ (dashdotted), inexact CG stopping tolerance $1/2\sqrt{\epsilon_{icg}}||b||_{A^{-1}}$ (dotted), bound for the error norm $\omega_j$ (densely dashdotted), residual gap norm $||r(x_j) - r_j||_{A^{-1}}$ (loosely dotted), and the cost function change $|J(x_j) - J(x_*)|$ (densely dashed) which is bounded by the quantity $\epsilon_{icg}|J(x_*)|$ (dashed). The right image shows the evolution of corresponding Parareal iterations k with icg-iterations j. No approximations are used.</figcaption></figure>

**Remark 5.2.** Throughout the manuscript we will adhere to the same line style coding for quantities used in Figure 7 across different inexact CG configurations.

As a reminder from Theorem 4.1 inexact CG terminates when both

$$ (5.15) \qquad ||\mathbf{r}_j||_{A^{-1}} \le \frac{\sqrt{\epsilon_{icg}}}{2} ||\mathbf{b}||_{A^{-1}} \quad \text{and} \quad ||\mathbf{r}(\mathbf{x}_j) - \mathbf{r}_j||_{A^{-1}} \le \frac{\sqrt{\epsilon_{icg}}}{2} ||\mathbf{b}||_{A^{-1}} $$

and we have $|J(\mathbf{x}_j) - J(\mathbf{x}_*)| \le \epsilon_{icg} |J(\mathbf{x}_*)|$. From Figure 7 that at the end of the minimisation the conditions of Theorem 4.1 are well satisfied. The minimisation ends after 24 icg-iterations involving a total of 436 Parareal iterations. On average, 18.2 Parareal iterations per icg-iteration are used which is almost equal to the total number of time windows ($N = 20$). This implies negligible speedup since convergence is achieved only during the second last or last iteration. The use of $||E_j||_{A^{-1},A}$ for deriving the Parareal stopping criterion leads to a suboptimal algorithm.

**5.3.2 Illustration of the modified Parareal criterion** Here we demonstrate the significant reduction in the number of Parareal iterations when switching from the original to the modified Parareal stopping criterion for inexact CG. Figure 8 below compares the two criteria for a particular icg-iteration j. The left image on uses the original criterion (involving $||E_j||_{A^{-1},A}$), while the right image uses the modified criterion (involving $||E_j p_j||_{A^{-1}}$). A substantial difference is evident: only 6 Parareal iterations are required with the modified criterion compared to 19 Parareal iterations with the original criterion. If we look at the definition,

$$ (5.16) \qquad ||\mathbf{E}_j||_{A^{-1},A} = \sup_{\delta \mathbf{x} \neq 0} \frac{||\mathbf{E}_j \delta \mathbf{x}||_{A^{-1}}}{||\delta \mathbf{x}||_{A}} $$

we observe that $||E_j||_{A^{-1},A}$ does not involve the conjugate direction vector $p_j$ as a product with $E_j$. Theoretically, as j increases we would expect $||p_j||_A \to 0$, which should lead to decreasing $||E_j p_j||_{A^{-1}}$ values. However, since $p_j$ is absent in the original criterion, the relatively high values present the most pessimistic case.
<figure><img src="image_11.png" /><figcaption>FIG. 5.5. Modified practical inexact CG</figcaption></figure>

Looking at all the previous results, we observe that if we take the minimum of all the values of $\xi_j$ which turns out to be approximately $10^{-6}$, running the exact CG with now the matrix-vector multiplication coming through the parareal method, it should work and it does. For $\epsilon_p = 10^{-6}$ as a fixed parareal stopping criterion, we end up using 193 p-iterations. However on average we use 7.72 p-iterations per cg-iteration which is slightly higher due to a fixed parareal tolerance.

<figure><img src="image_12.png" /><figcaption>FIG. 5.6. Exact CG with matrix-vector product from parareal</figcaption></figure>

Some of the obtained parameter values

* courant number: 0.55

* diffusion courant number: 0.0075
<figure><img src="image_9.png" /><figcaption>FIG. 8. Comparison of original (left) and modified (right) Parareal stopping criterion. The norm values (solid) satisfy the corresponding tolerances $\omega_j$ (dashed) and $\xi_j$ (dotted) for a given icg-iteration $j$. A big difference can be seen with the original criterion taking a lot more Parareal iterations than with the modified criterion.</figcaption></figure>

The benefits of using the modified Parareal stopping criterion for inexact CG are shown in Figure 9. While the minimisation again ends in 24 iterations, the total number of Parareal iterations is significantly reduced to 132, compared to 436 Parareal iterations in subsection 5.3.1). This corresponds to an average of 5.5 Parareal iterations per icg-iteration, a big improvement.

<figure><img src="image_10.png" /><figcaption>FIG. 9. Inexact CG following Gratton et al. but now with modified Parareal stopping criterion $||\mathbf{E}_j\mathbf{p}_j||_{\mathbf{A}^{-1}} \le \xi_j$. The new quantity $\xi_j$ (dashdotdotted) on the left image is the bound for $||\mathbf{E}_j\mathbf{p}_j||_{\mathbf{A}^{-1}}$. Still no approximations are used. The number of Parareal iterations $k$ required at a given icg-iteration $j$ is plotted on the right image.</figcaption></figure>

Note that we no longer have to fix the Parareal tolerance $\epsilon_p$ a priori. The modified Parareal criterion allows inexact CG to use Parareal adaptively since its tolerance ($\xi_j$) adjusts with every icg-iteration. Additionally, Figure 9 includes the quantity $\omega_j$ to show that while the permissible matrix-vector product error ($||\mathbf{E}_j||_{\mathbf{A}^{-1},\mathbf{A}}$) gradually increases with icg-iterations, then method still converges and maintains the same level of accuracy as CG.
• expected speedup:

$$ \psi = \frac{\text{number of time windows}}{\text{average p-iterations per cg-iteration}} = \frac{20}{6.36} = 3.145 $$

The initial condition retrieved by 4D-Var by using exact and inexact CG is almost identical to the initial state $x_0$

<figure><img src="image_13.png" /></figure>

As one would expect, the free surface and velocity values also coincide with the values obtained from the reference solution.

<figure><img src="image_14.png" /></figure>

**6. Conclusion and future work.** In an attempt to couple a parallel-in-time method (parareal) with 4D-Var we see that if we use the proposed modified inexact CG method one would expect the optimum use of p-iterations. Even though the p-iterations do not decrease monotonically which could be accounted for by the values of $\xi_j$, we still get an expected speedup of 3. We can at least conclude that the inexact CG by Gratton et al. [13] provides a well defined stopping criterion for the parareal algorithm at each minimisation iteration. This way we don't have to run parareal with a very fine precision each
5.3.3 Inexact conjugate gradient with practical estimates In this sub-section we are going to use the derived estimate of $||\mathbf{E}_j\mathbf{p}_j||_{\mathbf{A}^{-1}}$ for the modified Parareal stopping criterion. Recall that the approximation is given as

$$ (5.17) \qquad \|\mathbf{E}_j\mathbf{p}_j\|_{\mathbf{A}^{-1}} \approx \|\mathbf{P}(k+1)\mathbf{p}_j - \mathbf{P}(k)\mathbf{p}_j\|_2 $$

To assess the accuracy of our approximation, Figure 10 compares the exact and approximate values of $||\mathbf{E}_j\mathbf{p}_j||_{\mathbf{A}^{-1}}$ for some icg-iterations. It can be seen that both exact and approximate values align, providing a nice approximation of $||\mathbf{E}_j\mathbf{p}_j||_{\mathbf{A}^{-1}}$ by the end of the Parareal iteration (where the Parareal stopping criterion is satisfied by $\xi_j$). Figure 11 further illustrates the results of using approximate $||\mathbf{E}_j\mathbf{p}_j||_{\mathbf{A}^{-1}}$ in the mod-

<figure><img src="image_11.png" /><figcaption>FIG. 10. Comparison of the exact (solid) and approximate (dashed) values of $||\mathbf{E}_j\mathbf{p}_j||_{\mathbf{A}^{-1}}$ for icg-iterations 1, 2, 10, 15, 20, 24. The relatively close values at each Parareal iteration suggest that the approximation is a good one.</figcaption></figure>

ified Parareal criterion for inexact CG. The minimisation performs 24 icg-iterations and a total of 153 Parareal iterations, resulting in an average of 6.375 Parareal iterations per icg-iterations. If we compare Figure 11 and Figure 9 (result with exact $||\mathbf{E}_j\mathbf{p}_j||_{\mathbf{A}^{-1}}$), minimal differences can be observed on the left images among various quantities. The only difference is noticed on the right images where inexact CG with approximated $||\mathbf{E}_j\mathbf{p}_j||_{\mathbf{A}^{-1}}$ requires one more Parareal iteration per icg-iteration due to the construction of our approximation (the use of $\mathbf{P}(k+1)$).

Finally we present the main result of our work in Figure 12 which uses practical estimates for all the remaining norms (from subsection 4.3), approximate $||\mathbf{E}_j\mathbf{p}_j||_{\mathbf{A}^{-1}}$, and the inaccuracy budget within inexact CG. This result follows from Algorithm 4.2 which also requires value of an integer $d$ ($d=2$ in this case) used for the inexact CG termination criterion (4.5). Even after using all approximations, we see that inexact CG gives similar results to the other variants of inexact CG discussed above. The minimisation ends after 24 iterations and needing a total of 154 Parareal iterations. On average, the number of Parareal iterations per icg-iteration remains 6.4. The maximum possible speedup is bounded by a factor of

$$ (5.18) \qquad S = \frac{\text{number of time windows}}{\text{average Parareal iterations per icg-iteration}} = \frac{20}{6.4} \approx 3.12 $$
time but only till where the bound is satisfied.

As a remark, it is a well known fact that the parareal algorithm shows instabilities when it is applied to hyperbolic problems and second order ODEs [25, 28]. One way to considerably increase the parareal's performance is by constructing a coarse solver from the subspace of the initial conditions at each time sub-domain [8]. This parareal variant is known as the Krylov enhanced parareal method. As a side experiment, we ran our modified inexact CG method with the Krylov enhanced version and for 25 minimisation iterations the average p-iterations per cg-iteration almost get halved (3.24).

Clearly our task becomes much simpler. However, at each parareal iteration one has to store all the fine evaluations for constructing the subspace leading to computational inefficiency. Here it works well because the problem is a simple 1-D case but it would be really difficult to implement for higher dimensional problems. Choosing appropriate basis vectors which contribute to the subspace could be one of the solutions which can be seen as a topic of future study.

<figure><img src="image_15.png" /><figcaption>FIG. 6.1. Modified practical inexact CG with krylov enhanced parareal</figcaption></figure>

One thing to note is that the adjoint in the cost function (2.7) still uses the transpose of the matrix. We have not introduced time parallelisation for the backward model integration which could be the next step. This could be applied in exactly similar manner and can double the degree of parallelism in our algorithm. Finally parareal is just one of the techniques which can be employed for data assimilation. There are several other parallel-in-time methods for example the PFASST algorithm [6] or the ParaOpt algorithm [11] where the forward and backward model are solved together as a single optimality system.

We are currently investigating the idea of introducing asynchronous parallel-in-time algorithms which makes the processors independent of one another and therefore we can cut the communication time and obtain further speedup. The idea would be to make these p-iterations asynchronous [18, 19] both for the forward and the backward model integration.

## REFERENCES

[1] L. BAFFICO, S. BERNARD, Y. MADAY, G. TURINICI, AND G. ZÉRAH, *Parallel-in-time molecular-dynamics simulations*, Physical Review E, 66 (2002), p. 057701.

[2] G. BAL, *Parallelization in time of (stochastic) ordinary differential equations*, Math. Meth. Anal. Num. (submitted), (2003).

[3] F. BOUTTIER AND P. COURTIER, *Data assimilation concepts and methods march 1999*, Meteorological training course lecture series. ECMWF, 718 (2002), p. 59.

[4] P. COURTIER, J.-N. THÉPAUT, AND A. HOLLINGSWORTH, *A strategy for operational implementation of 4d-var, using an incremental approach*, Quarterly Journal of the Royal Meteorological Society, 120 (1994), pp. 1367–1387.

[5] B. CUSHMAN-ROISIN AND J.-M. BECKERS, *Introduction to geophysical fluid dynamics: physical and numerical aspects*, Academic press, 2011.
<figure><img src="image_12.png" /><figcaption>FIG. 11. Inexact CG using modified Parareal stopping criterion but now with approximation for $||E_j p_j||_{A^{-1}}$ obtained in subsection 4.4. The minimisation result and the corresponding required Parareal iterations $k$ are shown on left and right images respectively. We end up using one more Parareal iteration per icg-iteration due to the approximation used.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 12. Inexact CG using: approximation for $||E_j p_j||_{A^{-1}}$, practical estimates from subsection 4.3, and the inaccuracy budget from subsection 4.5. The results are almost identical to Figure 11 showing that even after all the approximations the algorithm performs well.</figcaption></figure>

To check the accuracy of the obtained initial state from inexact CG, we compare the contour plots of the solution from inexact CG with the exact solution for the whole time domain as shown in Figure 13. By inexact CG solution and exact solution, we mean the solutions obtained by applying the fine solver F to the initial state retrieved from inexact CG and to the given initial state respectively. The free surface values $\eta$ and velocity values $u$ are plotted separately on the top and bottom row of Figure 13 respectively. Visual inspection of Figure 13 shows minimal differences, indicating that the retrieved initial state is sufficiently accurate. Additionally, Figure 14 presents a contour plot of the error between the exact solution and the inexact CG solution, confirming that the errors are indeed small.
[6] M. EMMETT AND M. MINION, *Toward an efficient parallel in time method for partial differential equations*, Communications in Applied Mathematics and Computational Science, 7 (2012), pp. 105–132.

[7] M. FISHER AND S. GÜROL, *Parallelization in the time dimension of four-dimensional variational data assimilation*, Quarterly Journal of the Royal Meteorological Society, 143 (2017), pp. 1136–1147.

[8] M. GANDER AND M. PETCU, *Analysis of a krylov subspace enhanced parareal algorithm for linear problems*, in ESAIM: Proceedings, vol. 25, EDP Sciences, 2008, pp. 114–129.

[9] M. J. GANDER, *50 years of time parallel time integration*, in Multiple shooting and time domain decomposition methods, Springer, 2015, pp. 69–113.

[10] M. J. GANDER AND E. HAIRER, *Nonlinear convergence analysis for the parareal algorithm*, in Domain decomposition methods in science and engineering XVII, Springer, 2008, pp. 45–56.

[11] M. J. GANDER, F. KWOK, AND J. SALOMON, *Paraopt: A parareal algorithm for optimality systems*, SIAM Journal on Scientific Computing, 42 (2020), pp. A2773–A2802.

[12] M. J. GANDER AND S. VANDEWALLE, *Analysis of the parareal time-parallel time-integration method*, SIAM Journal on Scientific Computing, 29 (2007), pp. 556–578.

[13] S. GRATTON, E. SIMON, D. TITLEY-PELOQUIN, AND P. L. TOINT, *Minimizing convex quadratics with variable precision conjugate gradients*, Numerical Linear Algebra with Applications, 28 (2021), p. e2337.

[14] I. C. IPSEN AND C. D. MEYER, *The idea behind krylov methods*, The American mathematical monthly, 105 (1998), pp. 889–899.

[15] A. LAWLESS, S. GRATTON, AND N. NICHOLS, *An investigation of incremental 4d-var using non-tangent linear models*, Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography, 131 (2005), pp. 459–476.

[16] F.-X. LE DIMET AND O. TALAGRAND, *Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects*, Tellus A: Dynamic Meteorology and Oceanography, 38 (1986), pp. 97–110.

[17] J.-L. LIONS, Y. MADAY, AND G. TURINICI, *Résolution d'edp par un schéma en temps pararéel*, Comptes Rendus de l'Académie des Sciences-Series I-Mathematics, 332 (2001), pp. 661–668.

[18] F. MAGOULES AND G. GBIKPI-BENISSAN, *Asynchronous parareal time discretization for partial differential equations*, SIAM Journal on Scientific Computing, 40 (2018), pp. C704–C725.

[19] F. MAGOULÈS, G. GBIKPI-BENISSAN, AND Q. ZOU, *Asynchronous iterations of parareal algorithm for option pricing models*, Mathematics, 6 (2018), p. 45.

[20] M. MINION, *A hybrid parareal spectral deferred corrections method*, Communications in Applied Mathematics and Computational Science, 5 (2011), pp. 265–301.

[21] J. NIEVERGELT, *Parallel methods for integrating ordinary differential equations*, Communications of the ACM, 7 (1964), pp. 731–733.

[22] B. W. ONG AND J. B. SCHRODER, *Applications of time parallelization*, Computing and Visualization in Science, 23 (2020), pp. 1–15.

[23] J. RANTAKOKKO, *Strategies for parallel variational data assimilation*, Parallel Computing, 23 (1997), pp. 2017–2039.

[24] V. RAO AND A. SANDU, *A time-parallel approach to strong-constraint four-dimensional variational data assimilation*, Journal of Computational Physics, 313 (2016), pp. 583–593.

[25] D. RUPRECHT, *Wave propagation characteristics of parareal*, Computing and Visualization in Science, 19 (2018), pp. 1–17.

[26] D. RUPRECHT AND R. KRAUSE, *Explicit parallel-in-time integration of a linear acoustic-advection system*, Computers & Fluids, 59 (2012), pp. 72–83.

[27] V. SIMONCINI AND D. B. SZYLD, *Theory of inexact krylov subspace methods and applications to scientific computing*, SIAM Journal on Scientific Computing, 25 (2003), pp. 454–477.

[28] G. A. STAFF AND E. M. RONQUIST, *Stability of the parareal algorithm*, in Domain decomposition methods in science and engineering, Springer, 2005, pp. 449–456.

[29] Y. TRÉMOLET AND F.-X. LE DIMET, *Parallel algorithms for variational data assimilation and coupling models*, Parallel Computing, 22 (1996), pp. 657–674.

[30] J. VAN DEN ESHOF AND G. L. SLEIJPEN, *Inexact krylov subspace methods for linear systems*, SIAM Journal on Matrix Analysis and Applications, 26 (2004), pp. 125–153.
<figure><img src="image_14.png" /><figcaption>FIG. 13. Contour plots of the iterative solution (left column) and exact solution (right column) for free surface (first row) and velocity (second row). The iterative and exact solution are obtained by running the fine solver on the initial state retrieved by inexact CG and on the given (exact) initial state respectively. A first glance suggests that the retrieved initial state is very close to the exact initial state.</figcaption></figure>

**6 Conclusion and Future Work** In this paper we focused on time parallelising the inner loop of the incremental 4D-Var by utilising Parareal for the tangent-linear model integration. This led to solving an approximate linear system using an inexact version of the conjugate gradient method. The advantage of inexact CG over the usual CG is that it allowed us to adaptively control the required Parareal's accuracy. We also proposed a modified precision criterion based on the original criterion which significantly reduced the required number of required Parareal iterations. Numerical results on the 1D linearised shallow water equations showed that inexact CG with all approximations outperformed CG with Parareal (with a fixed Parareal tolerance) by giving a speedup factor of 3.12.

Areas of future work include extending these developments to a more realistic 2D shallow water model. In that case the inclusion of parameters such as Coriolis force and additional dimension will make the problem more challenging in terms of size and complexity. Another point of interest is the introduction of potential time parallelism in the adjoint direction, i.e. time parallelism for the backward model ($F^N)^T$ in equation (3.5). This aspect has not been explored in our work since the existing approximation for $\|E_j p_j\|_{A^{-1}}$ is only valid when exact adjoint is used. If the adjoint could be parallelised, a natural extension of this idea would involve using asynchronous Parareal iterations [31, 30] where the algorithm runs without a synchronisation
<figure><img src="image_15.png" /><figcaption>FIG. 14. Contour plots depicting the error between the inexact CG solution and exact solution for free surface (left) and velocity (right).</figcaption></figure>

point. In that case the forward and adjoint parallel runs could be done simultaneously without needing to wait for any exchange of information.

Finally, some open questions remain which are worth investigating in the future. One question is whether the modified Parareal stopping criterion derived here is also suitable for other types of PDEs or data assimilation setup. In the numerical experiments, the residual mismatch norm $\|r(x_j) - r_j\|_{A^{-1}}$ in Figure 9 is noticeably larger than other examples. This may explain the better Parareal performance, as less computational effort is spent on making the mismatch unnecessarily small. Additionally, in Figure 10 the exact and approximate values of $\|E_j p_j\|_{A^{-1}}$ agree very closely. Is this level of agreement generally expected or is there a specific reason that makes experiments considered particularly amenable.

**Acknowledgement** We would like to thank the two anonymous reviewers for their valuable comments which greatly improved our paper.

**Data availability statement** The python code for the numerical experiments shown in Section 5 can be found in a public GitHub repository. The link is as follows: https://github.com/rbcfc/parareal_data_assimilation.

## REFERENCES

[1] M. ASCH, M. BOCQUET, AND M. NODET, *Data assimilation: methods, algorithms, and applications*, SIAM, 2016.

[2] L. BAFFICO, S. BERNARD, Y. MADAY, G. TURINICI, AND G. ZÉRAH, *Parallel-in-time molecular-dynamics simulations*, Physical Review E, 66 (2002), p. 057701.

[3] G. BAL, *On the convergence and the stability of the parareal algorithm to solve partial differential equations*, in Domain decomposition methods in science and engineering, Springer, 2005, pp. 425–432.

[4] A. BOURAS AND V. FRAYSSÉ, *Inexact matrix-vector products in krylov methods for solving linear systems: a relaxation strategy*, SIAM J. Matrix Anal. Appl., 26 (2005), pp. 660–678.

[5] P. COURTIER, E. ANDERSSON, W. HECKLEY, D. VASILJEVIC, M. HAMRUD, A. HOLLINGSWORTH, F. RABIER, M. FISHER, AND J. PAILLEUX, *The ecmwf implementation of three-dimensional variational assimilation (3d-var). i: Formulation*, Q. J. R. Meteorol. Soc., 124 (1998), pp. 1783–1807.

[6] P. COURTIER, J.-N. THÉPAUT, AND A. HOLLINGSWORTH, *A strategy for operational implementation of 4d-var, using an incremental approach*, Q. J. R. Meteorol. Soc., 120 (1994),
pp. 1367–1387.

[7] B. CUSHMAN-ROISIN AND J.-M. BECKERS, *Introduction to geophysical fluid dynamics: physical and numerical aspects*, Academic press, 2011.

[8] R. DALEY, *Atmospheric data analysis*, no. 2, Cambridge university press, 1993.

[9] X. DU, M. SARKIS, C. E. SCHAERER, AND D. B. SZYLD, *Inexact and truncated parareal-in-time krylov subspace methods for parabolic optimal control problems*, Electron. Trans. Numer. Anal., 40 (2013), pp. 36–57.

[10] D. R. DURRAN, *Numerical methods for wave equations in geophysical fluid dynamics*, vol. 32, Springer Science & Business Media, 2013.

[11] C. FARHAT AND M. CHANDESRIS, *Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid-structure applications*, Internat. J. Numer. Methods Engrg., 58 (2003), pp. 1397–1434.

[12] M. FISHER AND S. GÜROL, *Parallelization in the time dimension of four-dimensional variational data assimilation*, Q. J. R. Meteorol. Soc., 143 (2017), pp. 1136–1147.

[13] M. FISHER, Y. TREMOLET, H. AUVINEN, D. TAN, AND P. POLI, *Weak-constraint and long-window 4D-Var*, ECMWF Reading, UK, 2012.

[14] M. J. GANDER, *50 years of time parallel time integration*, in Multiple shooting and time domain decomposition methods, Springer, 2015, pp. 69–113.

[15] M. J. GANDER AND E. HAIRER, *Nonlinear convergence analysis for the parareal algorithm*, in Domain decomposition methods in science and engineering XVII, Springer, 2008, pp. 45–56.

[16] M. J. GANDER, F. KWOK, AND J. SALOMON, *Paraopt: A parareal algorithm for optimality systems*, SIAM J. Sci. Comput., 42 (2020), pp. A2773–A2802.

[17] M. J. GANDER AND S. VANDEWALLE, *Analysis of the parareal time-parallel time-integration method*, SIAM J. Sci. Comput., 29 (2007), pp. 556–578.

[18] A. E. GILL, *Atmosphere-ocean dynamics*, vol. 30, Academic press, 1982.

[19] G. H. GOLUB AND Q. YE, *Inexact preconditioned conjugate gradient method with inner-outer iteration*, SIAM J. Sci. Comput., 21 (1999), pp. 1305–1320.

[20] S. GRATTON, E. SIMON, D. TITLEY-PELOQUIN, AND P. L. TOINT, *Minimizing convex quadratics with variable precision conjugate gradients*, Numer. Linear Algebra Appl., 28 (2021), p. e2337.

[21] J. R. HOLTON, *An introduction to dynamic meteorology*, American Journal of Physics, 41 (1973), pp. 752–754.

[22] K. IDE, P. COURTIER, M. GHIL, AND A. C. LORENC, *Unified notation for data assimilation: Operational, sequential and variational (gtspecial issueltdata assimilation in meteology and oceanography: Theory and practice)*, Journal of the Meteorological Society of Japan. Ser. II, 75 (1997), pp. 181–189.

[23] L. ISAKSEN, *Data assimilation on future computer architectures*, in Proc. ECMWF Seminar on Data Assimilation for Atmosphere and Ocean, 2012, pp. 301–322.

[24] E. KALNAY, *Atmospheric modeling, data assimilation and predictability*, Cambridge university press, 2003.

[25] A. LAWLESS, S. GRATTON, AND N. NICHOLS, *An investigation of incremental 4d-var using non-tangent linear models*, Q. J. R. Meteorol. Soc., 131 (2005), pp. 459–476.

[26] J.-L. LIONS, Y. MADAY, AND G. TURINICI, *Résolution d'edp par un schéma en temps pararéel*, Comptes Rendus de l'Académie des Sciences-Series I-Mathematics, 332 (2001), pp. 661–668.

[27] A. LORENC, S. BALLARD, R. BELL, N. INGLEBY, P. ANDREWS, D. BARKER, J. BRAY, A. CLAYTON, T. DALBY, D. LI, ET AL., *The met. office global three-dimensional variational data assimilation scheme*, Q. J. R. Meteorol. Soc., 126 (2000), pp. 2991–3012.

[28] Y. MADAY, M.-K. RIAHI, AND J. SALOMON, *Parareal in time intermediate targets methods for optimal control problems*, Control and optimization with PDE constraints, (2013), pp. 79–92.

[29] G. MADEC, R. BOURDALLÉ-BADIE, P.-A. BOUTTIER, C. BRICAUD, D. BRUCIAFERRI, D. CALVERT, J. CHANUT, E. CLEMENTI, A. COWARD, D. DELROSSO, ET AL., *Nemo ocean engine*, (2017).

[30] F. MAGOULÈS AND G. GBIKPI-BENISSAN, *Asynchronous parareal time discretization for partial differential equations*, SIAM J. Sci. Comput., 40 (2018), pp. C704–C725.

[31] F. MAGOULÈS, G. GBIKPI-BENISSAN, AND Q. ZOU, *Asynchronous iterations of parareal algorithm for option pricing models*, Mathematics, 6 (2018), p. 45.

[32] T. P. MATHEW, M. SARKIS, AND C. E. SCHAERER, *Analysis of block parareal preconditioners for parabolic optimal control problems*, SIAM J. Sci. Comput., 32 (2010), pp. 1180–1200.

[33] M. MINION, *A hybrid parareal spectral deferred corrections method*, Commun. Appl. Math.
Comput. Sci., 5 (2011), pp. 265–301.

[34] B. W. ONG AND J. B. SCHRODER, *Applications of time parallelization*, Comput. Vis. Sci., 23 (2020), pp. 1–15.

[35] F. RABIER, H. JÄRVINEN, E. KLINKER, J.-F. MAHFOUF, AND A. SIMMONS, *The ecmwf operational implementation of four-dimensional variational assimilation. i: Experimental results with simplified physics*, Q. J. R. Meteorol. Soc., 126 (2000), pp. 1143–1170.

[36] D. RUPRECHT, *Wave propagation characteristics of parareal*, Comput. Vis. Sci., 19 (2018), pp. 1–17.

[37] D. RUPRECHT AND R. KRAUSE, *Explicit parallel-in-time integration of a linear acoustic-advection system*, Comput. & Fluids, 59 (2012), pp. 72–83.

[38] V. SIMONCINI AND D. B. SZYLD, *Theory of inexact krylov subspace methods and applications to scientific computing*, SIAM J. Sci. Comput., 25 (2003), pp. 454–477.

[39] J.-N. THEPAUT AND P. COURTIER, *Four-dimensional variational data assimilation using the adjoint of a multilevel primitive-equation model*, Q. J. R. Meteorol. Soc., 117 (1991), pp. 1225–1254.

[40] J.-N. THÉPAUT, R. N. HOFFMAN, AND P. COURTIER, *Interactions of dynamics and observations in a four-dimensional variational assimilation*, Mon. Weather Rev., 121 (1993), pp. 3393–3414.

[41] G. K. VALLIS, *Atmospheric and oceanic fluid dynamics*, Cambridge University Press, 2017.

[42] J. VAN DEN ESHOF AND G. L. SLEIJPEN, *Inexact krylov subspace methods for linear systems*, SIAM J. Matrix Anal. Appl., 26 (2004), pp. 125–153.

[43] D. C. WILCOX ET AL., *Turbulence modeling for CFD*, vol. 2, DCW industries La Canada, CA, 1998.