# A new global ocean ensemble system at the Met Office: Assessing the impact of hybrid data assimilation and inflation settings

Daniel J. Lea, James While, Matthew J. Martin, Anthony T. Weaver, Andrea Storto, Martin Chrust

► To cite this version:

Daniel J. Lea, James While, Matthew J. Martin, Anthony T. Weaver, Andrea Storto, et al.. A new global ocean ensemble system at the Met Office: Assessing the impact of hybrid data assimilation and inflation settings. Quarterly Journal of the Royal Meteorological Society, 2022, 148 (745), pp.1996-2030. 10.1002/qj.4292. hal-04740255

HAL Id: hal-04740255

https://cnrs.hal.science/hal-04740255v1

Submitted on 23 Oct 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
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

# A new global ocean ensemble system at the Met Office:
## Assessing the impact of hybrid data assimilation and
inflation settings

Daniel J. Lea<fnref n="1" /> ID

James While<fnref n="1" /> ID

Matthew J. Martin<fnref n="1" /> ID

Anthony Weaver<fnref n="2" /> ID

Andrea Storto<fnref n="3" /> ID

Marcin Chrust<fnref n="4" />

**Correspondence**

Daniel J. Lea, Met Office, Exeter, EX1 3PB,
UK.

Email: daniel.lea@metoffice.gov.uk

**Funding information**

Copernicus Climate Change Service,
Grant/Award Number: C3S 321b; WCSSP
India

**Abstract**

We have developed a global ocean and sea-ice ensemble forecasting system based on the operational forecasting ocean assimilation model (FOAM) system run at the Met Office. The ocean model Nucleus for European Modelling of the Ocean (NEMO) and the community ice code (CICE) sea-ice model are run at 1/4° resolution and the system assimilates data using a three-dimensional variational assimilation (3DVar) version of NEMOVAR. This data assimilation (DA) system can perform hybrid ensemble/variational assimilation. A 36-member ensemble of hybrid ensemble variational assimilation systems with perturbed observations (values and locations) has been set up, with each member forced at the surface by a separate member of the Met Office Global-Regional Ensemble Prediction System (MOGREPS-G). The unperturbed member is forced by atmospheric fields from the Met Office operational numerical weather prediction (NWP) deterministic system. The system includes stochastic model perturbations and a relaxation to prior spread (RTPS) inflation scheme. A control run of the system using an ensemble of 3DVars is shown to be generally reliable for Sea-Level Anomaly (SLA), temperature, and salinity (the ensemble spread being a good representation of the uncertainty in the ensemble mean), although the ensemble is underspread in eddying regions. The ensemble mean gives a 4% reduction in error in SLA compared with the deterministic 3DVar system currently used operationally. The system was tested with different weights for the ensemble component of the hybrid background-error covariance matrix and different inflation factors. The best results, in terms of short-range forecast error and ensemble reliability statistics, were obtained with hybrid three-dimensional ensemble variational DA (3DEnVar). The RTPS inflation scheme is shown to be beneficial in producing an appropriate ensemble spread in response to hybrid DA. 3DEnVar with an ensemble hybrid weight of 0.8 leads to a reduction of 20% (5%) in the ensemble mean error for SLA (profile temperature and salinity) compared with an ensemble of standard 3DVars.

**KEYWORDS**

data assimilation, ensembles, global, ocean, variational

<footnote n="1">¹Met Office, Exeter, UK</footnote>

<footnote n="2">²CERFACS/CECI CNRS UMR 5318,
Toulouse, France</footnote>

<footnote n="3">³Italian National Research Council,
Institute of Marine Sciences, Rome, Italy</footnote>

<footnote n="4">⁴ECMWF, Reading, UK</footnote>
# 1 | INTRODUCTION

Most operational atmospheric forecasting systems have used ensemble methods for many years. In contrast, most global short-range ocean forecasting systems have generally not used ensembles. With a push towards coupled ocean-atmosphere numerical weather prediction (NWP) at many forecasting centres, there is increasing interest in developing global ocean ensemble forecasting capability in coupled systems. The ensemble approach is appealing, partly due to the fact that it provides a real-time estimate of the uncertainty in the model forecasts. This information can be used to improve the results of the data assimilation (DA), and thereby the forecasts, by providing a more accurate representation of the background-error covariances than the static or parametrised estimates used in non-ensemble DA. The spread of the ensembles can also be used with forecasts to quantify their uncertainty, which is useful to many users of forecast products. In coupled systems, the uncertainty in one fluid (e.g., the near-surface ocean) can also propagate into other components of the coupled system (e.g., the atmospheric boundary layer), thereby giving improved uncertainty estimates for the full coupled system. In the ocean, unlike in the atmosphere, the mesoscale (the synoptic-scale equivalent), while resolved by models, is not well constrained by observations (see e.g., Thoppil et al., 2021). Taking an ensemble approach can allow these unconstrained scales to be filtered out in initial conditions and forecasts. Another major benefit of ensemble methods is that they are highly parallelisable and therefore make effective use of computers with increasing numbers of processors.

A number of centres are developing, or are currently using, ensemble-based ocean and coupled ocean-atmosphere forecasting systems. The US Navy have developed a coupled ensemble DA (EDA) system using 3DVar in the ocean (Barton et al., 2020). The European Centre for Medium-Range Weather Forecasts (ECMWF) has a weakly coupled atmosphere-ocean-wave-ice-land ensemble system (Browne et al., 2019) using 3DVar in the ocean. Penny et al. (2015) developed a hybrid data assimilation approach for coupled seasonal forecasting at the National Center for Environmental Prediction (NCEP). Development of a global ensemble Kalman filter (EnKF) has also been undertaken in the Australian Bluelink system (Schiller et al., 2020), demonstrating significant improvements to the accuracy of the ocean forecasts. A number of regional ocean forecasting systems also use ensembles: for example, the TOPAZ system has been running an EnKF for many years (Sakov et al., 2012).

The Met Office will soon be implementing a major upgrade to its global NWP system. This will include a coupled model with evolving ocean and sea-ice components

as well as the existing atmosphere and land components. Both the high-resolution deterministic NWP system and the lower resolution Met Office Global-Regional Ensemble Prediction System (MOGREPS-G: Bowler et al., 2008) will be run coupled to the ocean. The initialisation of the deterministic coupled forecasts will use the weakly coupled data assimilation (DA) approach described by Lea et al. (2015). The ocean and sea-ice components of that system are based on those used in the operational global ocean forecasting system called the Forecasting Ocean Assimilation Model (FOAM: Blockley et al., 2014).

The aim of the work described here is to develop a global ocean and sea-ice ensemble forecasting system that provides a good representation of the uncertainty in short-range forecasts. We also aim to make use of ensemble forecast information to improve the forecast-error covariance information used by the ocean DA. The plan is to implement the developments in an upgrade to the coupled MOGREPS-G system. They should provide improved uncertainty information about sea-surface temperature (SST) and sea-ice to the atmosphere. This should also mean that the ocean/sea-ice ensemble in the coupled MOGREPS-G can provide realistic uncertainty information about ocean forecasts (including below the surface) to users. The improved coupled ensemble could also provide a better set of initial ocean/sea-ice conditions for ensemble seasonal forecasts compared with the existing lagged ensemble approach used in Met Office's seasonal prediction system (GloSea, described in MacLachlan et al., 2015).

The ocean data assimilation system used at the Met Office is based on NEMOVAR, which is a collaborative development with Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique (CERFACS), ECMWF, and Institut National de Recherche en Information et en Automatique (INRIA). The current operational implementation of NEMOVAR uses a first guess at appropriate time (3DVar-FGAT) method and is described by Waters et al. (2015). The background-error covariances in this scheme determine how much weight is given to new observations in the analysis, and how their information is spread in the horizontal and vertical directions and between variables. These background-error covariances are currently based on a combination of climatological estimates coming from previous long reanalyses and parametrisations based on the current background model state (see Waters et al. (2015) and Mirouze et al. (2016) for a detailed description of how NEMOVAR is configured in FOAM). However, a well-designed and validated ensemble could provide much better information about the “errors-of-the-day”, improving the background-error covariance specifications. This would allow us to make better use of ocean observations. This information could also
be useful for targeting observations, such as those from gliders, to regions of high uncertainty in order to improve the forecast accuracy.

The core ocean and sea-ice model and data assimilation system is based on the latest operational FOAM system and is described in Section 2, together with an overview of how it has been converted into an ensemble system. This includes details on the atmospheric forcing, the scheme for observation perturbations, and the stochastic model perturbation schemes. In Section 3 we describe how we use the ensemble information in a hybrid-3DEnVar scheme (whereby climatological background-error covariances are combined linearly with ensemble-based background-error covariances), as well as the ensemble inflation method. In Section 4 we describe the experimental setup and the methods used to assess the performance of the ensemble system. In this section we demonstrate that the baseline ensemble system using standard 3DVar with no inflation performs well. In the results, Section 5, we explore the impact of various weights for the ensemble part of the hybrid-3DEnVar compared with standard 3DVar and pure 3DEnVar (in which only ensemble-based background-error covariances are used). Finally, a summary and discussion can be found in Section 6.

## 2 | ENSEMBLE SYSTEM DESCRIPTION

### 2.1 | Deterministic system description

The developments described here are based on the latest version of the global FOAM system (Blockley et al., 2014). It is run on a daily cycle and comprises the Nucleus for European Modelling of the Ocean (NEMO) ocean model at version 3.6, the Community Ice Code (CICE) sea-ice model at version 8.1 and the NEMOVAR data assimilation system at version 6.2.6. A recent upgrade to the operational system (Aguiar et al., in preparation) has included the implementation of a 1/12° resolution model, but here we use a version of the model with an “ORCA” grid at approximately 1/4° resolution and 75 vertical levels with 1-m vertical resolution near the surface with increasing level thickness with depth. The ocean model is coupled to the CICE sea-ice model as described by Ridley et al. (2018), and is run in the GO6 configuration described by Storkey et al. (2018). The model is forced at the surface by atmospheric fields from the Met Office operational NWP deterministic system with approximately 10-km resolution at 3-hr frequency for the heat and freshwater inputs and 1-hr frequency for the winds. The freshwater runoff from land is specified as climatological monthly varying estimates.

The data assimilation component of FOAM is NEMOVAR and is implemented at the Met Office as a multivariate, incremental 3DVar-FGAT scheme. Waters et al. (2015) give a description of the scheme as implemented in an earlier FOAM version. The main updates to the scheme since then (apart from the implementation of the ensemble capability) include the implementation of a two length-scale background-error covariance (Mirouze et al., 2016) and the use of a two-dimensional by one-dimensional implicit diffusion operator for modelling the three-dimensional background-error correlations (Weaver et al., 2016; 2021).

The data assimilation in FOAM uses a one-day time window. During the nonlinear trajectory, the model counterparts of the observations are calculated at the nearest model time step to the observation time. The analysis increments generated by NEMOVAR are subsequently used in a run of NEMO over the same day. The increment is converted into a forcing and added throughout the one day period to the model tendencies using the technique of incremental analysis updates (IAU: Bloom et al., 1996). The pressure-correction scheme of Bell et al. (2004) is used to reduce spurious circulations that can be generated when assimilating data near the Equator.

The observation types assimilated include in situ and satellite SST data, satellite altimeter sea-level anomaly (SLA) data, satellite sea-ice concentration data, and in situ profiles of temperature and salinity from various platforms. A variational observation bias correction scheme is used, as described for SST by While and Martin (2019) and for sea-surface height (SSH) by Lea et al. (2008). The SSH bias-correction scheme was developed to account for errors in the mean dynamic topography. A scheme was recently introduced to deal with discrepancies in the impact of the atmospheric forcing on the model’s barotropic SSH and the way this is dealt with in the observation preprocessing (see Aguiar et al., in preparation, for more details). The observations used here were quality-controlled by the FOAM operational suite. This is done using the procedure described by Ingleby and Huddleston (2007), which includes a background check. The implementation of this in FOAM is as described in Storkey et al. (2010).

### 2.2 | Atmospheric forcing ensemble

The ensemble system runs 36 perturbed members in addition to a deterministic run. Each of the ensemble members differs from the deterministic run in the atmospheric forcing, the observations that are assimilated, and ocean model differences. The way in which these three sources of uncertainty are represented in our
ensemble system is described in this and the following
subsections.

A major source of uncertainty in ocean forecasts comes
from errors in the atmospheric forcing. These uncertain-
ties could be represented by perturbing the determin-
istic atmospheric forcing using methods such as those
described by Mirouze and Storto (2019). However, since
we have access to an atmospheric ensemble and the aim
of this work is to implement the ocean ensemble within a
coupled ensemble system, we use the atmospheric ensem-
ble directly here. Each ocean/sea-ice ensemble mem-
ber is forced by a particular member of the operational
MOGREPS-G system. MOGREPS-G is at N640 resolu-
tion (about 20 km) during the period of the experiments
described in this article, compared with the determinis-
tic NWP resolution at N1280 (about 10 km). MOGREPS-G
runs long forecasts, suitable for a 24-hr cycling ocean, for
18 members on each 6-hr cycle. As this work is moti-
vated to prepare for work to include an ocean ensemble in
the Met Office coupled NWP suite (which has 44 ensem-
ble members), we want to run more than 18 members
here. To provide forcing for 36 ocean ensemble members,
we use the MOGREPS-G forecasts initialised at 0000 and
1800 UTC.

## 2.3 | Observation perturbations

The overall aim of our work is to develop an ensem-
ble system such that the spread in short-range forecasts
generated from the ensemble provides a good approxi-
mation to the uncertainty in the ensemble mean fore-
cast. The spread in initial conditions is an important part
of such a system, especially since we would like to use
the short-range (1-day) forecast ensemble for specifying
the background-error covariance matrix in a hybrid data
assimilation scheme. The surface forcing perturbations
alone, while an important component of the ensemble
generation, are unlikely to generate ensemble spread that
is representative of the uncertainty everywhere, particu-
larly in the deep ocean.

There are various methods described in the literature
for generating an ensemble of initial conditions. Given that
the aim is for the system developed here to become part of
a coupled ensemble prediction system, we have followed
the overall design of the atmospheric ensemble genera-
tion system. In the Met Office NWP system, Bowler *et al.*
(2017) describe an ensemble of 4DEnVar data assimila-
tion systems with observation perturbations generated by
randomly sampling observation uncertainty using the data
assimilation system’s observation-error covariance matrix
**R**. Here we also implement an ensemble of data assimila-
tion systems with perturbed observations.

In the data assimilation context, observation errors are
made up of measurement errors and representation errors
(which account for the fact that models do not repre-
sent all the physical processes and scales that a particular
observation measures; Janjic *et al.*, 2018). In the global
FOAM system used here, the observation-error covariance
matrix used in the data assimilation, **R**, is assumed to
be diagonal (no spatial, temporal, or interplatform error
correlations are represented). The observation-error vari-
ances are specified as spatially (horizontally and vertically)
and seasonally varying estimates, which were estimated
from the innovations in previous reanalysis integrations
using the method described in Waters *et al.* (2015). That
method estimates the total observation errors, including
both measurement and representation errors.

Rather than following the approach in the atmospheric
system and directly using the statistics of the **R** matrix,
we generate observation perturbations aimed at mimick-
ing the measurement and representation errors separately.
The measurement error is included by randomly perturb-
ing the observation value with a standard deviation given
by the known measurement errors, as given in Table 1.
For in situ data, this is a very small component of the
total error, except in the deep ocean. The representation
error is included by randomly perturbing the observa-
tion location, following a similar approach to that of (Zuo
et al., 2017; 2018). In that system, observation locations
are perturbed by selecting different observations in the
horizontal and vertical during the observation-thinning
preprocessing step. Here we perturb the location of the
observations explicitly with a standard deviation of 0.1°
in latitude and longitude. This relates to the scales that
the model can represent, and so for a higher resolu-
tion model generally the perturbation distances should
be reduced. During the application of the observation
operator, the model values are interpolated to these per-
turbed locations, thereby generating perturbations in the
innovations (observation-minus-model background val-
ues), which represent the subgrid-scale variability.

For profile data, simply adding random perturbations,
without any vertically correlated component, results in
unrealistic profiles, which are likely to be statically unsta-
ble. The data assimilation will act to filter out much of
the noise imposed by such schemes, and, when applied
to the model, much of the remaining noise in the vertical
is likely to be removed by vertical mixing processes. Per-
turbing observation locations as we do here has the advan-
tage over simple random perturbations of the observed
profile values that the resulting observed profiles are real-
istic representations of the true ocean. For surface ocean
observations, the measurement errors are a bigger frac-
tion of the total observation error, but we apply the same
scheme of perturbing observation locations and values for
<figure><img src="image_3.png" /></figure>

<table><caption>TABLE 1 Measurement-error standard deviations used to perturb observation values</caption><thead><tr><th>Variable (unit)</th><th>Platform</th><th>Measurement-error standard deviation</th></tr></thead><tbody><tr><td>SLA (m)</td><td>Jason-3, Sentinel-3A, Cryosat-2, Altika</td><td>0.04</td></tr><tr><td rowspan="4">SST (°C)</td><td>AVHRR/NOAA*</td><td>0.39/0.42/1.08</td></tr><tr><td>AVHRR/MetOp*</td><td>0.31/0.39/0.55</td></tr><tr><td>VIIRS/Suomi-NPP*</td><td>0.13/0.32/1.04</td></tr><tr><td>Ships, drifting and moored buoys</td><td>0.3</td></tr><tr><td>Sea-ice concentration</td><td>SSMI/S</td><td>0.05</td></tr><tr><td rowspan="5">Temperature profiles (°C)</td><td>Argo</td><td>0.002</td></tr><tr><td>Marine mammals</td><td>0.002</td></tr><tr><td>CTDs</td><td>0.002</td></tr><tr><td>XBTs</td><td>0.15</td></tr><tr><td>Moored buoys</td><td>0.05</td></tr><tr><td rowspan="3">Salinity profiles (PSS)</td><td>Argo</td><td>0.01</td></tr><tr><td>CTDs</td><td>0.005</td></tr><tr><td>Moored buoys</td><td>0.02</td></tr></tbody></table>

Note: *For satellite SST data, per-measurement standard deviations (single-sensor error statistics (SSES)) are used, as provided with the data from the Group for High Resolution SST (GHRSSST), the minimum/mean/maximum of which are indicated in the table for a 10-day period at the start of March 2018.

all data types. For sea-ice concentration, we cap the perturbed value to ensure it is not greater than one or less than zero. Note that the SST and SSH observation bias correction schemes are run for each ensemble member. As the model state varies in each ensemble member and the observations are perturbed, there is also some spread of observation bias estimates in the ensemble.

The recent article of van Leeuwen (2020) develops a theoretical framework on observation perturbations for the first time and demonstrates that the correct approach is in fact to perturb the modelled observations (output from the observation operator), not the observations themselves (as we have done in the work here). The perturbation of the observation locations can be seen as a perturbation of the model equivalent value, so is unchanged by adopting the van Leeuwen (2020) approach. For the random value perturbations part, we should instead be applying these perturbations to the model equivalent value. In most situations this would not change the results, except where the error distribution is asymmetric, such as where a cap is imposed. This change would thus have most effect on sea-ice concentration and for salinity when the salinity is very low. We intend to implement this in future versions of the ocean ensemble system.

To test the impact of the different observation perturbation strategies, three ensemble runs were carried out: one with only the atmospheric forcing ensemble (called no_obperts); one with the atmospheric forcing ensemble and observation value perturbations based directly on the

R matrix (called Rval_perts); and one with the atmospheric forcing ensemble, observation values perturbed based on Table 1, and locations perturbed with a standard deviation of 0.1° in latitude and longitude (called valloc_perts). The impact of the perturbed observations over the first month is shown in Figures 1–3 for SST, SLA, and profiles, respectively. Note that, as the spin-up time for the ensemble system is longer than one month, these results are indicative only. In these runs, the initial state of all ensemble members was the same, so these figures show the impact of the observation perturbations during the spin-up of the ensemble.

Figure 1 shows that the observation value and location perturbations result in a significant increase in SST ensemble spread in areas of greater variability, such as the western boundary currents and the Antarctic circumpolar current (ACC), where we expect the forecasts to be more uncertain. At the same time, the observation value and location perturbations slightly reduce the root-mean-square (RMS) error in the ensemble mean, particularly in those same highly variable regions. When only value perturbations are applied based on the R matrix, there is a general increase in the SST ensemble spread over most areas (with less spread in the Tropics), but the RMS error in the ensemble mean is not affected significantly. The globally averaged standard deviation specified in R for SST is 0.63 °C and there is some spatial variability in its definition, with higher values in the more variable regions (not shown). This is not reflected in the change in
<figure><img src="image_4.png" /><figcaption>FIGURE 1 Plots (a)–(e) show the RMS difference between the ensemble mean and the observations, plots (f)–(j) show the ensemble spread at observed locations. The left plots, (a) and (f), are for the *no_obperts* run, plots (b) and (g) are for the *Rval_perts* run, and plots (d) and (i) are for the *valloc_perts* run (see the main text for a description of these experiments). The right plots show the difference between the observation perturbation runs and the one with no observation perturbations. These are calculated over March 2018 in 2° × 2° latitude/longitude bins for SST (°C) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

ensemble spread, perhaps due to vertical smoothing of perturbations by the assimilation and model, or perhaps due to inappropriate specifications in R for the other observation types assimilated (which will also affect the spread in SST). However, the measurement and location perturbation method shows good performance, so we focus on that method in the rest of the article.

The SLA spread is increased by the measurement and location perturbation method in similar regions to the SST spread, as shown in Figure 2. There is also a small

reduction in the RMS error in the ensemble mean SLA in those regions. The temperature and salinity spread is also significantly increased by the measurement and location perturbation method, as shown in Figure 3. The temperature spread is approximately doubled over most of the top 1000 m, with smaller but still significant increases below that. The increases in temperature and salinity spread result in only minor changes to the RMS error in the ensemble mean (not shown). It should be noted that these results are only for the first month of the ensemble spin-up,
<figure><img src="image_5.png" /></figure>

FIGURE 2 Plots (a)–(c) show the RMS difference between the ensemble mean and the observations, plots (d)–(f) show the ensemble spread at observed locations. The left plots, (a) and (d), are for the *no_obperts* run, and plots (b) and (e) are for the *valloc_perts* run (see the main text for a description of these experiments). The right plots show the difference between the observation perturbation runs and the one with no observation perturbations. These are calculated over March 2018 in 2° × 2° latitude/longitude bins for SLA (m) [Colour figure can be viewed at wileyonlinelibrary.com]

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

FIGURE 3 Globally averaged ensemble spread at observed locations for potential temperature (°C) and salinity (psu) calculated over the first month of ensemble spin-up (March 2018). The red line is for the *no_obperts* run and the blue line is for the *valloc_perts* run [Colour figure can be viewed at wileyonlinelibrary.com]
and we might expect larger impacts once the ensemble is
fully spun up, as shown later in the article.

**2.4 | Stochastic model perturbations**

Models contain many approximations, due to the resolu-
tion at which the physical equations are discretised, the
parametrisation of unresolved processes, and the numeri-
cal schemes. These various approximations result in model
error, which contains both random and systematic com-
ponents. The aim of stochastic model perturbations is to
represent uncertainties in the model due to these approxi-
mations, thereby giving an ensemble that reflects the true
uncertainty in the forecast. The inclusion of such schemes
can also reduce the error in the ensemble mean (Leut-
becher et al., 2017).

Three common schemes used for perturbing models
in ensemble systems have been implemented in NEMO
by Storto and Andriopoulos (2021) and are used here.
The stochastically perturbed parametrisation tendencies
(SPPT) scheme (Buizza et al., 1999; Sanchez et al., 2016)
adds a perturbation to the total parametrised physics
tendencies that is the product of the unperturbed ten-
dencies by a random field. The stochastically perturbed
parameters (SPP) scheme (Ollinaho et al., 2017) intro-
duces stochastic perturbations to model parameters within
parametrised processes. The stochastic kinetic energy
backscatter (SKEB) scheme (Tennant et al., 2011) adds per-
turbations to the barotropic streamfunction proportionally
to some sinks of energy on small scales in the model, which
in a higher resolution model would otherwise have been
backscattered to the scales we can resolve. Below we give
a brief overview of each of these three schemes, and the
reader is referred to Storto and Andriopoulos (2021) for
more information.

In the SPPT scheme, model tendencies due to various
model processes are accumulated at each time step. At the
end of the tracer and dynamics updates for each time step,
the tendencies are updated, adding a term that is the prod-
uct of the unperturbed tendencies associated with subgrid
processes, by a Gaussian field with zero mean and a spec-
ified standard deviation, which is correlated in space and
time.

The random field (which is the same for all variables)
is specified with a standard deviation of 0.1, and 150 passes
of a Shapiro filter are used to generate spatial correla-
tions in the field (giving smoothing with a length-scale of
~600 km). The random field is set to be zero at the sur-
face, half the specified standard deviation at the second
level, then the full standard deviation at the third level and
below. A similar tapering of the standard deviation occurs
at the bottom.

The scheme allows for many different processes to con-
tribute to the SPPT perturbations. In our implementation
we include the lateral and vertical tracer diffusion, lat-
eral and vertical viscosity, and tracer damping. A test of
the impact of the SPPT scheme on the ensemble spread
was carried out (not shown), in which the maximum dif-
ferences were of the order of 0.2°C, 0.02 practical salinity
scale (PSS), 20 m, 4 cm·s⁻¹ and 0.5 cm for SST, sea surface
salinity (SSS), mixed-layer depth, surface current speed,
and SSH, respectively.

In the SPP scheme, the parameter values used in
specific parametrised process are perturbed. Instead of
a single number parameter, used currently in NEMO, a
two-dimensional stochastic version of the parameter is
used, which embeds spatial and temporal correlations. The
two-dimensional random field follows a log-normal distri-
bution that can be specified so that the mean of the result-
ing parameters is the same as the original (unperturbed)
parameter value (Ollinaho et al., 2017). The random field
is generated using a Shapiro filter as for SPPT, but with 300
iterations (giving a smoothing length-scale of ~ 800 km),
and a time decorrelation scale of 10 days.

A list of the model parameters that are perturbed is
given in Table 2, together with their unperturbed value and
the relative standard deviation used to generate the ran-
dom field of perturbations to the parameters. The param-
eters we perturb include the vertical mixing coefficients
for tracers and dynamics, diffusivity, viscosity, mixing near
river mouths, solar penetration, surface salinity relaxation,
and various parameters of the turbulent kinetic energy
(TKE) vertical mixing scheme.

In some short tests, the largest impacts of the SPP
scheme (not shown) are of the order of 0.1°C, 0.005 PSS,
5 m, and 1 cm·s⁻¹ for temperature, salinity, mixed-layer
depth, and current speed respectively. There was very little
difference in the sea-surface height fields.

Finally, the SKEB scheme aims to backscatter energy
that is dissipated at small scales to the scales resolved by
the model. The model processes that dissipate the energy
considered in the current implementation are numerical
advection schemes and the parametrisation of convection.
Numerical dissipation occurs due to explicit diffusion and
also due to approximations in the advection scheme. It is
expected to be larger when the shear and tension are large,
and it is therefore formulated as a function of these.

The convective dissipation component of SKEB is
meant to represent the energy that should be injected
by convection, which is missing from the convective
parametrisation. The numerical and convective dissipa-
tion terms are combined to calculate a perturbation to
the streamfunction that is the product of such dissipated
energy, scaled by the backscatter ratio, and a stochastic
forcing, generated using the same spatially and temporally
<table><caption>TABLE 2 A list of the active components of SPP used here. All the parameters listed are perturbed with a log-normal distribution where the mean of the resulting parameter set is the same as the unperturbed parameter</caption><thead><tr><td>Description</td><td>Unperturbed value</td><td>Relative standard deviation</td></tr></thead><tbody><tr><td>Vertical eddy diffusivity</td><td>1.2 × 10<sup>-5</sup> m<sup>2</sup>·s<sup>-1</sup></td><td>0.2</td></tr><tr><td>Vertical eddy viscosity</td><td>1.2 × 10<sup>-4</sup> m<sup>2</sup>·s<sup>-1</sup></td><td>0.2</td></tr><tr><td>Horizontal diffusivity (T, U, V, W)</td><td>150 m<sup>2</sup>·s<sup>-1</sup></td><td>0.2</td></tr><tr><td>Horizontal bi-Laplacian viscosity</td><td>-1.5 × 10<sup>11</sup> m<sup>4</sup>·s<sup>-1</sup></td><td>0.2</td></tr><tr><td>Mixing close to river mouth</td><td>2 × 10<sup>-3</sup> m<sup>2</sup>·s<sup>-1</sup></td><td>0.1</td></tr><tr><td>Solar penetration depth of extinction</td><td>0.35 m</td><td>0.1</td></tr><tr><td>SSS relaxation coeff.</td><td>-33.33 mm·day<sup>-1</sup></td><td>0.2</td></tr><tr><td>TKE Langmuir cell coefficient</td><td>0.15</td><td>0.1</td></tr><tr><td>TKE Kolmogoroff dissipation coeff.</td><td>0.7</td><td>0.1</td></tr><tr><td>TKE fraction of surface TKE below thermocline</td><td>5 × 10<sup>-2</sup></td><td>0.1</td></tr></tbody></table>

correlated random-field generation method as described above for the SPPT and SPP schemes. Comparisons of NEMO surface current outputs at 1/12° and 1/4° show that the 1/4° resolution configuration has much less energy below scales of about 250 km (Kjellsson and Zanna, 2017). We use 40 iterations of the Shapiro filter (giving a smoothing length-scale of ~300 km) to ensure these scales are attenuated by 50%.

Perturbations to the model's velocity fields are then calculated from the stochastic streamfunction and added to the model's velocity fields at each time step after the other dynamics updates have been applied. We set the backscatter ratio to 50% in our implementation: namely, half of the energy from numerical and convective dissipation is re-injected at the scales resolved by the model. In tests during March 2018 (not shown), the impact of SKEB was greatest on the currents in the North Atlantic, particularly in the Gulf Stream and the Labrador Sea. These are areas where the velocity gradients are particularly large, and where there is deep convection occurring at this time of year. There are also some smaller magnitude perturbations near the Equator. The maximum differences in the velocity fields were less than 0.5 cm·s⁻¹.

### 2.4.1 | Impact of combining the stochastic model perturbation schemes

To demonstrate the impact of all three of the stochastic model perturbation schemes together, we ran two ensemble experiments starting from 1st March 2018. The first included the atmospheric forcing ensemble and observation perturbations described in the previous two subsections. The second also applied the stochastic model perturbations. Figure 4 shows the impact of the

stochastic model perturbations on the ensemble spread in temperature and salinity, averaged over the last 10 days of the second month of the runs. The ensemble spread of SST is increased in regions where it was previously underspread, such as the Gulf Stream and regions in the Tropics. The SSS ensemble spread is increased in the Amazon outflow region and in the tropical Atlantic more generally. There is also a large impact of the stochastic model perturbations on the SSS near sea-ice, particularly in the Arctic. The maximum globally averaged increase in temperature spread occurs at about 100-m depth, with less impact from stochastic model perturbations below that depth. The maximum globally averaged increase in salinity spread occurs at the surface, with the schemes having less impact below the surface. The 108-m temperature spread is particularly large in the Amazon outflow region in the central tropical Pacific, and in the Brazil-Falkland Confluence Zone.

The SSH spread (not shown) is largely similar in the two runs, with only small increases in the Gulf Stream, Labrador Sea, and ACC. Mixed-layer depth spread (not shown) is increased over large areas of the North Atlantic, particularly in the storm-track region and in the Labrador Sea. The two velocity components (not shown) have increased spread in the western boundary currents, the ACC, and near the Equator, particularly in the zonal velocity in the equatorial Pacific.

## 3 | HYBRID DATA ASSIMILATION

A brief description of the implementation of NEMOVAR in the Met Office FOAM system is provided in Section 2.1. Weaver et al. (2018) describe in detail the developments made to NEMOVAR to allow the use of ensemble information in a hybrid ensemble/variational framework.
<figure><img src="image_7.png" /><figcaption>FIGURE 4 Difference in ensemble standard deviation between an ensemble with stochastic physics and an ensemble run without stochastic physics, averaged over 21–30 April, 2018. The left panel shows temperature (°C) and the right panel shows salinity (PSS); the top row is at the surface, the middle row is at 108-m depth, and the bottom row shows the globally averaged spread at each depth down to 2000 m. Positive values indicate a larger spread with stochastic physics [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<figure><img src="image_8.png" /></figure>

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

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

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

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

FIGURE 5 A time series of global departure and spread, at observation locations, in the spin-up control run for (a) SLA in m, (b) all SST in °C, (c) profile temperature in °C, and (d) profile salinity in PSS. The blue lines show the departure and the orange lines the ensemble spread [Colour figure can be viewed at wileyonlinelibrary.com]
<figure><img src="image_13.png" /><figcaption>FIGURE 7 Reliability statistics of the control ensda00_info00 for in situ SST (°C 2) calculated over Oct–Dec 2018. Note that a uniform observation uncertainty of 0.3 °C is assumed [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

where δ**w** = (K<sub>w</sub><sup>b</sup>)<sup>-1</sup>δ**x**<sub>w</sub>; δ**x**<sub>w</sub> = **x**<sub>w</sub> - **x**<sub>w</sub><sup>f</sup>; **x**<sub>w</sub><sup>f</sup> is the model forecast of temperature, salinity, and SSH; **d** = **y** - H(**x**<sub>w</sub><sup>f</sup>) is the innovation vector; **y** is the observation vector; **H** is the observation operator; and **H** is the observation operator linearised about **x**<sub>w</sub><sup>f</sup>. The matrices **B** and **R** are estimates of the background- and observation-error covariance matrices, respectively. The matrix K<sub>w</sub><sup>b</sup> is the invertible, lower-triangular block matrix element of the complete linearised balance operator:

$$K_x^b = \begin{pmatrix} K_w^b \\ K_u^b \end{pmatrix},$$

where K<sub>u</sub><sup>b</sup> is that component of the linearised balance operator that acts on **w** to produce a balanced velocity increment. If δ**w**<sup>a</sup> is the minimising solution of Equation (1), then δ**x**<sup>a</sup> = K<sub>x</sub><sup>b</sup>δ**w**<sup>a</sup> is the complete analysis increment (which includes an update to horizontal velocity) that is used to correct the forecast model state. In what follows, we will only consider the component K<sub>w</sub><sup>b</sup>, so we will drop the index **w** to ease the notation.

The quadratic cost function is minimised iteratively using a **B**-preconditioned conjugate gradient (BCG algorithm: Gürol *et al.*, 2014) where the first guess for the minimisation is a zero increment. A practical advantage of this particular CG algorithm is that it can be

applied without having to specify **B**<sup>-1</sup> (only **B** itself is required).

## 3.1 | Modelled error covariances

Historically, the background-error covariance matrix in NEMOVAR has been modelled as

$$\mathbf{B}^m = \mathbf{\Sigma}^m \mathbf{C}^m \mathbf{\Sigma}^m, \quad (2)$$

where the superscript m denotes “modelled” or parametrised covariances. The matrix Σ<sup>m</sup> is a diagonal matrix of the estimated background-error standard deviations of the control variables. The matrix C<sup>m</sup> is a correlation matrix that is defined to be block-diagonal with respect to the different control variables. In FOAM, the correlations are modelled using an implicit diffusion operator, formed as a symmetric product of a 2D horizontal and 1D vertical diffusion operator (Weaver *et al.*, 2016; 2021). For the 2D horizontal correlations, the diffusion operator is defined as a linear combination of two diffusion operators, where each operator represents the error correlations for a different band of spectral scales (Mirouze *et al.*, 2016). The length-scales used to model C<sup>m</sup> and the background-error standard deviations specified in Σ<sup>m</sup> are based on a combination of parametrisations
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 8 Reliability statistics of the control ensda00_inf00 for potential temperature (denoted POTM) profiles (°C 2) calculated over Oct–Dec 2018. Observation uncertainty as used in the assimilation [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

and statistical estimates calculated from outputs of previous long runs, as described in Section 2.1. The vertical length-scale and variance parametrisations in $B^m$ depend on the background state, which means that the error covariance specifications contain some flow dependence. This can be expected to be beneficial to the analysis if the background state is a reasonable estimate of the true state, but could exacerbate errors if it is a poor estimate of the true state.

## 3.2 | Using ensembles in the background-error covariances

Another possibility for estimating $B$ is to use information from an ensemble. Providing the ensemble accurately represents a sample of the true background error, a covariance based on ensemble information can account for the evolution of background errors from one assimilation cycle to the next. For example, if a new observing network is introduced, the ensemble will automatically adapt and produce new background errors. Certain locations (e.g., near fronts or at the ice edge) have larger background errors, which should be captured by the ensemble. While we attempt to account for some of these variations in $B^m$, an ensemble should provide additional information on the background errors.

Two options are available in NEMOVAR to use ensemble information. The first option uses the ensemble to adjust the parameters of the diffusion-modelled correlations and the background-error variances described in Equation (2). This can be done in a hybrid scheme, whereby these parameters are a combination of climatological and ensemble-based estimates. This approach is being explored at ECMWF. The second option uses the ensemble directly in a sample covariance estimate of $B$. This option allows us to capture the rich spatial patterns of the ensemble covariances, but admits a significant amount of sampling error due to the limited size of the ensemble. Covariance localisation is required to mitigate the effects of sampling error. This second option is the one we explore in this study. Details of both approaches are found in Weaver et al. (2018).

Consider the centred perturbation

$$ \epsilon_k^f = \mathbf{x}_k^f - \bar{\mathbf{x}}, \qquad (3) $$

where $\mathbf{x}_k^f$ is the forecast of the model state vector for ensemble member $k$, and

$$ \bar{\mathbf{x}} = \frac{1}{N_e} \sum_{n=1}^{N_e} \mathbf{x}_n^f $$

is the ensemble mean, which is estimated using all $N_e$ perturbed ensemble members. The unperturbed member
<figure><img src="image_16.png" /><figcaption>FIGURE 9 Reliability statistics of the control ensda00_info0 for salinity (denoted PSAL) profiles (PSS²) calculated over Oct–Dec 2018. Observation uncertainty as used in the assimilation [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

(k = 0) is excluded from the ensemble mean, since it uses higher resolution atmospheric forcing and thus may have different underlying statistics compared with the perturbed members.

The forecast of the control vector for ensemble member k, **w**<sub>k</sub><sup>f</sup>, is related to **x**<sub>k</sub><sup>f</sup> via the nonlinear balance operator:

$$ \mathbf{x}_{k}^{\mathrm{f}}=\mathrm{K}^{\mathrm{b}}\left(\mathbf{w}_{k}^{\mathrm{f}}\right). \qquad (4) $$

For estimating the sample covariance matrix for ensemble member $\ell$, we exclude the perturbed member $\mathbf{x}_{\ell}^{\mathrm{f}}$ to avoid the effects of “inbreeding” on the ensemble spread (Bowler *et al.*, 2017). The perturbed member $\ell$ is retained for linearising the balance operator in Equation (4), as well as for computing the mean. Substituting Equations (4) in (3) and linearising $\mathrm{K}^{\mathrm{b}}$ about $\mathbf{x}_{\ell}^{\mathrm{f}}=\mathrm{K}^{\mathrm{b}}(\mathbf{w}_{\ell}^{\mathrm{f}})$ yields

$$ \epsilon_k^f \approx \mathbf{K}_\ell^b \mathbf{e}_k^f, \qquad (5) $$

where $\mathbf{K}_\ell^b$ is the linearised balance operator and $\mathbf{e}_k^f = \mathbf{w}_k^f - \bar{\mathbf{w}}$. The sample covariance matrix for ensemble member $\ell$ is then defined as

$$ \mathbf{B}_{\ell}=\frac{1}{N_{\mathrm{e}}^{\prime}-1} \sum_{k=1; k \neq \ell}^{N_{\mathrm{e}}} \mathbf{e}_{k}^{\mathrm{f}}\left(\mathbf{e}_{k}^{\mathrm{f}}\right)^{\mathrm{T}}, \qquad (6) $$

where $N_e' = N_e$ for ensemble member 0, $N_e' = N_e - 1$ for the other (perturbed) ensemble members, and, using the approximation in Equation (5),

$$ \mathbf{e}_{k}^{\mathrm{f}}=\left(\mathbf{K}_{\ell}^{\mathrm{b}}\right)^{-1} \epsilon_{k}^{\mathrm{f}}. \qquad (7) $$

In compact form, we can express the sample covariance matrix for ensemble member $\ell$ as

$$ \mathbf{B}_{\ell}=\mathbf{X}_{\ell}^{\prime} \mathbf{X}_{\ell}^{\prime \mathrm{T}}, \qquad (8) $$

where

$$ \mathbf{X}_{\ell}^{\prime}=\frac{1}{\sqrt{N_{\mathrm{e}}^{\prime}-1}}\left(\mathbf{e}_{1}^{\mathrm{f}} \cdots \mathbf{e}_{\ell-1}^{\mathrm{f}} \mathbf{e}_{\ell+1}^{\mathrm{f}} \cdots \mathbf{e}_{N_{\mathrm{e}}}^{\mathrm{f}}\right), \qquad (9) $$

the prime indicating that these are ensemble perturbations centred about the ensemble mean and the subscript $\ell$ indicating that ensemble member $\ell$ has been excluded.

We aim to use the sample covariance of the ensemble directly in the cost function, so that we can make use of the anisotropic spatial structures present in the ensemble covariances, since these are not easily modelled using $\mathbf{B}^m$. However, we only have a limited-size ensemble, which leads to a large amount of sampling error in the ensemble covariance statistics. The two standard methods for reducing the effects of sampling error in ensembles are to localise the covariances (Hamill *et al.*, 2001; Houtekamer
<figure><img src="image_17.png" /><figcaption>FIGURE 10 Reliability statistics of the hybrid DA experiment ensda08_inf08 for SLA (m²) calculated over Oct–Dec 2018. Observation uncertainty as used in the assimilation (except no coastal amplification) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 11 Time series of the RMS (solid lines) and mean (dashed lines) in the ensemble mean one-day forecast errors compared with observations. Typical 95th percentile uncertainties for the RMS values are estimated using the bootstrapping method of Efron (1979) with 1000 resamples and are as follows: (a) SLA 0.0002 m, (b) SST 0.002 °C, (c) profile temperature 0.004°C, (d) profile salinity 0.001 PSS [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>FIGURE 12 Time series of the ensemble spread in observation space [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

and Mitchell, 2001; Ménétrier and Auligné, 2015), since the sampling error forms a larger part of the estimate at large spatial distances, and to combine the ensemble covariance estimates with the existing modelled covariance in a hybrid B scheme (Clayton et al., 2013; Bannister, 2017).

## 3.3 | Ensemble localisation

The estimate provided by Equation (8) is severely low-rank in practice, since the number of ensemble members is much smaller than the size of the control vector. If used directly, it would in turn severely limit the degrees of freedom of the assimilation and thus lead to problems fitting the observations. For this reason, and because the sampling error is larger than the useful signal from the ensemble at larger separation distance, a localisation procedure is applied. In variational assimilation, this is usually done by performing an element-by-element multiplication (Schur product) of the sample covariance matrix with a spatially limiting localisation correlation matrix L.

In NEMOVAR, localisation is performed on the sample correlation matrix of the control variables, $\tilde{X}'_e \tilde{X}'_e^T$, where

$\tilde{X}'_e = \Sigma_e^{-1} X'_e$ is the matrix of ensemble perturbations normalised by the diagonal matrix of sample standard deviations $\Sigma_e$ (the square root of the diagonal elements of $B_e$ in Equation 8). The form of the localised sample covariance matrix is then

$$B_e^e = \Sigma_e^e (L \circ \tilde{X}'_e \tilde{X}'_e^T) \Sigma_e^e, \quad (10)$$

where $\circ$ denotes the Schur product and $\Sigma^e$ is a diagonal matrix of standard deviations. Since the diagonal elements of both terms of the Schur product are all equal to one, the diagonal elements of their Schur product are also all equal to one. This formulation provides useful flexibility, as it allows, for example, separate filtering or hybridisation of the sample variances in the specification of $\Sigma^e$. If $\Sigma^e = \Sigma_e$, as in this study, then Equation (10) is equivalent to localising the sample covariance matrix directly.

Localisation can disrupt the balance properties of the analysis (Kepert, 2009). This problem is reduced in our formulation, since we apply localisation to the sample correlation matrix of the control vector and include the linearised balance operator as part of the generalised observation operator in the cost function (Equation (1)). In
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>FIGURE 13 Ensemble departure (standard deviation of the difference between ensemble mean and observations before assimilation) and ensemble spread (at observation locations) for SLA (in m), SST (in °C), profile temperature (in °C), and profile salinity (PSS) as a function of hybrid weight $\beta_e^2$ for experiments ensda00_inf00, ensda02_inf08, ensda05_inf08, ensda08_inf08, and ensda10_inf08. The unfilled circles indicate the values for the non-inflation experiment ensda08_inf00. Also plotted is the standard deviation of the innovations for the unperturbed member (for comparison with the ensemble mean departure). The statistics are for the global ocean from October–December 2018 [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

this case, the effective localised covariance matrix in state space is balanced and given by $\mathbf{K}_e^b \mathbf{B}_e^c (\mathbf{K}_e^b)^T$.

In practice, applying a Schur product directly would be prohibitively expensive. In variational assimilation, we can use the equivalent form

$$ \mathbf{L} \circ \tilde{\mathbf{X}}'_e \tilde{\mathbf{X}}'_e^T = \sum_{k=1; k \neq e}^{N_e} \Lambda_k \mathbf{L} \Lambda_k, \quad (11) $$

where $\Lambda_k$ is a diagonal matrix whose diagonal elements are the normalised perturbations for ensemble member $k$. The right-hand side of Equation (11) involves standard matrix–matrix multiplications and thus can be applied in the CG minimisation algorithm (which requires matrix–vector multiplications) using an operator form of the matrices. In NEMOVAR, localisation via $\mathbf{L}$ is done using a diffusion-based correlation operator. The cost of the localisation is proportional to the ensemble size, as

applying Equation (11) to a vector requires applying a diffusion-based correlation operator $N'_e$ times. While this is an expensive operation, it is not prohibitive. Furthermore, for NEMOVAR, we have made significant progress in reducing the cost of localisation by allowing the diffusion operator to be applied on a coarser grid, which is appropriate, since localisation length-scales are typically broad compared with the resolution of the model grid. In this study, however, we have not applied the coarse-grid feature.

The block structure of $\mathbf{L}$ can be formulated in a way to accommodate multivariate ensemble covariances and separate localisation length-scales for each variable, while maintaining the requirement that $\mathbf{L}$ be symmetric and positive definite. In this study, we consider a block-diagonal form of $\mathbf{L}$, which allows us to specify different localisation length-scales for each control variable (temperature, unbalanced salinity, unbalanced SSH) but
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>FIGURE 14 Ensemble departure of SLA in m. Panel (a) shows the ensemble departure for the control experiment ensda00_inf00. The difference of the control from experiments ensda02_inf08, ensda08_inf08, and ensda10_inf08 is shown in panels (b), (c), and (d), respectively. For these panels, the blue areas are where the departure of the experiment is less than the control ensda00_inf00 [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

eliminates the sample cross-correlations between these variables. In this case, the multivariate response of the analysis is determined entirely by the balance operator. Finally, localisation is only applied to the horizontal correlations in the experiments presented here. This was mainly a practical decision, due to the difficulty in determining suitable vertical localisation length-scales for different locations.

## 3.4 | Hybrid EnVar

In most practical cases, the ensemble is relatively small and thus the sample covariance is still noisy. A more robust background-error covariance matrix can be obtained by combining the ensemble-based and standard modelled covariance matrices linearly to form a hybrid B matrix, which in our case is defined by

$$ \mathbf{B}_{\ell}^{h} = (1 - \beta_{e}^{2}) \mathbf{B}^{m} + \beta_{e}^{2} \mathbf{B}_{\ell}^{e}, \quad (12) $$

where $\beta_e^2$ is a weight that controls how much the modelled covariance matrix $\mathbf{B}^m$ and the ensemble covariance matrix $\mathbf{B}_\ell^e$ contribute to the overall covariance matrix $\mathbf{B}_\ell^h$. Note that we have defined the weight such that its overall value for both components always adds up to 1. Therefore, if the variances in $\mathbf{B}^m$ and $\mathbf{B}_\ell^e$ are the same, then different values of $\beta_e^2$ will preserve that variance. We call the case where $\beta_e = 0$ standard 3DVar, the case where $\beta_e = 1$ 3DEnVar, and the case where $0 < \beta_e < 1$ hybrid-3DEnVar, following Lorenc and Jardak (2018).

## 3.5 | Ensemble inflation method

After the assimilation step, we apply inflation to reduce any underspread in the ensemble. Underspread is expected to occur for a number of reasons, principally sampling issues (see Sacher and Bartello, 2008), deficiencies in the data assimilation, and model error. Numerous inflation techniques have been proposed in the literature (see Houtekamer and Zhang, 2016, for a detailed review). For
<figure><img src="image_25.png" /><figcaption>FIGURE 15 As Figure 14, but ensemble spread SLA per m [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

our purposes, it was decided to use a return to prior spread (RTPS: Whitaker and Hamill, 2012) multiplicative inflation scheme similar to the RTPS scheme used in the MOGREPS-G atmospheric ensemble system (Jardak et al., 2018). This scheme is expected to compensate for underspread caused by the data assimilation and sampling issues, but not underspread caused by model error. The model perturbation scheme described in Section 2.3 is designed to deal with the effects of model error.

In an RTPS scheme, the ensemble analysis perturbation matrix $X'^a$ (i.e., the matrix of the perturbations after assimilation without self-exclusion) is replaced by a scaled version of itself:

$$X'^a \leftarrow YX'^a, \qquad (13)$$

where $Y$ is a diagonal matrix of multiplication factors.

This is a general formula for all multiplicative inflation methods. However, in an RTPS scheme, inflation is applied to the analysis perturbations, with $Y$ chosen so that the analysis spread is relaxed back towards the background spread. This means that

$$Y = \alpha \left( \frac{\mathrm{diag}[\mathbf{B}^e]^{\frac{1}{2}} - \mathrm{diag}[\mathbf{A}]^{\frac{1}{2}}}{\mathrm{diag}[\mathbf{A}]^{\frac{1}{2}}} \right) + \mathbf{I}, \qquad (14)$$

where $\mathbf{B}^e$ is the covariance of the full ensemble background ($\mathbf{B}^e = X'X'^T$: similar to $\mathbf{B}_e$ defined in Equations (8) and (9), but all ensemble members are included in its computation), and $\mathbf{A} = X'^aX'^aT$ is the full analysis ensemble covariance matrix.

Equation (14) contains a relaxation parameter $\alpha$, which controls how strongly the inflation pulls the analysis spread towards the background spread. From inspection of Equations (13) and (14), if $\alpha = 0$ the inflation will do nothing. Furthermore, taking the expectation of $YX'^a(YX'^a)^T$ and after some rearrangement, $E[YX'^a(YX'^a)^T] = (\mathrm{diag}[\mathbf{B}^e]/\mathrm{diag}[\mathbf{A}])\mathbf{A}$ when $\alpha = 1$; i.e., the spread of the analysis is inflated to have the same variance as the background ensemble (although the individual ensemble members will differ). Values of $\alpha$ between 0 and 1 will relax the analysis ensemble spread part way towards the background ensemble spread.

Care must be taken when inflating ocean variables that are correlated with each other, due to physical relationships. Inflation can potentially break these relationships and cause unwanted mixing and turbulence. This is a problem common in data assimilation, and is normally solved by applying inflation to the control variables, so that when they are converted back to full state variables using
<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /><figcaption>FIGURE 16 Ensemble departure of in situ SST (in °C). Panel (a) shows the ensemble departure for the control experiment ensda00_inf00. The difference of the control from experiments ensda02_inf08, ensda08_inf08, and ensda10_inf08 is shown in panels (b), (c), and (d), respectively. For these panels, any blue areas are where the departure of the experiment is less than the control ensda00_inf00 [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

the balance operator the physical relationships between parameters are respected.

A key aspect of our inflation implementation is that we do not inflate the analysis ensemble spread of the state variables directly. Instead, inflation is applied by modifying the analysis increments for the control variables for the ensemble members according to

$$ \Gamma \leftarrow \bar{\Gamma} + \sqrt{N_e - 1} ((\Upsilon - I)\mathbf{X}' + \Upsilon\Gamma'). \quad (15) $$

In Equation 15, $\Gamma$, the matrix of the analysis increments (where the increments of each ensemble member make up the columns of the matrix), is replaced by its own mean $\bar{\Gamma}$ plus a weighted combination of its own perturbations $\Gamma'$ and the background perturbations $\mathbf{X}'$. As the perturbations have, by definition, zero mean, it is easy to see that the mean of the increments, and thus of the analysis, is unchanged. A derivation of Equation (15) is given in the Appendix A. The inflated analysis increments for the control variables are then converted into analysis increments for the state variables temperature, salinity, SSH, and velocity, and the appropriate increment

is applied to the model for each ensemble member using IAU.

# 4 | EXPERIMENTAL SETUP AND ASSESSMENT METHODS

## 4.1 | Experimental setup

With our experiments, we aimed to address two main questions. First, we wished to assess how well the ocean ensemble system performs, that is, that we have a good ensemble (in some sense). Second, we wanted to see how sensitive the results are to the settings of our system, in particular the hybrid weight for the ensemble covariance, and to demonstrate the impact of inflation on the system. To address both questions, we ran a series of reanalysis experiments for the four-month period 1 September, 2018–31 December, 2018.

All of the ensemble runs assimilated the same observations with the same perturbations, as described in Section 2.3. They all started from a spun-up ensemble,
<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /><figcaption>FIGURE 17 As Figure 16, but for ensemble spread SST (in °C) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

which was generated by running the ensemble system described earlier, with the standard 3DVar, for six months starting on 1 March, 2018. On that date, the ocean and sea-ice fields from a previous run of the deterministic FOAM system were used to initialise all the ensemble members. Six months was long enough for the ensemble spread to stabilise in all variables, as shown in Figure 5.

The basic experiment details are listed in Table 3. Each experiment is named in the form ensda[x₀x₁]_inf[y₀y₁], where the hybrid weight on the ensemble component of the background-error covariance matrix is x₀.x₁ (denoted β²_e in Equation (12)) and y₀.y₁ is the multiplicative inflation factor α (see Equation (14)). All the experiments are identical apart from using different settings for β_e and α.

## 4.2 | Ensemble assessment methods

One of the basic assessment methods for a data assimilation system is to look at the innovations, that is, the difference between the observations and the model's short-range forecast before assimilation of the observations. The observations in the assessment are

the same ones assimilated, which were quality-control checked and have the previous day's observation-bias correction applied (see Section 2.1 for more details on the observations). For an ensemble system, not only do we want to have smaller innovations, we also need to have an ensemble spread that correctly the true forecast uncertainty represents. A method to assess this is given by Rodwell et al. (2016). In that article, the following equation is given as a diagnostic of the reliability of an ensemble:

$$Depar^2 = Bias^2 + EnsVar + ObsUnc^2 + Resid. \quad (16)$$

Each of the terms of this equation is an important diagnostic of an ensemble in its own right. The departure (Depar) is the standard deviation of the observations around the ensemble mean. The bias here refers to the mean of the observations minus the ensemble mean. The ensemble variance (EnsVar) is the variance of the ensemble members around their own mean. The observation uncertainty (ObsUnc) is the standard deviation of the observation errors. This value is often not well-known. For a reliable ensemble, the residual (Resid) should be statistically indistinguishable from zero.
<figure><img src="image_30.png" /><figcaption>FIGURE 18 Profile temperature (in °C) and salinity (in PSS) statistics of departure (solid line) and bias (dashed line) are shown in the top panels. The bottom panels show the ensemble spread. 95th percentile uncertainties are shown as shaded areas in the RMS departures. These are estimated using the bootstrapping method of Efron (1979) with 1000 resamples [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

It should be noted that the interpretation of the ensemble reliability depends on the observation errors, which are difficult to estimate. Here we assume the observation errors for altimeter, temperature and salinity profiles are the same as those used in *R* in the assimilation. For SST, we assume a uniform 0.3°C observation error—the values used in the assimilation for SST are artificially inflated to reduce the impact of observation-error correlations, so are not appropriate for use in this context. If these are not reasonable assumptions, this will affect the residual and our assessment of the ensemble reliability may not be correct.

The system used here is a coupled ocean/sea-ice model with assimilation of both ocean and sea-ice data. However, the ensemble aspects of the sea-ice have not yet received much attention: we have not implemented stochastic model perturbations in the sea-ice component of the model, and, while we do include a perturbed observation ensemble of 3DVars for the sea-ice concentration assimilation, we have not yet explored the use of hybrid-3DEnVar for the sea-ice component here.

# 5 | RESULTS

## 5.1 | Ensemble system performance: control run

We first look at the spatial distribution of the reliability for the control run ensda00_inf00. Results are shown in Figures 6–9 for the period October–December 2018, in which negative (positive) values in the residual plots indicate that the ensemble is overspread (underspread). The results for SLA (Figure 6) show that the spatial distribution of the ensemble spread reflects well the distribution of the uncertainty in the ensemble mean. The main areas where this is not true are the ACC, northern North Atlantic, and Agulhas regions, where the ensemble appears to be overspread. Note this assessment is dependent on our assumptions about the observation errors. Negative residuals in these regions could also result from assumed observation errors being too high (for this observation type). For SST (Figure 7), the residual is generally quite small, indicating that the ensemble spread is consistent with the
<figure><img src="image_31.png" /><figcaption>FIGURE 19 SSH increments in m in the unperturbed member on the first day of the experiments in the Gulf Stream region [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

uncertainty. The residual is negative in the Gulf Stream and ACC. This could result from an overspread ensemble in these areas. For the three-dimensional temperature of the ensemble in Figure 8, the residual is generally quite small, but it shows that the ensemble is generally overspread. The salinity residuals (Figure 9) follow a similar pattern to the temperature ones, although the salinity ensemble seems quite reliable in the Tropics except for the Indian Ocean, and particularly the Bay of Bengal, where the ensemble appears underspread. The land freshwater runoff has not been perturbed in the ensemble and this is something that could be included to increase the spread in the Bay of Bengal, the River Plate outflow region, and other large rivers.

Comparing the control with an experiment with hybrid DA ensda08_inf08 in Figure 10 shows that the hybrid

DA run has very good ensemble reliability for SLA (and other variables, not shown). In fact, the results for SLA (in Figure 10) are markedly improved over the control (Figure 6), particularly where the ensemble was previously apparently overspread in the ACC.

## 5.2 | Assessment of impact of changing hybrid weights and inflation

In Section 4 we introduced a series of experiments to test the impact of different hybrid weights and inflation factors. Table 3 summarises the settings of each experiment.

Figure 11 gives the evolution of the RMS and mean departures (the difference of the observations from the ensemble mean one-day forecast). For SLA, we see that
<figure><img src="image_32.png" /><figcaption>**FIGURE 20** Surface-temperature increments in °C in the unperturbed member on the first day of the experiments in the Gulf Stream region [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

all the hybrid-3DEnVar experiments and the pure 3DEn-Var experiment reduce the SLA errors markedly compared with the standard 3DVar experiment. These results are statistically significant: for example, when comparing five-day periods of hybrid-3DEnVar with the standard 3DVar experiment using the paired *t*-test, 0 or negligible *p*-values (probability of the two experiments being the same) are found. The differences between hybrid-3DEnVar experiments with hybrid weights of 0.5 compared with 0.8 are often less significant, with larger *p*-values sometimes close to 1 depending on the time period. Note that all the experiments start from the same initial conditions, but the hybrid-3DEnVar improves rapidly compared with the control, with most of the improvement established in the first month. An interesting result is that bias for SLA is close to zero except for

the control. This is only seen in the perturbed ensemble members and not in the unperturbed member (equivalent plot not shown, but compare Figures 6–10). This suggests that it is caused by an aspect of the ensemble perturbation strategy that is ameliorated by hybrid-3DEnVar. Most likely it is the SLA observation perturbations, but this is not confirmed.

For SST there is very little difference for most of the experiments, except for pure 3DEnVar, which has approximately 0.35 °C departure compared with around 0.30°C for the hybrid-3DEnVar and control experiments. Also, we can see the mean error increasing for pure 3DEnVar, whereas it remains close to 0 for the other experiments. Again, this difference is seen relatively fast, within a month of the start of the experiments. For profile temperature and salinity, both the 3DVar and 3DEnVar
<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /><figcaption>FIGURE 21 Temperature increments in °C (section 27°N–47°N at 68°W) in the unperturbed member on the first day of the experiments in a section in the Gulf Stream [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

give worse results than the hybrid-3DEnVar experiments, with 3DEnVar worse than the control. The spin up of the hybrid-3DEnVar for profile temperature is less rapid than SLA, taking between 1 and 2 months. An experiment with hybrid-3DEnVar but no inflation is included (*ensda08_inf00*). This shows that the impact of inflation on the errors with respect to observations is relatively small.

Looking at the ensemble spread time series in Figure 12, we see that hybrid DA causes a large fall in the ensemble spread, particularly for SLA. This is reduced the most at about 0.8 ensemble hybrid weight. We can also see that inflation brings the ensemble spread up compared with no inflation, and indications from assessing the ensemble reliability in Figure 10 suggest that without inflation the ensemble would be underspread with hybrid DA, but with inflation the ensemble spread is good in most locations. Again, a good deal of the spin-up of the experiments is complete in the first month and it may not be too far from complete after four months of running.

Now we compare time-averaged statistics for all experiments. We exclude the first month of the experiments from

this assessment, in order to remove most of the spin-up from the comparisons.

Figure 13 shows plots of the departure and ensemble spread for the period October–December 2018 over the global ocean for all the experiments as a function of the hybrid weight $\beta_e^2$. First, it should be noted that for most variables the ensemble mean departure is smaller than the departure for the unperturbed member. For SLA, for example, the ensemble mean departure is 4% smaller. This means that, even without hybrid DA, we achieve reductions in the short-term forecast error using the ensemble. This is the case for other variables too, except SST, where there is no reduction in error from using the ensemble mean. Perhaps the ensemble perturbation strategy introduces too much noise to SST in the ensemble members. The other key result shown in Figure 13 is that hybrid DA gives the optimal results. Particularly striking is the 20% reduction in SLA departure with hybrid DA compared with the control. The spread is generally reduced compared with the control, but so also is the departure. Going to pure 3DEnVar generally increases the departure.
<figure><img src="image_35.png" /><figcaption>FIGURE 22 Surface-temperature increments in °C in the unperturbed member on the first day of the experiments in the eastern tropical Pacific [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

A hybrid ensemble weight of either 0.5 or 0.8 gives the smallest departure, depending on the variable and region concerned. It is notable that the biggest changes in the departure come with either 3DVar or pure 3DEnVar, and it is relatively insensitive to the hybrid weights selected. Consistent with the time-series plots, we see the impact of inflation on increasing the ensemble spread. It has little impact on the departures.

A comparison of spatial maps of the time-average SLA departure (difference of the ensemble mean and observations before assimilation) for some of the experiments is shown in Figure 14. The departure is highest in active eddy regions like the Gulf Stream and the ACC. This is expected, as these are chaotic regions where the ocean forecast errors grow fastest. Comparing the control ensda00_inf00

with the experiments with hybrid or pure 3DEnVar shows reduced departures in all these regions. There is a slight increase in error that becomes apparent in parts of the Tropics and is seen particularly in the pure ensemble DA experiment ensda10_inf08.

The ensemble spread, shown in Figure 15, is also lower in the active eddy regions. This should not be a concern, since it reflects the reduced forecast errors due to an improved fit to the observations. Particularly for the pure ensemble DA experiment ensda10_inf08, there is an increase in spread seen in parts of the Tropics, particularly the tropical Atlantic and western Indian Ocean.

For SST (Figure 16), there is very little difference between the experiments, apart from an almost global increase in departure for the pure ensemble
<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /><figcaption>FIGURE 23 Temperature increments in °C (equatorial section 110°W–90°W) in the unperturbed member on the first day of the experiments in the eastern tropical Pacific [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

DA experiment ensda00_inf00. There are some signs of reduced errors along the Gulf Stream path, but there is no obvious signal in other eddying regions. Limited sampling, due to the short residency time, makes seeing such signals difficult in in situ data. The results for SST spread are somewhat similar to those seen with SLA (Figure 17). This suggests a not unsurprising association between the ensemble spread in SLA and the ensemble spread in SST associated with the position of eddies and fronts.

Another way to compare the experiments is to look at the impact on the ocean subsurface. Figure 18 shows the ensemble departure and spread for the experiments for temperature and salinity as a function of depth. The key results are consistent with those seen for surface variables. Hybrid-3DEnVar gives the lowest departures in temperature and salinity, with a hybrid ensemble weight of 0.5 or 0.8 giving the lowest ensemble mean departures, while 3DEnVar and standard 3DVar give higher departures. These differences are statistically significant when tested with the paired *t*-test. Using a paired *t*-test to compare the profile temperature RMS statistics of the hybrid-3DEnVar with the standard 3DVar experiment

(by looking at five-day periods), *p*-values are negligible. Comparing the hybrid-3DEnVar experiments with hybrid weights of 0.8 and 0.5, the differences in statistics are not so significant, with *p*-values sometimes close to 1 depending on the paired observations assessed. The spread of the ensemble is reduced, the higher the hybrid ensemble weight, but inflation acts to increase the spread somewhat without affecting the departure statistics significantly. One subtle difference is that 3DEnVar gives worse temperature departures than the 3DVar control near the surface and this peaks at about 200-m depth. Below about 400 m, the results for these experiments (3DEnVar and 3DVar) are the same. The DA experiments with some hybridisation still have lower departures than either 3DVar or 3DEnVar.

## 5.3 | Data assimilation increments

In order to try to understand the statistical results better, it is useful to look at the assimilation increments. In our experiments, the settings are the same, apart from the assimilation and inflation. On the first day of the runs, all
the inputs to the assimilation are identical, including the background fields used as the linearisation states for each member and the innovations (the observation and model perturbations are reproducible, so a run for the same cycle and ensemble member has the same perturbations). Any differences on the first day are due entirely to changes we have made to the assimilation itself. In Figure 19 we show the SSH increments on the first day in the Gulf Stream region. The impact of the tracks of the altimeter data is very clear, particularly in the standard 3DVar experiment ensda00\_inf00. The increments are “blobby” in form: for example, a large innovation at 69°W, 37°N gives an isotropic increment with a correlation length-scale equal to the first Rossby radius of deformation. There are also larger length-scale features in the increment, which come from the longer 400-km correlation length-scale. As the hybrid ensemble weight is increased, the (smooth) longer length-scale features reduce. The increments become less isotropic and the gaps in the tracks become less obvious as they are filled in by increments.

The surface-temperature increments on the first day (Figure 20) in the Gulf Stream region differ less for most of the experiments, particularly near the Gulf Stream. This is likely because the high density of the SST data means the background-error length-scales/structures affect the increments less than for other variables and at other depths. Nevertheless, for the 3DEnVar experiment (ensda10\_inf08) we do see much smaller increments away from the Gulf Stream. The ensemble spread (not shown) varies much more with location than the equivalent modelled background-error covariance (not shown). All indications are that the ensemble spread is reasonable, but it does not represent the total error that the assimilation is trying to correct; in particular, there are likely to be significant biases due to forcing and model errors that are not accounted for in the ensemble. One possibility may be that the ocean is too damped and is insufficiently sensitive to atmospheric forcing variability. This might be addressed in future with higher resolution models (or other model physics improvements). For now, the hybrid DA approach is effective. For this to work optimally, the error covariances for the modelled B may need retuning to ensure that sufficiently large corrections/increments are made in locations away from the boundary currents.

A 27°N–47°N temperature section (Figure 21) at 68°W shows a strong positive increment at the location of the Gulf Stream front, likely due to assimilation of SLA data. In the control experiment, the increments are vertically confined as a result of the multivariate balance relationships in the assimilation. As hybrid-3DEnVar is introduced, the increments are found to slope with depth, which reflects the slope in the Gulf Stream front. The increments are somewhat larger in value in the region of the front. As

the increments seem more physically plausible, they are expected to cause less shock to the model, which could explain the performance improvement of hybrid-3DEnVar (and 3DEnVar) in eddying locations.

A similar set of temperature increment plots is shown for the tropical Pacific in Figure 22. Some of the results are similar to those seen in the Gulf Stream. In addition, here the features in the SST increments exhibit markedly smaller scales as the hybrid ensemble weight increases. In this part of the ocean, there are strong tropical instability waves and these exhibit fine filaments in SST. The positioning of the filaments differs in each ensemble member, so they determine the main correlation patterns. This means these patterns will be seen in the 3DEnVar and hybrid-3DEnVar increments. Again, it is notable that the SST increments, particularly to the south of the Equator, are weak in the 3DEnVar experiment. A 110°W–90°W section along the Equator (Figure 23) shows that the increments are smaller in scale at all depths as the hybrid ensemble weight is increased. It is likely that the true background errors in the Tropics have a significant large-scale component, particularly away from the most variable regions. These are not being represented in the 3DEnVar experiment and this may be leading to worse SST performance in the Tropics.

Time-average surface-temperature increments for the unperturbed ensemble member in the Gulf Stream region are shown in Figure 24 for selected experiments. As hybrid DA is introduced and the hybrid ensemble weight is increased, the surface-temperature mean increments reduce, particularly away from the Gulf Stream. Time-mean increments such as these are often taken to identify locations where there is model bias. The limitations of such an assumption can be clearly seen here, as, with the same model, the mean increments are different depending on whether we use 3DVar, 3DEnVar, or hybrid DA.

At least for 3DEnVar, some biases in the SST remain uncorrected. For example, south of the Gulf Stream the mean SST innovations are seen to be significantly larger in 3DEnVar than in 3DVar or hybrid-3DEnVar. We also get an indication of how much work is being done by the assimilation by looking at the standard deviation of the increments. Ideally, we would like to achieve a good fit to the data with the minimum standard deviation of the increments, since applying strong increments to a model may result in initialisation shock. It is quite encouraging (in Figure 25) that the assimilation performs well when assessed against observations, as the hybrid ensemble weight increases even as the area of high standard deviation of the increments is increasingly restricted to locations close to the Gulf Stream front. However, for full ensemble DA, 3DEnVar, the increments may be too small,
<figure><img src="image_38.png" /></figure>

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

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

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

<figure><img src="image_42.png" /><figcaption>FIGURE 24 Time mean (October 2018) of the surface-temperature increments in °C for the unperturbed member near the Gulf Stream [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

away from the strongly eddying regions, to keep the model
close to the observed SST.

6 | SUMMARY AND DISCUSSION

This article describes the new Met Office ocean ensemble
system. It uses an ensemble of atmospheric forcing from
MOGREPS-G to force an ensemble of 3DVars in the ocean
and sea-ice with perturbed observation locations and val-
ues following a similar approach to that of Zuo et al. (2018).
The system includes stochastic model perturbations based
on the schemes developed by Storto and Andriopoulos
(2021). Given our best estimates of the observation errors,
we have found the ensemble to be reliable, that is, it has

a spread that is representative of the uncertainty in the
ensemble mean.

We focussed on a series of four-month experiments to
test the use of ensemble information in the specification of
the background-error covariances in an ensemble of data
assimilations in the ocean, and also the impact of includ-
ing an RTPS inflation scheme. Each run started from a
previously spun up ensemble and assimilated the same
observations. The experiments included standard 3DVar,
hybrid-3DEnVar with different hybrid weights, and 3DEn-
Var. In each experiment, all the ensemble members used
the same DA approach with the same settings. We found
that a hybrid ensemble weight of 0.8 and a RTPS inflation
factor of 0.8 together were able to improve the ensemble
reliability, as well as reducing the errors in the ensemble
<figure><img src="image_43.png" /><figcaption>FIG URE 25 Time-averaged standard deviation (October 2018) of the surface temperature increments in °C for the unperturbed member near the Gulf Stream [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

mean. The improvements (reduced ensemble spread and ensemble mean error) are particularly noticeable in the Gulf Stream and ACC, where the ensemble was previously overdispersive.

The departure (the difference between the ensemble mean and the observations) of the ensemble system and the innovations for the unperturbed member are reduced significantly with hybrid-3DEnVar (demonstrated by paired t-tests and bootstrap confidence limits). In particular, a hybrid ensemble weight of 0.8 (or perhaps 0.5 for some variables) produces a reduced departure compared with a 3DVar ensemble run or a 3DEnVar ensemble. The highlight is the SLA RMS error, which is 20% lower in hybrid-3DEnVar compared with 3DVar. Since it is expensive to test many different values, it is fortunate that the

performance is found to be relatively independent of the exact hybrid ensemble weight used.

We examined the increments from our experiments to attempt to understand why hybrid-3DEnVar works so well. Two main elements seem to be important: first, that increments are in a form that is more consistent with the model, particularly in terms of the vertical and horizontal structure of the increments, and second, that hybrid-3DEnVar avoids adding too large increments in more quiescent regions where the background error is lower than assumed in standard 3DVar. A 3DEnVar run was also tested and this did not perform as well as the hybrid-3DEnVar. This seemed to be because it was unable to correct biases in regions where the ensemble spread was low.
The main drawback of hybrid-3DEnVar is the additional computational expense (at least a factor of two more expensive). The assimilation step must read in the ensemble perturbations, which needs more memory. The localisation of the ensemble also requires a costly application of a correlation operator. In addition, NEMOVAR must still perform all the same calculations as the standard 3DVar. Since it is so effective at improving the accuracy of the analysis and short-range forecasts, future work will explore whether the efficiency can be improved while still retaining as many of the benefits as possible. Options to explore are employing a cheaper but more sensitive multivariate localisation, use of a lower resolution grid for the ensemble localisation, or performing some aspects of the computations at lower precision. It may be that, even after exploring these options, operationally hybrid-3DEnVar is still found to be too expensive for initial operational implementation, so it is therefore worth restating that the ensemble of 3DVars performs well and is already an improvement over the current operational system.

The ultimate aim of this work is to implement improvements to an operational global coupled NWP atmosphere-land-ocean-ice ensemble system. The main differences between the frameworks used to develop the ocean and sea-ice aspects described in this article are that they have been made in an uncoupled system, and with a one-day assimilation time window as opposed to the 6-hr window of the coupled NWP system being prepared for operational implementation. In order to be made ready for operational use, these two aspects will be addressed. The impact of the improved ocean component of the coupled ensemble on the atmospheric ensemble will be assessed in future work.

## AUTHOR CONTRIBUTIONS

**Daniel J. Lea:** conceptualization; formal analysis; investigation; methodology; project administration; software; validation; visualization; writing – original draft; writing – review and editing. **James While:** conceptualization; formal analysis; investigation; methodology; software; validation; visualization; writing – original draft; writing – review and editing. **Matthew J. Martin:** conceptualization; formal analysis; investigation; methodology; software; supervision; validation; visualization; writing – original draft; writing – review and editing. **Anthony Weaver:** conceptualization; methodology; software; writing – original draft; writing – review and editing. **Andrea Storto:** conceptualization; methodology; software; writing – original draft; writing – review and editing. **Marcin Chrust:** conceptualization; methodology; software; writing – review and editing.

## ACKNOWLEDGEMENTS

The contributions by D.J.L., M.J.M., and J.W. to this work were conducted through the Weather and Climate Science for Service Partnership (WCSSP) India, a collaborative initiative between the Met Office, supported by the UK Government 2019 Newton Fund, and the Indian Ministry of Earth Sciences (MoES). The contribution of A.W. to this work has been funded by the European Union through the Copernicus Climate Change Service (C3S) implemented by ECMWF. This is part of the work contributing to contract C3S_321b: "Technical preparations for C3S Seasonal Initialization and Global Reanalysis: Enabling an Ensemble of Data Assimilations for the Ocean".

## ORCID

Daniel J. Lea https://orcid.org/0000-0003-1736-628X

James While https://orcid.org/0000-0001-9870-8725

Matthew J. Martin https://orcid.org/0000-0003-0293-3106

Anthony Weaver https://orcid.org/0000-0002-8952-338X

Andrea Storto https://orcid.org/0000-0003-3856-8905

## REFERENCES

Balmaseda, M.A., Mogensen, K. and Weaver, A.T. (2013) Evaluation of the ECMWF ocean reanalysis system ORAS4. *Quarterly Journal of the Royal Meteorological Society*, 139, 1132–1161. https://doi.org/10.1002/qj.2063.

Bannister, R.N. (2017) A review of operational methods of variational and ensemble-variational data assimilation. *Quarterly Journal of the Royal Meteorological Society*, 143, 607–633. https://doi.org/10.1002/qj.2982.

Barton, N., Metzger, E.J., Reynolds, C.A., Ruston, B., Rowley, C., Smedstad, O.M., Ridout, J.A., Wallcraft, A., Frolov, S., Hogan, P., Janiga, M.A., Shriver, J., McLay, J., Thoppil, P., Huang, A., Crawford, W., Whitcomb, T., Bishop, C., Zamudio, L. and Phelps, M. (2020) The Navy’s Earth System Prediction Capability: a new global coupled atmosphere-ocean-sea ice prediction system designed for daily to subseasonal forecasting. *Earth and Space Science*, 8, e2020EA001199. https://doi.org/10.1029/2020EA001199.

Bell, M.J., Martin, M.J. and Nichols, N.K. (2004) Assimilation of data into an ocean model with systematic errors near the equator. *Quarterly Journal of the Royal Meteorological Society*, 130, 873–893. https://doi.org/10.1256/qj.02.109.

Blockley, E.W., Martin, M.J., McLaren, A.J., Ryan, A.G., Waters, J., Guiavarc'h, C., Lea, D.J., Mirouze, I., Peterson, K.A., Sellar, A., Storkey, D. and While, J. (2014) Recent development of the Met Office operational ocean forecasting system: an overview and assessment of the new Global FOAM forecasts. *Geoscientific Model Development*, 6, 6219–6278. https://doi.org/10.5194/gmdd-6-6219-2013.

Bloom, S.C., Takacs, L.L., da Silva, A.M. and Ledvina, D. (1996) Data assimilation using incremental analysis updates. *Monthly Weather Review*, 124, 1256–1271. https://doi.org/10.1175/1520-0493(1996)124<1256:DAUIAU>2.0.CO;2.
Bowler, N.E., Arribas, A., Mylne, K.R., Robertson, K.B. and Beare, S.E. (2008) The MOGREPS short-range ensemble prediction system. *Quarterly Journal of the Royal Meteorological Society*, 134, 703–722. https://doi.org/10.1002/qj.234.

Bowler, N.E., Clayton, A.M., Jardak, M., Lee, E., Lorenc, A., Pring, C.P.S.R., Wlasak, M.A., Barker, D.M., Inverarity, G.W. and Swinbank, R. (2017) Inflation and localization tests in the development of an ensemble of 4D-ensemble variational assimilations. *Quarterly Journal of the Royal Meteorological Society*, 143, 1280–1302. https://doi.org/10.1002/qj.3004.

Browne, P.A., de Rosnay, P., Zuo, H., Bennett, A. and Dawson, A. (2019) Weakly coupled ocean–atmosphere data assimilation in the ECMWF NWP System. *Remote Sensing*, 11, 234. https://doi.org/10.3390/rs11030234.

Buizza, R., Milleer, M. and Palmer, T.N. (1999) Stochastic representation of model uncertainties in the ECMWF ensemble prediction system. *Quarterly Journal of the Royal Meteorological Society*, 125, 2887–2908. https://doi.org/10.1002/qj.49712556006.

Clayton, A.M., Lorenc, A.C. and Barker, D.M. (2013) Operational implementation of a hybrid ensemble/4D-Var global data assimilation system at the Met Office. *Quarterly Journal of the Royal Meteorological Society*, 139, 1445–1461. https://doi.org/10.1002/qj.2054.

Efron, B. (1979) Bootstrap methods: Another look at the jackknife. *The Annals of Statistics*, 7, 1–26.

Gürol, S., Weaver, A.T., Moore, A.M., Piacentini, A., Arango, H.G. and Gratton, S. (2014) B-preconditioned minimization algorithms for variational data assimilation with the dual formulation. *Quarterly Journal of the Royal Meteorological Society*, 140, 539–556. https://doi.org/10.1002/qj.2150.

Hamill, T.M., Whitaker, J.S. and Snyder, C. (2001) Distance-dependent filtering of background-error covariance estimates in an ensemble Kalman filter. *Monthly Weather Review*, 129, 2776–2790.

Houtekamer, P.L. and Mitchell, H.L. (2001) A sequential ensemble Kalman filter for atmospheric data assimilation. *Monthly Weather Review*, 129, 123–137. https://doi.org/10.1175/1520-0493(2001)129<0123:ASEKFF>2.0.CO;2.

Houtekamer, P.L. and Zhang, F. (2016) Review of the ensemble Kalman filter for atmospheric data assimilation. *Monthly Weather Review*, 144, 4489–4532. https://doi.org/10.1175/MWR-D-15-0440.1.

Ingleby, B. and Huddleston, M. (2007) Quality control of ocean temperature and salinity profiles – historical and real-time data. *Journal of Marine Systems*, 65, 158–175. https://doi.org/10.1016/j.jmarsys.2005.11.019.

Janjic, T., Bormann, N., Bocquet, M., Carton, J., Cohn, S., Dance, S., Losa, S., Nichols, N., Potthast, R., Waller, J. and Weston, P. (2018) On the representation error in data assimilation. *Quarterly Journal of the Royal Meteorological Society*, 144, 1257–1278. https://doi.org/10.1002/qj.3130.

Kepert, J.D. (2009) Covariance localisation and balance in an ensemble Kalman filter. *Quarterly Journal of the Royal Meteorological Society*, 135, 1157–1176. https://doi.org/10.1002/qj.443.

Kjellsson, J. and Zanna, L. (2017) The impact of horizontal resolution on energy transfers in global ocean models. *Fluids*, 2, 45. https://doi.org/10.3390/fluids2030045.

Lea, D.J., Drecourt, J.P., Haines, A.K. and Martin, M.J. (2008) Ocean altimeter assimilation with observational and model-bias

correction. *Quarterly Journal of the Royal Meteorological Society*, 134, 1761–1774. https://doi.org/10.1002/qj.320.

Lea, D.J., Mirouze, I., Martin, M.J., King, R.R., Hines, A., Walters, D. and Thurlow, M. (2015) Assessing a new coupled data assimilation system based on the Met Office coupled atmosphere–land–ocean–sea ice model. *Monthly Weather Review*, 143, 4678–4694. https://doi.org/10.1175/MWR-D-15-0174.1.

van Leeuwen, P.J. (2020) A consistent interpretation of the stochastic version of the ensemble Kalman filter. *Quarterly Journal of the Royal Meteorological Society*, 146, 2815–2825. https://doi.org/10.1002/qj.3819.

Leutbecher, M., Lock, S.-J., Ollinaho, P., Lang, S.T.K., Balsamo, G., Bechtold, P., Bonavita, M., Christensen, H.M., Diamantakis, M., Dutra, E., English, S., Fisher, M., Forbes, R.M., Goddard, J., Haiden, T., Hogan, R.J., Juricke, S., Lawrence, H., MacLeod, D., Magnusson, L., Malardel, S., Massart, S., Sandu, I., Smolarkiewicz, P.K., Subramanian, A., Vitart, F., Wedi, N. and Weisheimer, A. (2017) Stochastic representations of model uncertainties at ECMWF: state of the art and future vision. *Quarterly Journal of the Royal Meteorological Society*, 143, 2315–2339. URL: https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.3094.

Lorenc, A.C. and Jardak, M. (2018) A comparison of hybrid variational data assimilation methods for global NWP. *Quarterly Journal of the Royal Meteorological Society*, 144, 2748–2760.

MacLachlan, C., Arribas, A., Peterson, K.A., Maidens, A., Fereday, D., Scaife, A.A., Gordon, M., Vellinga, M., Williams, A., Comer, R.E., Camp, J., Xavier, P. and Madec, G. (2015) Global seasonal forecast system version 5 (GloSea5): a high-resolution seasonal forecast system. *Quarterly Journal of the Royal Meteorological Society*, 141, 1072–1084. https://doi.org/10.1002/qj.2396.

Ménétrier, B. and Auligné, T. (2015) Optimized localization and hybridization to filter ensemble-based covariances. *Monthly Weather Review*, 143, 3931–3947. https://doi.org/10.1175/MWR-D-15-0057.1.

Mirouze, I., Blockley, E.W., Lea, D.J., Martin, M.J. and Bell, M.J. (2016) A multiple length scale correlation operator for ocean data assimilation. *Tellus A: Dynamic Meteorology and Oceanography*, 68. https://doi.org/10.3402/tellusa.v68.29744.

Mirouze, I. and Storto, A. (2019) Generating atmospheric forcing perturbations for an ocean data assimilation ensemble. *Tellus A: Dynamic Meteorology and Oceanography*, 71, 1624459. https://doi.org/10.1080/16000870.2019.1624459.

Ollinaho, P., Lock, S.-J., Leutbecher, M., Bechtold, P., Beljaars, A., Bozzo, A., Forbes, R.M., Haiden, T., Hogan, R.J. and Sandu, I. (2017) Towards process-level representation of model uncertainties: stochastically perturbed parametrizations in the ECMWF ensemble. *Quarterly Journal of the Royal Meteorological Society*, 143, 408–422. https://doi.org/10.1002/qj.2931.

Penny, S.G., Behringer, D., Carton, J. and Kalnay, E. (2015) A hybrid global ocean data assimilation system at NCEP. *Monthly Weather Review*, 143, 4660–4677. https://doi.org/10.1175/MWR-D-14-00376.1.

Ridley, J.K., Blockley, E.W., Keen, A.B., Rae, J.G.L., West, A.E. and Schroeder, D. (2018) The sea ice model component of HadGEM3-GC3.1. *Geoscience Model Development*, 11, 713–723. https://doi.org/10.5194/gmd-11-713-2018.
Rodwell, M.J., Lang, S.T.K., Ingleby, B.N., Bormann, N., Hólm, E., Rabier, F., Richardson, D.S. and Yamaguchi, M. (2016) Reliability in ensemble data assimilation. *Quarterly Journal of the Royal Meteorological Society*, 142, 443–454. https://doi.org/10.1002/qj.2663.

Sacher, W. and Bartello, P. (2008) Sampling errors in ensemble Kalman filtering. Part I: Theory. *Monthly Weather Review*, 136, 3035–3049. https://doi.org/10.1175/2007MWR2323.1.

Sakov, P., Counillon, F., Bertino, L., Lisæter, K., Oke, P. and Korablev, A. (2012) TOPAZ4: an ocean-sea ice data assimilation system for the North Atlantic and Arctic. *Ocean Sci.*, 8, 633–656. https://doi.org/10.5194/os-8-633-2012.

Sanchez, C., Williams, K.D. and Collins, M. (2016) Improved stochastic physics schemes for global weather and climate models. *Quarterly Journal of the Royal Meteorological Society*, 142, 147–159. https://doi.org/10.1002/qj.2640.

Schiller, A., Brassington, G., Oke, P., Cahill, M., Divakaran, P., Entel, M., Freeman, J., Griffin, D., Herzfeld, M., Hoeke, R., Huang, X., Jones, E., King, E., Parker, B., Pitman, T., Rosebrock, U., Sweeney, J., Taylor, A., Thatcher, M., Woodham, R. and Zhong, A. (2020) Bluelink ocean forecasting Australia: 15 years of operational ocean service delivery with societal, economic and environmental benefits. *Journal of Operational Oceanography*, 13, 1–18. https://doi.org/10.1080/1755876X.2019.1685834.

Storkey, D., Blaker, A.T., Mathiot, P., Megann, A., Aksenov, Y., Blockley, E.W., Calvert, D., Graham, T., Hewitt, H.T., Hyder, P., Kuhlbrodt, T., Rae, J.G.L. and Sinha, B. (2018) UK Global Ocean GO6 and GO7: a traceable hierarchy of model resolutions. *Geoscientific Model Development*, 11, 3187–3213. https://doi.org/10.5194/gmd-11-3187-2018.

Storkey, D., Blockley, E.W., Furner, R., Guiavarc'h, C., Lea, D., Martin, M.J., Barciela, R.M., Hines, A., Hyder, P. and Siddorn, J.R. (2010) Forecasting the ocean state using NEMO: The new FOAM system. *Journal of Operational Oceanography*, 3, 3–15. https://doi.org/10.1080/1755876X.2010.11020109.

Storto, A. and Andriopoulos, P. (2021) A new stochastic oceanphysics package and its application to hybrid-covariance data assimilation. *Quarterly Journal of the Royal Meteorological Society*, 147, 1691–1725. https://doi.org/10.1002/qj.3990.

Tennant, W.J., Shutts, G.J., Arribas, A. and Thompson, S.A. (2011) Using a stochastic kinetic energy backscatter scheme to improve MOGREPS probabilistic forecast skill. *Monthly Weather Review*, 139, 1190–1206. https://doi.org/10.1175/2010MWR3430.1.

Thoppil, P.G., Frolov, S., Rowley, C.D., Reynolds, C.A., Jacobs, G.A., Metzger, E.J., Hogan, P.J., Barton, N., Wallcraft, A.J., Smedstad, O.M. and Shriver, J.F. (2021) Ensemble forecasting greatly expands the prediction horizon for ocean mesoscale variability. *Communications Earth and Environment*, 2, 89. https://doi.org/10.1038/s43247-021-00151-5.

Waters, J., Lea, D., Martin, M., Mirouze, I., Weaver, A. and While, J. (2015) Implementing a variational data assimilation system in an operational 1/4 degree global ocean model. *Quarterly Journal of the Royal Meteorological Society*, 141, 333–349. https://doi.org/10.1002/qj.2388.

Weaver, A.T., Crust, M., Ménétrier, B. and Piacentini, A. (2021) An evaluation of methods for normalizing diffusion-based covariance operators in variational data assimilation. *Quarterly Journal of the Royal Meteorological Society*, 147, 289–320. https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.3918.

Weaver, A.T., Crust, M., Ménétrier, B., Piacentini, A., Tshimanga, J., Yang, Y., Gurol, S. and Zuo, H. (2018) *Using ensemble-estimated background error variances and correlation scales in the NEMOVAR system*. Toulouse, France: CERFACS. Technical Report TR/PA/18/15.

Weaver, A.T., Deltel, C., Machu, E., Ricci, S. and Daget, N. (2005) A multivariate balance operator for variational ocean data assimilation. *Quarterly Journal of the Royal Meteorological Society*, 131, 3605–3625. https://doi.org/10.1256/qj.05.119.

Weaver, A.T., Tshimanga, J. and Piacentini, A. (2016) Correlation operators based on an implicitly formulated diffusion equation solved with the Chebyshev iteration. *Quarterly Journal of the Royal Meteorological Society*, 142, 455–471. https://doi.org/10.1002/qj.2664.

While, J. and Martin, M.J. (2019) Variational bias correction of satellite sea-surface temperature data incorporating observations of the bias. *Quarterly Journal of the Royal Meteorological Society*, 145, 2733–2754. https://doi.org/10.1002/qj.3590.

Whitaker, J.S. and Hamill, T.M. (2012) Evaluating methods to account for system errors in ensemble data assimilation. *Monthly Weather Review*, 140, 3078–3089. https://doi.org/10.1175/MWR-D-11-00276.1.

Zuo, H., Alonso-Balmaseda, M., de Boisseson, E., Hirahara, S., Crust, M. and de Rosnay, P. (2017) *A generic ensemble generation scheme for data assimilation and ocean analysis*. Reading, UK: ECMWF. ECMWF Technical Memorandum No. 795.

Zuo, H., Alonso-Balmaseda, M., Mogensen, K. and Tietsche, S. (2018) *OCEAN5: The ECMWF ocean reanalysis system and its real-time analysis component*. ECMWF. ECMWF Technical Memorandum No. 823.

**How to cite this article:** Lea, D.J., While, J., Martin, M.J., Weaver, A., Storto, A. & Crust, M. (2022) A new global ocean ensemble system at the Met Office: Assessing the impact of hybrid data assimilation and inflation settings. *Quarterly Journal of the Royal Meteorological Society*, 148(745), 1996–2030. Available from: https://doi.org/10.1002/qj.4292

## APPENDIX A. INFLATING AN ENSEMBLE BY MODIFYING THE ANALYSIS INCREMENTS

It is convenient to combine ensemble inflation with the analysis increments, so that they can be applied to the model together. The key requirement is to apply the multiplicative inflation as an addition operation. Here we give the derivation used to modify the increments applied to each ensemble member to achieve this.

We define X and Γ as the matrices whose columns are the state vectors of the individual ensemble members and their increments, respectively. It is important to note that these matrices are non-exclusive and contain all of the ensemble members. Consequently they do not have the
subscript 'e' seen in similar matrices in Section 3.2. With these definitions,

$$ \mathbf{X}^a = \mathbf{X} + \mathbf{\Gamma}, \qquad (A1) $$

where the superscript “a” denotes analysis. We can split $\mathbf{X}^a$ into the sum of a mean matrix (denoted by an overbar) and a scaled perturbation matrix (as in the main text, this is denoted by a prime). A mean matrix has the same shape as the perturbation matrix, but has identical columns that are equal to the ensemble mean. With this decomposition, we obtain

$$ \overline{\mathbf{X}}^a + \sqrt{N_e - 1} \mathbf{X}'^a = (\overline{\mathbf{X}} + \overline{\mathbf{\Gamma}}) + \sqrt{N_e - 1} (\mathbf{X} + \mathbf{\Gamma})', \quad (A2) $$

where the $\sqrt{N_e - 1}$ scaling factor is due to our definition of the perturbation matrix (see Equations (8) and (9)).

Using Equation (13) to replace the perturbations by their inflated values, multiplicative inflation applied to the analysis ensemble is given by

$$ \overline{\mathbf{X}}^a + \sqrt{N_e - 1} \mathbf{Y} \mathbf{X}'^a = (\overline{\mathbf{X}} + \overline{\mathbf{\Gamma}}) + \sqrt{N_e - 1} \mathbf{Y} (\mathbf{X} + \mathbf{\Gamma})'. \quad (A3) $$

Since all averages are over the same number of ensemble members, we can split the averaged and perturbation matrices on the right-hand side, so that

$$ \overline{\mathbf{X}}^a + \sqrt{N_e - 1} \mathbf{Y} \mathbf{X}'^a = \overline{\mathbf{X}} + \overline{\mathbf{\Gamma}} + \sqrt{N_e - 1} \mathbf{Y} (\mathbf{X}' + \mathbf{\Gamma}')'. \quad (A4) $$

Adding and subtracting $\sqrt{N_e - 1} \mathbf{X}'$ from the right-hand side and noting that $\mathbf{X} = \overline{\mathbf{X}} + \sqrt{N_e - 1} \mathbf{X}'$, we obtain

$$ \overline{\mathbf{X}}^a + \sqrt{N_e - 1} \mathbf{Y} \mathbf{X}'^a = \mathbf{X} + \sqrt{N_e - 1} (\mathbf{Y} - \mathbf{I}) \mathbf{X}' + \mathbf{Y} \mathbf{\Gamma}' + \overline{\mathbf{\Gamma}}. \quad (A5) $$

Then, by inspection, multiplicative inflation can be applied by replacing the increments according to

$$ \mathbf{\Gamma} \leftarrow \overline{\mathbf{\Gamma}} + \sqrt{N_e - 1} (\mathbf{Y} - \mathbf{I}) \mathbf{X}' + \mathbf{Y} \mathbf{\Gamma}'. $$