Comparison of 4-dimensional variational and ensemble optimal interpolation data assimilation systems using a Regional Ocean Modeling System (v3.4) configuration of the eddy-dominated East Australian Current system

Ocean models must be regularly updated through the assimilation of observations (data assimilation) in order to correctly represent the timing and locations of eddies. Since initial conditions play an important role in the quality of short-term ocean forecasts, an effective data assimilation scheme to produce accurate state estimates is key to improving prediction. Western boundary current regions, such as the East Australia Current system, are highly variable regions, making them particularly challenging to model and predict. This study assesses the performance of two ocean data assimilation systems in the East Australian Current system over a 2-year period. We compare the time-dependent 4-dimensional variational (4D-Var) data assimilation system with the more computationally efficient, time-independent ensemble optimal interpolation (EnOI) system, across a common modelling and observational framework. Both systems assimilate the same observations: satellite-derived sea surface height, sea surface temperature, vertical profiles of temperature and salinity (from Argo floats), and temperature profiles from expendable bathythermographs. We analyse both systems' performance against independent data that are withheld, allowing a thorough analysis of system performance. The 4D-Var system is 25 times more expensive but outperforms the EnOI system against both assimilated and independent observations at the surface and subsurface. For forecast horizons of 5 d, root-mean-squared forecast errors are 20 %–60 % higher for the EnOI system compared to the 4D-Var system. The 4D-Var system, which assimilates observations over 5 d windows, provides a smoother transition from the end of the forecast to the subsequent analysis field. The EnOI system displays elevated low-frequency (>1 d) surface-intensified variability in temperature and elevated kinetic energy at length scales less than 100 km at the beginning of the forecast windows. The 4D-Var system displays elevated energy in the near-inertial range throughout the water column, with the wavenumber kinetic energy spectra remaining unchanged upon assimilation. Overall, this comparison shows quantitatively that the 4D-Var system results in improved predictability as the analysis provides a smoother and more dynamically balanced fit between the observations and the model's time-evolving flow. This advocates the use of advanced, time-dependent data assimilation methods, particularly for highly variable oceanic regions, and motivates future work into further improving data assimilation schemes.

1 Introduction

Data assimilation (DA), the combination of numerical modelling and observations, is essential to produce accurate forecasts of the atmosphere or ocean circulation. The goal of any DA scheme is to combine observations and a numerical model such that the result is a better estimate of the ocean circulation than either alone. Observations provide sparse data points, while the model provides context. Since initial conditions play an important role in forecast quality, accurate and dynamically consistent state estimates are key to improving prediction. This study focuses on the comparison of two DA techniques applied to forecasting the ocean mesoscale circulation in a highly dynamic oceanic region.

Mesoscale eddies exist throughout the global ocean and contain more than half of the kinetic energy of the ocean circulation. Western boundary current (WBC) regions are hotspots of high eddy variability as eddies emerge due to instabilities in the strong boundary current flow. The high mesoscale eddy variability (Stammer, 1997; Mata et al., 2000) and the complexities of eddy shedding processes and evolution (Mata et al., 2006; Bull et al., 2017) make WBCs challenging to model and predict (Feron, 1995; Imawaki et al., 2013; Roughan et al., 2017). Due to the chaotic nature of the mesoscale circulation, ocean models must be regularly updated through the assimilation of observations in order to correctly represent the timing and locations of eddies (e.g. Kerry et al., 2016; Li and Roughan, 2023), and accurate forecasts of eddies as they shed, evolve, and interact in WBC regions are lacking.

The East Australian Current (EAC), the WBC of the South Pacific subtropical gyre (Fig. 1a), and its associated eddies dominate the circulation along the southeastern coast of Australia. The southward-flowing current is most coherent off 27° S (Sloyan et al., 2016) and intensifies at around 31° S (Kerry and Roughan, 2020a). The current typically separates from the coast between 31 and 32.5° S (Cetina Heredia et al., 2014) and turns eastward to form the EAC eastern extension, shedding large warm-core eddies in the Tasman Sea (Oke and Middleton, 2000; Cetina Heredia et al., 2014; Oke et al., 2019). In the EAC, eddies can directly influence shelf circulation (Schaeffer et al., 2014; Schaeffer and Roughan, 2015; Malan et al., 2023) and often intensify as the jet separates from the coast. After shedding, eddies propagate and evolve (Pilo et al., 2015b, a) and can display a complex vertical structure including tilting and stacking (Oke and Griffin, 2011; Macdonald et al., 2013; Roughan et al., 2017; Pilo et al., 2018). As such, the EAC is a challenging region to predict and provides an ideal test bed for comparison of DA methods.

There are various DA techniques, by which a model estimate of the ocean state can be combined with ocean observations, that vary in complexity. Simpler, computationally efficient, time-independent methods such as 3-dimensional variational data assimilation (3D-Var) and ensemble optimal interpolation (EnOI) centre the observations and model on a single time and are capable of resolving slowly evolving flows governed by simple balance relationships at synoptic scales. These methods have provided useful state estimates and predictions. For example, the European Centre for Medium-Range Weather Forecasts uses 3D-Var to produce initial conditions for its coupled ocean–atmosphere modelling system (Mogensen et al., 2012), and EnOI was effectively employed in Australia's Bluelink Ocean Data Assimilation System (Oke et al., 2008a). In Oke et al. (2010) a case was presented for the use of EnOI, weighing up the predictive skill against its computational efficiency. Specifically, EnOI is highly computationally efficient as it does not represent the errors of the day; rather it assumes that the background error covariances are well represented by a stationary or seasonally varying ensemble. More recent work has shown that combining flow-dependent background error covariances (from an ensemble of model solutions) with a static ensemble achieves improved predictive skill (Brassington et al., 2023).

With increasing computational capacity and the pursuit of more accurate weather and ocean forecasts over the last 2 decades, a shift has been made to more advanced, time-dependent DA techniques (Edwards et al., 2015; Moore et al., 2019). Advanced DA methods make use of the time-variable dynamics of the model, allowing the observations to be assimilated over a time interval given the temporal evolution of the circulation. In the atmosphere, these methods have provided considerable improvement compared to the earlier, time-independent DA techniques, particularly for forecasts (e.g. Lorenc and Rawlins, 2005; Brousseau et al., 2012) and for highly intermittent flows with irregularly sampled observations (e.g. Xu, 2013). Indeed, the two techniques that are the most promising in numerical weather prediction (NWP) are 4-dimensional variational data assimilation (4D-Var) and the ensemble Kalman filter (EnKF), and ocean DA is following suit (Moore et al., 2019).

In 4D-Var the model and observations are combined using subsequent iterations of the tangent linear and adjoint models to compute increments in the forecast model (initial conditions, boundary conditions, and surface forcing) such that the difference between the new model solution and the observations is minimised over a time window (Moore et al., 2004). With 4D-Var, a continual and full estimate of the ocean over the assimilation window is created. This is ideal for both accuracy and timeliness of current state estimates and future predictions, as a continuous field evolves by the nonlinear primitive equations. The Kalman filter (KF) can be formally posed in the same way as 4D-Var (Lorenc, 1986) and in practice uses an ensemble of perturbed model simulations to approximate the model error covariances and their temporal evolution, and the ensemble mean is considered the best estimate of the state of the system (Evensen, 2002). An advantage of generating an ensemble of forecasts is that probabilistic forecasts can be derived from the ensemble spread.

Indeed, with the shift to more advanced DA techniques in ocean forecasting, it is important to quantify the improvements gained. Here we use a Regional Ocean Modeling System (ROMS) configuration of a dynamic WBC (the EAC) to compare two DA methods in a quantifiable manner. We compare the time-independent DA technique (EnOI) with the time-dependent technique (4D-Var) using the same numerical model configuration and suite of observations. We quantify the differences in predictive skill achieved by the two systems against assimilated and independent observations at the surface and subsurface. We focus our analysis on the performance of the short-range (5 d) forecasts. After presenting the experiments (Sect. 2), we begin by comparing forecast performance against assimilated observations (Sect. 3.1). Then we employ a suite of independent observations to assess the forecast skill of the two systems (Sect. 3.2). The model energetics (Sect. 4.2) and the temporal and spatial scales of variability (Sect. 4.3) are then compared to understand what may drive differences in predictive skill. Finally we summarise and discuss the way forward for improvements in Sect. 5.

Figure 1 (a) Mean kinetic energy from AVISO, with mean eddy kinetic energy contours, showing the circulation in the EAC system and the model domain. The cyan lines show the sections through 278 and 34° S. (b) Location of traditional observations used in the TRAD assimilation systems (SSH, SST, and SSS are not shown). (c) Location of additional observations used in the FULL assimilation system and for independent analysis herein. (d) Number of AVISO SSH and NAVOCEANO SST observations and (e) number of Argo, XBT, and SSS observations per 5 d assimilation window.

2 Model and data assimilation system configuration

2.1 The Regional Ocean Modeling System configuration

We use the Regional Ocean Modeling System (ROMS) to simulate the eddying ocean circulation off the southeastern coast of Australia between January 2012 and December 2013. This modelling suite is named the South East Australian Coastal Ocean Forecast System (SEA-COFS, Roughan and Kerry, 2023). ROMS is a widely used free-surface, hydrostatic, terrain-following, primitive equation ocean model and is described by Haidvogel et al. (2000), Marchesiello and Middleton (2000), and Shchepetkin and McWilliams (2005). The model configuration used in this study has been used in various past studies of the EAC and is described in detail in Kerry et al. (2016, 2020a) and Roughan and Kerry (2023).

The study domain covers SE Australia from 25.25 to 41.55° S and approximately 1000 km offshore (Fig. 1a). The domain covers the latitudinal extent of the EAC system from where the current jet is most coherent, the EAC separation region, the region of high eddy activity associated with the EAC eastern extension, and the EAC southern extension. The grid is rotated 20° clockwise such that the domain y axis is oriented roughly parallel with the coastline. The cross-shore horizontal resolution varies from 2.5 km over the continental shelf and gradually increases to 6 km offshore. The horizontal resolution is 5 km in the along-shore direction. Higher resolution over the shelf allows the steep topography to be maintained while minimising pressure gradient errors that emerge in terrain-following coordinate schemes, which otherwise may result in artificial along-slope flow for steep topography (Haney, 1991; Mellor et al., 1994). As such, less topographic smoothing is required to ensure low horizontal pressure gradient errors while still representing the shelf and seamount structures in the model. The model utilises 30 vertical s layers with higher resolution in the upper 500 m to resolve mesoscale dynamics and higher resolution near the seabed for improved representation of the bottom boundary layer. To better resolve surface currents, a near-constant-depth surface layer is provided by applying the vertical stretching scheme of De Souza et al. (2015).

Initial conditions and boundary forcing are derived from the Bluelink ReANalysis version 3 (BRAN3; Oke et al., 2013). The boundary forcing is applied daily, and misfits in baroclinic energy to the BRAN3 condition are absorbed at the boundary via a flow-relaxation scheme. The model is forced at the surface with realistic atmospheric forcing derived from the 12 km resolution Bureau of Meteorology (BOM) Australian Community Climate and Earth-System Simulation (ACCESS) analysis (Puri et al., 2013). The atmospheric forcing fields are applied every 6 h and used to compute the surface wind stress and surface net heat and freshwater fluxes using the bulk flux parameterisation of Fairall et al. (1996).

The free-running configuration, while unable to reproduce the temporal evolution of the mesoscale eddies, has been shown to accurately represent the mean dynamical features of the EAC and both the surface and subsurface (0–2000 m) variability (Kerry and Roughan, 2020a). Specifically, they show that the model accurately represents the mesoscale eddy-related variability in sea surface height (SSH), the frequency in occurrence of EAC separation latitude, the seasonal cycle in sea surface temperature (SST), the ocean's subsurface structure based on data from Argo profiling floats, EAC transport, and the temperature depth structure across the EAC. Thus, using data assimilation, we aim to constrain the model to reproduce the temporal evolution of the mesoscale eddies and examine the forecast skill achieved.

2.2 Observations

The same set of observations are assimilated into the ROMS model configuration using the two DA systems (EnOI and 4D-Var) for comparison in this study. These include satellite-derived SSH, SST, sea surface salinity (SSS), vertical profiles of temperature and salinity from profiling Argo floats, and vertical profiles of temperature from expendable bathythermographs (XBTs) (refer to Fig. 1b). The number of processed observations assimilated for each 5 d assimilation window is shown in Fig. 1d and e. These observations are referred to as the “traditionally” available observations (TRAD) (Siripatana et al., 2020). We describe the observations used and the observation uncertainties specified below. For a detailed description of the observations, the processing performed prior to assimilation, and the prescribed observation uncertainties, the reader is referred to Kerry et al. (2016).

2.2.1 Satellite-derived sea surface height

Archiving, Validation and Interpretation of Satellite Oceanographic Data (AVISO), France, produces global, daily, gridded ($\mathrm{1}/\mathrm{4}$° × $\mathrm{1}/\mathrm{4}$°) mean sea level anomaly (SLA) data by merging of all available along-track satellite altimetry data, computed with respect to a 7-year mean. We add the AVISO SLA data to the dynamic SSH mean from a long free run such that the sea level data are consistent with the ROMS model configuration. The AVISO delayed-time global SLA product error for the region is estimated at 2 cm (AVISO, 2015). We prescribe an additional 4 cm of uncertainty to account for imbalances between this statistical field and a dynamically balanced SSH field required by the model, as well as the smaller spatial-scale processes resolved by the model compared to the gridded product. As such, we prescribe an observation uncertainty of 6 cm. As the AVISO gridded product poorly resolves continental shelf processes, we exclude SSH observations over water depths less than 1000 m.

We use the gridded AVISO product to constrain SSH, rather than the along-track altimetry, for this comparison study. Current work including the development of a high-resolution coastal ocean forecast system (Roughan and Kerry, 2023) is now making use of along-track SSH data successfully with 4D-Var.

2.2.2 Satellite-derived sea surface temperature

SST data from the US Naval Oceanographic Office's Global Area Coverage Advanced Very High Resolution Radiometer level-2 product (NAVOCEANO's GAC AVHRR L2P SST) are used for this study. Data are available 2–3 times per day. We remove day-time SST observations and any night-time observations when wind speed <2 m s⁻¹ (Donlon et al., 2002). The percentage of SST observations removed per 5 d cycle is 0.33 %–54.3 % (mean of 20.77 %). As the resolution of the data is similar to the resolution of the model, the observation uncertainty for the assimilation is chosen to be equal to the specified product error (Andreu-Burillo et al., 2010), which is 0.4–0.5 °C.

2.2.3 Satellite-derived sea surface salinity

We use the Level-3 gridded sea surface salinity (SSS) product derived from the National Aeronautics and Space Administration (NASA) Aquarius satellite (http://www.aquarius.umaine.edu/, last access: 6 March 2024). This product provides daily fields at a 1° resolution. We set the observation uncertainty to 0.4. The specified Aquarius SSS product error is ∼ 0.2, and 0.4 is chosen to account for representation errors. The value is considerably higher than the uncertainties specified for other in situ salinity observations, so SSS provides little constraint to the system (Kerry et al., 2016, 2018).

2.2.4 Argo floats

Argo (free-drifting profiling) floats measure temperature and salinity of the upper 2000 m of the global ocean (http://www.argo.ucsd.edu, last access: 6 March 2024, Fig. 1b). The Argo data points are averaged to the model grid (in the horizontal and vertical) and a 5 min time step. Uncertainty profiles are defined to specify the nominal minimal uncertainties for subsurface temperature and salinity (method described in Kerry et al., 2016). The profiles provide greater uncertainties in the depth ranges of greatest variability where representation errors are likely to be the largest. The observation error variance is specified as the maximum of this nominal minimum error variance and the variance of the observations from the same model cell.

2.2.5 Expendable bathythermographs

Expendable bathythermographs (XBTs) collect temperature profiles along repeat lines sampled by merchant ships; the Sydney–Wellington (PX34) and the Brisbane–Fiji (PX30) routes intersect our model domain (Fig. 1b). Four PX30 lines and seven PX34 lines took place over the assimilation period (2012–2013; Fig. 1e). XBT casts are performed at 10 km intervals along the sections, and the XBT data points are averaged to the model grid and a 5 min time step. The nominal minimal uncertainty variance profiles used for the Argo temperature observations are doubled for the XBT observations, and the observation error variance is specified as the maximum of the nominal minimum error variance and the variance of the observations from the same model cell.

2.2.6 Independent observations used for system assessment

A suite of additional observations were also available over the simulation period (2012–2013) that were collected as part of Australia's Integrated Marine Observing System (IMOS). These include surface velocity measurements from high-frequency coastal radar (HF radar); temperature, salinity, and velocity observations from continental-shelf moorings off the coast of New South Wales (NSW) and South East Queensland (SEQ); temperature, salinity, and velocity observations from five deep-water moorings across the core of the EAC at 28° S (EAC array); and temperature and salinity observations from ocean gliders (refer to Fig. 1c). These products provide independent observations against which we assess the performance of the two systems. Furthermore, these observations were assimilated into the ROMS model (along with the TRAD observations) using 4D-Var (Kerry et al., 2016, 2018; Siripatana et al., 2020). Given the full suite of available observations were assimilated, this system is referred to as the FULL system and considered the “best estimate” of the ocean state over the 2012–2013 period. As such, the FULL system is also used in this paper as a benchmark against which to compare the performance of the two systems presented in this study (4D-Var and EnOI systems that assimilate TRAD observations).

2.3 Data assimilation experiments

In this paper, we refer to three different configurations of the SEA-COFS model which differ in DA type and/or the observations assimilated. Each case is performed over the 2-year period from January 2012 and December 2013 and is described below.

1. 4D-Var TRAD refers to the 4D-Var system that assimilates “traditionally” available observations (SSH, SST, SSS, Argo, and XBT). This system is similar to the system described in Kerry et al. (2016) expect that it only assimilates the TRAD observations.

2. EnOI TRAD refers to the system that assimilates the same observations as the 4D-Var TRAD but using the EnOI DA method described in Sect. 2.4.1 below.

3. 4D-Var FULL refers to the 4D-Var system that assimilates all available observations (SSH, SST, SSS, Argo, XBT, HF radar, shelf and deep moorings, and glider data). It is similar to the system described in detail in Kerry et al. (2016, 2020b, 2018).

   A detailed comparison of the 4D-Var TRAD and the FULL systems was presented in Siripatana et al. (2020). The purpose of this paper is to compare the 4D-Var TRAD and the EnOI TRAD systems, in order to provide a comparison of the two DA schemes using a common suite of traditionally available observations. We introduce the 4D-Var FULL system as a benchmark when comparing against observations that are independent to the TRAD experiments in Sect. 3.2.

2.4 Data assimilation methods

The classic state estimation problem can be given by

$$\begin{array}{}\text{(1)}& {\mathbf{X}}_{\mathrm{a}}={\mathbf{X}}_{\mathrm{f}}+\mathbf{K}(\mathit{y}-H({\mathbf{X}}_{\mathrm{f}}\left)\right),\end{array}$$

where X is the state estimate; subscripts f and a refer to forecast and analysis, respectively; K is the Kalman gain; y is the observation vector; H is the observation operator that samples the background circulation to observation points in space and time. The y−H(X_f) term is referred to as the innovation vector and describes the difference between the observations and the forecast model mapped to observation space. The difference in DA techniques lies in the formulation of K, which determines how the forecast innovations are mapped into model space to produce the new state estimate (X_a). For the standard analysis equation that is solved by the Kalman filter and the dual form of 4D-Var, K can be expressed as

$$\begin{array}{}\text{(2)}& \mathbf{K}={\mathbf{BG}}^{\mathrm{T}}({\mathbf{GBG}}^{\mathrm{T}}+\mathbf{R}{)}^{-\mathrm{1}},\end{array}$$

where B is the background covariance, R is the observation error covariance, and G performs the mapping from model space to observation space.

For time-dependent methods (4D-Var and EnKF), observations are assimilated over a time window respecting the dynamics of the model. The observation operator H samples the nonlinear forecast model X_f at the observation locations in space and time over an assimilation cycle time interval. In 4D-Var, the background error covariance matrix B is typically assumed to be unchanging in time, so there is no explicit flow dependence of the B. Flow dependence is implicit via the terms BG^T and GBG^T in Eq. (2), since G is the operator that maps the tangent linear model solution to the observation points and G^T is the adjoint ocean model forced at observation points (Moore et al., 2011c, 2020). In EnKF, the background error covariance matrix and its evolution in time is estimated from an ensemble of nonlinear model solutions (Houtekamer and Zhang, 2016). For 3D-Var and EnOI, observations are all centred at a single time and, rather than using the model physics to constrain the model versus observation error, time-invariant covariances are prescribed.

2.4.1 EnOI

Ensemble methods (which include the time-dependent EnKF and the time-independent EnOI) use an ensemble of model anomalies to estimate the background error covariances. The EnKF allows for the time-varying statistics by using a fixed number of nonlinear model members (ensembles) to provide a statistical representation of K. The ensembles are generated for every assimilation period so as to capture the state-dependent “errors of the day”. For EnOI, the ensemble of model anomalies is generated from a long non-assimilating model run. This makes the assumption that the background error covariances are not state-dependent and are well represented by a stationary or seasonally varying ensemble. This method is considerably less expensive than the time-dependent EnKF or 4D-Var methods as, once the stationary ensemble is generated, EnOI requires only a single integration of the nonlinear model to generate a background state and only a single solution of the analysis equations to update the background. In contrast, to generate an analysis field using EnKF, the forward nonlinear model must be integrated m times (where m is the number of ensemble members) to represent the time-varying background error covariances and a background state (often based on the ensemble mean). All ensemble members are then updated, requiring m solutions of the analysis equations. Therefore EnOI is m times less expensive than EnKF.

A challenge of ensemble methods is to determine the sufficient number of ensemble members to capture the entirety of the state space, and techniques such as localisation and inflation are used to ensure unrealistic covariances are not applied (Houtekamer and Zhang, 2016). Specifically, localisation is used for three reasons: it reduces the fictitious large covariances at large distance due to sampling error; it improves the rank of the matrix inversion; and, with the use of a parametric form to taper to zero over the localisation distance, the inversions become perfectly parallel, improving computational efficiency (Gaspari and Cohn, 1999). Inflation is only applied to EnKF, not EnOI, with inflation of 5 % being typical. The localisation and inflation techniques however remove some dynamical consistency from the solution. Recent work by the Australian BOM uses a hybrid ensemble transform Kalman filter (Sakov and Oke, 2008) based on 48 dynamic and 96 stationary ensemble members (Brassington et al., 2023). With EnOI, there is less constraint on the number of ensemble members, as the ensembles are only performed once to generate the stationary or seasonally varying ensemble.

For EnOI, there is no time dependence in K (Eq. 2). The mapping from model space to observation space performed by G is time-independent; all observations are co-located at a single time, and the analysis equation (Eq. 1) is considered only at that time. The background error covariance matrix is estimated from a static ensemble of model state anomalies and is given by

$$\begin{array}{}\text{(3)}& \mathbf{B}={\displaystyle \frac{\mathrm{1}}{m-\mathrm{1}}}{\mathbf{AA}}^{\mathrm{T}},\end{array}$$

where A is the matrix of background ensemble anomalies, and m is the ensemble size.

In the EnOI system used in this study, we use a stationary ensemble to represent the intraseasonal model anomalies. Each member is calculated as a difference between a 2-week model average and a 2 d average, centred at the same time. This is repeated every 30 d to ensure the anomalies are independent, generating 266 ensemble members. The DA system is run with a 1 d cycle and centred observation window, so an analysis is generated every day. For SSH, temperature, and salinity, the observation time is assumed to coincide with the analysis time, and innovations are calculated as the difference between observation and model state at the analysis time. The localisation method applied is based on local analysis (Ott et al., 2004); that is, an analysis of a local region is produced with a local background error covariance matrix that has lower dimension than the full state vector. The local analyses are then used to construct complete model states for advancement to the next forecast time. Performing the data assimilation analysis locally is convenient for parallelising the solver. In addition to this, a polynomial taper function is applied to bring the covariance to exactly zero on a specified length scale (Gaspari and Cohn, 1999). The localisation radius is set to 250 km for SSH, temperature, and salinity observations and to 100 km for SST observations. The observation errors are set equal to those described in Sect. 2.2 (identical for both EnOI and 4D-Var systems), except for SST for which the error variance is increased by a factor of 2 for the EnOI system to prevent overfitting to SST. The observation impact was moderated with an adaptive quality control procedure via the so-called K factor (Sandery and Sakov, 2017) with the value of K=2.

For comparison with the 4D-Var system we perform 5 d forecasts based on the EnOI analyses every 4 d. Initial conditions for each subsequent 5 d forecast are taken from the EnOI analysis. In this paper we focus on the forecast skill between the 4D-Var and EnOI systems (not the analysis skill).

2.4.2 4D-Var

The 4D-Var system uses variational calculus to solve for increments in model initial conditions, boundary conditions, and forcing such that the differences between the observations and the new model trajectory are minimised – in a least-squares sense – over a specific assimilation window. The goal is for the model to represent all of the observations in time and space using the physics of the model and accounting for the uncertainties in the observations and background model state, producing a description of the ocean state that is dynamically balanced and a complete solution of the nonlinear model equations.

This is achieved by minimising an objective cost function, J, that measures normalised deviations of the modelled ocean state (given the increment adjustments to model initial conditions, boundary conditions, and forcing) from the observations as well as from the modelled background state (the model prior). The cost function is a function of the increment vector

$$\begin{array}{}\text{(4)}& \mathit{\delta }\mathit{z}=\left(\mathit{\delta }\mathbf{X}\right(t_{\mathrm{0}}{)}^{\mathrm{T}},\mathit{\delta }{\mathit{f}}^{\mathrm{T}}\left(t_{\mathrm{1}}\right),\mathrm{\dots },\mathit{\delta }{\mathit{f}}^{\mathrm{T}}\left(t_n\right),\mathit{\delta }{\mathit{b}}^{\mathrm{T}}\left(t_{\mathrm{1}}\right),\mathrm{\dots },\mathit{\delta }{\mathit{b}}^{\mathrm{T}}\left(t_n\right){)}^{\mathrm{T}}\end{array}$$

representing the increments to the initial conditions (time t₀) and the surface forcing and boundary conditions for model times t₁ to t_n. The cost function can then be written as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}J\left(\mathit{\delta }\mathit{z}\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{0}}^n(\mathbf{G}\mathit{\delta }\mathit{z}-{\mathit{d}}_i{)}^{\mathrm{T}}{\mathbf{R}}_i^{-\mathrm{1}}(\mathbf{G}\mathit{\delta }\mathit{z}-{\mathit{d}}_i)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathit{\delta }\mathit{z}{)}^{\mathrm{T}}{\mathbf{B}}^{-\mathrm{1}}\left(\mathit{\delta }\mathit{z}\right)\\ & =J_o+J_{\mathrm{b}}\end{array},\end{array}$$

where $\mathbf{G}=H_i\mathbf{M}(t_i,t_{\mathrm{0}})$, and M(t_i,t₀) represents the tangent linear version of the nonlinear model equations ℳ, integrated from t₀ to t_i. The difference between the modelled background state and the observations is represented by the innovation vector, introduced above, given at each time t_i by ${\mathit{d}}_i={\mathit{y}}_i-H_i\left({\mathbf{X}}_{\mathrm{f}}\right(t_i\left)\right)$, where y are the observations and H_i is the operator that samples the background circulation to observation points in space and time. As such, the Gδz−d_i term represents the difference between the model and the observations given the increment adjustment integrated through the tangent linear model. R is the observation error covariance matrix, and B is the background error covariance matrix.

We seek to minimise the cost function by equating the gradient to zero. The gradient of the cost function is given by

$$\begin{array}{}\text{(6)}& {\mathrm{\nabla }}_{\mathit{\delta }z}J=\sum _{i=\mathrm{0}}^n{\mathbf{G}}^{\mathrm{T}}{\mathbf{R}}_i^{-\mathrm{1}}(\mathbf{G}\mathit{\delta }\mathit{z}-{\mathit{d}}_i)+{\mathbf{B}}^{-\mathrm{1}}\left(\mathit{\delta }\mathit{z}\right),\end{array}$$

where G^T encompasses the adjoint of the tangent linear model equations. The desired analysis increment, δz_a, that minimises Eq. (5) corresponds to the solution of equation ∇_δzJ=0 and is given by

$$\begin{array}{}\text{(7)}& \mathit{\delta }{\mathit{z}}_{\mathrm{a}}={\mathbf{BG}}^{\mathrm{T}}({\mathbf{GBG}}^{\mathrm{T}}+\mathbf{R}{)}^{-\mathrm{1}}\mathit{d}\end{array}$$

for the dual form (in observation space).

In practice, with 4D-Var, subsequent integrations of the adjoint and tangent linear models are performed to solve for an increment vector that minimises (or acceptably reduces) J. This is performed in the inner loops. After the last inner loop, the final increment is applied to the initial conditions and boundary and surface forcing, and the new integration of the nonlinear model is performed. The integration of the nonlinear model given the increment adjustments that were solved for in the inner loops is referred to as the outer loop. The analysis field is given by the final integration of the nonlinear model (the final outer loop), which provides a model state estimate that is constrained to satisfy the nonlinear model equations (strong constraint) and better represent the observations over the assimilation window. The analysis provides an improved estimate of the initial conditions for the next assimilation window. In this study we find that 15 inner loops and a single outer loop give an acceptable reduction in J (rather than a true minimum).

To solve for the nonlinear ocean solution that better represents the observations, we must take into account the uncertainties in the system. As such, the background (prior model) error covariance matrix, B, and the observation error covariance matrix, R, are important scaling factors in the cost function, J (Eq. 5). The background error covariance matrix should represent the expected uncertainties in the model initial conditions and surface and boundary forcings. We estimate B by factorisation, as described in Weaver and Courtier (2001), such that

$$\begin{array}{}\text{(8)}& \mathbf{B}={\mathbf{K}}_{\mathrm{b}}\mathrm{\Sigma }\mathrm{\Lambda }{L}_v^{\mathrm{1}/\mathrm{2}}{L}_h{L}_v^{\mathrm{1}/\mathrm{2}}\mathrm{\Lambda }\mathrm{\Sigma }{\mathbf{K}}_{\mathrm{b}}^{\mathrm{T}},\end{array}$$

where K_b are the covariance operators of the balanced dynamics, Σ and Λ are the diagonal matrices of the background error standard deviations and normalisation factors respectively, and L_v and L_h are the univariate correlations in the vertical and horizontal directions. We prescribe K_b=I such that the dynamics are coupled through the use of the tangent linear and adjoint models but not in the statistics of B. The correlation matrices, L_v and L_h, and the normalisation factors, Λ, are computed as solutions to diffusion equations following Weaver and Courtier (2001). The characteristic length scales chosen for L_v and L_h are assumed to be homogeneous and isotropic (Table 1), and their choice is justified in Kerry et al. (2016). The specification of the observation error covariances is described in Sect. 2.2 above and in more detail in Kerry et al. (2016).

Table 1 Correlation lengths assumed for the control vector elements: 4D-Var system.

Because we use the linearised model equations, the assimilation window length is limited by the time over which the tangent linear assumption remains reasonable (although longer windows have been shown to produce useful results). For the 4D-Var system presented in this study, we find that a 5 d assimilation window is reasonable. We adjust the model initial conditions, boundary conditions, and surface forcing such that the new model solution (the analysis) better represents the observations over the assimilation interval. Open boundary conditions are adjusted every 12 h and surface forcing every 3 h. A 5 d analysis is generated every 4 d (that is, there is a 1 d overlap between the analyses). Initial conditions for the subsequent 5 d forecast are taken from day 4 of the previous analysis. The ROMS 4D-Var formulation and implementation is well described by Moore et al. (2011d, a, b), and it has been used successfully in many applications (e.g. Di Lorenzo et al., 2007; Powell et al., 2008; Powell and Moore, 2008; Broquet et al., 2009; Matthews et al., 2012; Zavala-Garay et al., 2012; Janeković et al., 2013; Souza et al., 2014; Kerry et al., 2016; Gwyther et al., 2022; Wilkin et al., 2022). This work adopts the same 4D-Var configuration as described in detail in Kerry et al. (2016).

2.4.3 System comparison

As discussed above, the way by which the observations and the model background are combined to generate the analysis is quite different for the 4D-Var and EnOI methods. Another significant difference is the computational expense. For the 15 inner loops and single outer loop used in this study, the 4D-Var data assimilation process is approximately 50 times more expensive than a single free run, making it 25 times more expensive than the EnOI system (once the stationary ensemble has been generated).

This is comparable to the expense of an EnKF using 50 ensembles. The advantage of EnKF (over 4D-Var) is that the tangent linear and adjoint models are not required, all calculations are performed in nonlinear space, and the ensemble members can be run simultaneously if sufficient computing resources are available. The drawback is underdispersion of the ensemble and the loss of dynamic consistency introduced through localisation and inflation. With a 4D-Var system, the use of the adjoint model can provide useful insight into the sensitivity of the ocean state to prior changes in state variables or forcings (e.g. Powell et al., 2013; Kerry et al., 2022) and the direct quantification of observation impacts (e.g. Powell, 2017; Kerry et al., 2018). Observation impacts can also be computed from ensemble methods (Liu and Kalnay, 2008).

Future work aims to compare the EnKF and 4D-Var methods and explore hybrid ensemble–4D-Var methods that capitalise on the advantages of both (i.e. the dynamical interpolation properties of the adjoint used in 4D-Var and the explicit flow-dependent error covariances of the EnKF (Lorenc et al., 2015; Lorenc and Jardak, 2018)). This paper sets a baseline for future work by first comparing the existing and commonly used EnOI method with the 4D-Var method, across a common modelling framework and observational network.

3 System performance: assessing predictive skill

3.1 Assimilated observations

We begin by assessing the performance of the EnOI and 4D-Var systems relative to the observations that the systems assimilate. The 5 d model forecast is compared to the observations that become available over those 5 d (that is, they have not yet been assimilated) to quantitatively assess the performance of the model forecasts over time. Comparing forecasts against observations provides objective assessment of the system performance.

Table 2 presents the mean innovation (mean absolute difference, MAD), innovation bias (mean difference, MD), and number of observations for the 2-year period. Both systems have an identical number of observations. Compared to the EnOI, the 4D-Var improves the SST forecast error from 0.42 to 0.36 °C, the SSH forecast error from 10.3 to 8.3 cm, in situ temperature from 0.90 to 0.71 °C, in situ salinity from 0.079 to 0.056 PSU, and SSS from 0.214 to 0.183 PSU. Overall, the improvement of the MAD for the 4D-Var over the EnOI is 9 %–21 %. The percentage differences in forecast error between the two systems are less for the surface observations (SLA, SST, and SSS) compared to the in situ observations, indicating that the advantages of 4D-Var extend through the water column. In WBC regions, the parent model displayed MADs between reanalysed and observed SST values on day 1 of each assimilation of 0.2–0.6 °C and MADs of 6–12 cm for SSH (Chamberlain et al., 2021b).

Table 2 Summary of performance of the EnOI and 4D-Var systems. Obs num refers to the average number of observations per 5 d assimilation window.

The performance of the two systems relative to SSH, SST, and Argo observations is presented in more detail using the root-mean-square difference (RMSD) between the model forecasts at the observation locations and the observation values. Figure 2a and b show the RMSD between the forecasts (4D-Var and EnOI, respectively) and observations for SSH across the model domain, averaged over the 2-year period. The EnOI forecasts display higher SSH errors across the model domain, with both systems showing higher errors in the eddy-dominated region compared to the rest of the domain. Figure 2c shows that the spatially averaged RMSD between the forecast and the observations is consistently higher for the EnOI forecasts over the 2-year period.

As each forecast is initialised from the previous analysis, forecast errors typically increase over the forecast horizon. SSH forecast errors are averaged across the model domain (Fig. 2d) and for the eddy-dominated region (Fig. 2e) for each day of the 5 d forecast horizon. With SSH, the forecast errors are consistently lower for the 4D-Var system due to lower errors in the initial conditions, while the rate of error increase is similar between the 4D-Var and EnOI systems. At day 5, the domain-averaged (eddy-dominated region averaged) root-mean-squared (rms) SSH forecast errors are 61 % (64 %) higher for the EnOI system compared to the 4D-Var system.

Figure 2 (a) RMSD between forecast and observed SSH for all 186 forecast cycles over the 2-year assimilation period for the TRAD 4D-Var system. (b) Same as (a) for the EnOI system. (c) Spatially averaged RMSD between forecast and observed SSH for each 5 d forecast window. (d) Spatially averaged RMSD between forecast and observed SSH for each day of the 5 d forecast window, averaged over the 186 forecast cycles. (e) Same as (d) but for the high eddy kinetic energy (EKE) region (shown in a and b).

In a similar manner to the SSH forecast errors in Fig. 2, the forecast errors relative to SST observations are presented in Fig. 3. Both systems display higher errors in the core of the EAC upstream of the typical separation region and in the eddy-dominated region. The EnOI forecasts display higher SST errors across the model domain, with the most pronounced difference in the eddy-dominated region (Fig. 3a, b). The time series of RMSD for EnOI and 4D-Var (Fig. 3c) are highly correlated as the statistics are sensitive to the number of observations and the coverage in the high variability area. While the EnOI analyses provide a slightly improved fit to SST (Fig. 3d, e at day 0), SST forecast errors grow more quickly than in the 4D-Var system and the 4D-Var system outperforms the EnOI system for SST forecasts after 1 d. At day 5, the domain-averaged (eddy-dominated region averaged) rms SST forecast errors are 21 % (29 %) higher for the EnOI system compared to the 4D-Var system.

Figure 3 Same as Fig. 2 but for SST observations.

To assess the subsurface predictive skill, we extract the 5 d model forecast values at the observation times and locations for all Argo floats that were observed in the region over the forecast window. Binning these observations with depth, we present profiles for temperature and salinity of the mean (Fig. 4a, e), bias (Fig. 4b, f), and the RMSD between the forecasts and the observations for all observations that fall on the first day of the forecasts (Fig. 4c, g) and all observations that fall on day 5 of the forecasts (Fig. 4d, h). The magnitude of the RMSDs can be compared to the root-mean-squared (rms) observation anomaly, which describes the variability of the observations within each depth bin. For in situ temperature, both the 4D-Var and EnOI forecasts display similar skill on the first day of the forecasts (Fig. 4c); however by day 5 the 4D-Var forecasts display lower errors compared to the EnOI forecasts over the upper 600 m, with a maximum difference in RMSD (bias-corrected RMSD) of 0.56 °C (0.34 °C) at 200 m (Fig. 4d). For salinity, forecast errors at day 5 are of similar magnitude throughout the water column for the two systems (Fig. 4h). Both systems have rms errors considerably less that the rms observation anomaly. Salinity bias dominates the RMSD deeper than 600 m, so bias-corrected RMSD values are less that the total RMSD (Fig. 4g, h).

Figure 4 (a) Mean temperature observed by Argo floats and mean modelled temperature extracted at all Argo observation locations and times for 4D-Var and EnOI systems. (b) Temperature bias. (c) RMSD between forecast and observed temperature at all Argo observation locations and times that fall on forecast day 1, averaged over the 186 forecast cycles. (d) Same as (c) but for forecast day 5. (e)–(h) Same as (a)–(d) but for salinity observed by Argo floats.

3.2 Independent observations

As described in Sect. 2.2, a number of observations were withheld from the 4D-Var and EnOI DA systems presented in this paper, allowing the system performances to be assessed against independent observations. In this section, forecasts from the 4D-Var and EnOI systems (which assimilate the traditional suite of observations, TRAD) are compared to the analyses and forecasts produced by assimilating the full suite of observations (FULL). Comparisons are made between the observations and the model solutions extracted at the observation times and locations, and predictive skill is assessed for days 1 to 5 of the forecast horizons (and analysis windows in the case of the FULL analysis).

Under the HF radar footprint at 30° S, surface radial velocity observations from two sources are combined to compute surface velocities to about 100 km offshore, covering the shelf and shelf slope circulation. This coverage typically includes the EAC as a coherent jet and the intermittent formation of cyclonic frontal eddies inshore of the EAC (Archer et al., 2017; Schaeffer et al., 2017; Kerry et al., 2020a). The complex correlations between the observed and model velocities are presented in Fig. 5. At forecast day 5, the 4D-Var TRAD displays similar predictive skill to the FULL forecasts. The EnOI forecasts are worse than the 4D-Var TRAD across the 5 d, showing that the 4D-Var system provides better representation of the circulation under the HF radar footprint in the analyses and forecasts.

Figure 5 Complex correlation of daily averaged surface velocities measured by the HF radar with FULL analysis (row a), FULL forecast (row b), TRAD forecast (row c), and EnOI forecast (row d), separated by window day (columns). Black lines show 0.9 complex correlation contour, and grey lines show the 70, 200, 1000, and 2000 m isobaths. Only grid cells with a minimum of 15 velocity values over the 2-year period are shown; the values inside the 50 m isobath are removed as the computed velocities are unreliable here due to geometric dilution of precision.

Glider data over the study period (2012–2013) were predominantly available over the NSW continental shelf in water depths <200 m; however, from May–July 2012, several glider missions extended offshore into eddies and sampled down to below 1000 m. These glider observations were shown to be particularly impactful in constraining transport and EKE estimates in the FULL simulation (Kerry et al., 2018). These observations represent independent data for the 4D-Var and EnOI TRAD systems, and Fig. 6 shows how the simulations represent temperature and salinity as measured by the gliders.

Errors are lowest near the surface compared to over the thermocline region due to the assimilation of SST and SSS data in all three systems (4D-Var TRAD, EnOI TRAD, and 4D-Var FULL). The 4D-Var TRAD has rms forecast errors for temperature of a similar magnitude and depth structure as the rms observation anomalies, and the errors do not considerably change from day 1 to day 5 of the forecast window. The EnOI errors are of similar magnitude to the 4D-Var near the surface (∼ 1 °C), but they are 20 % greater between 100–200 m for day 1 and 40 % greater for that depth range at day 5 (Fig. 6c, d). Temperature bias plays a considerable part in the EnOI RMSD values below 100 m, but the bias-corrected RMSD for EnOI still exceeds the bias-corrected RMSD for 4D-Var TRAD at both day 1 and day 5 (Fig. 6c, d).

For salinity, the 4D-Var and EnOI display similar forecast errors in the upper 200 m. This depth range corresponds to where the many shelf glider observations exist. Below 200 m (the off-shelf missions into the Tasman Sea), forecast errors peak at 300 m reaching 0.30 for EnOI at day 5, compared to 0.23 for 4D-Var. Similar to the Argo-observed salinity (Fig. 4f, g, h), salinity bias dominates the errors associated with glider observed salinity below 500 m for 4D-Var TRAD and below 200 m for EnOI.

Figure 6 (a) Mean temperature observed by gliders and mean modelled temperature extracted at all glider observation locations and times for 4D-Var and EnOI systems. (b) Temperature bias. (c) RMSD between forecast and observed temperature at all glider observation locations and times that fall on forecast day 1, averaged over the 186 forecast cycles. (d) Same as (c) but for forecast day 5. (e)–(h) Same as (a)–(d) but for salinity observed by gliders.

Subsurface velocities are measured by acoustic Doppler current profilers mounted on moorings in the EAC array, the SEQ shelf and slope, and on the NSW shelf (Fig. 1c). In Fig. 7 we present the complex correlation between the modelled and observed velocities for selected moorings extending from 28° S to 34° S. The mooring locations are shown in Fig. 1c, with EAC2 and SEQ400 being in 1500 and 400 m water depth at 28° S, CH100 being in 100 m water depth at 30° S and SYD100 being in 100 m water depth at 34° S. At EAC2 and SEQ400, the 4D-Var TRAD displays similar predictive skill to the FULL after 5 d and considerably outperforms the EnOI system throughout the water column. This indicates the benefit of 4D-Var including the northern boundary conditions in the cost function. On the shelf at 30° S (CH100) and 34° S (SYD100), the EnOI and 4D-Var systems show similar predictive skill.

As shown in both Figs. 5 and 7, the 4D-Var FULL complex correlations display a rapid reduction in correlation by day 3–5 of the forecast. As discussed in Siripatana et al. (2020), while the analysis fits the velocity observations along the continental shelf, the forecast model is unable to resolve the complexities of the shelf circulation such as the cyclonic vorticity inshore of the EAC. As such, the forecast skill of the TRAD system is similar to that of the FULL system for 5 d forecast horizons.

Figure 7 Complex correlations between observed and modelled velocities for the 4D-Var TRAD forecast, the EnOI TRAD forecast, the FULL analysis, and the FULL forecast, at selected mooring locations, separated by window days 1, 3, and 5 (columns). Each row represents a single mooring site: EAC2 (row a), SEQ400 (row b), CH100 (row c), and SYD140 (row d).

We have shown that the 4D-Var TRAD system outperforms the EnOI TRAD system at the surface and subsurface when compared against both assimilated and independent observations. Improvements to temperature forecasts with 4D-Var are more pronounced in the subsurface (the upper ∼ 400 m) compared to at the surface (Figs. 4 and 6). We now examine the model forecasts to elucidate the differences between the representation of the ocean state (in model space, rather than observation space) across the two DA systems.

4 Comparisons in model space

4.1 Initial condition increments

The model forecast, X_f, is adjusted by the assimilation of observations (as per Eq. 1) to produce an analysis, X_a. This model state estimate should provide a better representation of the observations and provides updated (improved) initial conditions for the subsequent model forecast. In the 4D-Var system used in this study we perform a 5 d forecast and a 5 d analysis every 4 d, such that the initial conditions for the subsequent forecast are taken from day 4 of the previous analysis. For the EnOI system, an analysis is generated every day. For consistent comparison across the two systems, we take the analysis every 4 d as initial conditions and perform a 5 d forecast. In both cases there are discontinuities in the ocean state between day 4 of the previous forecast and the beginning of the subsequent forecast (which correspond to concurrent times). This is illustrated in Fig. 8i, which shows a time series of temperature at the surface at 34° S. Assimilated (SST) and independent (SYD140 mooring near-surface temperature data) are shown for reference. The discontinuities between the forecasts are less pronounced for the 4D-Var system compared to the EnOI system. Over the entire 2-year test period, the RMSD between the initial conditions (from the analysis) and the previous forecast field at that time illustrates greater discontinuities for the EnOI system compared to the 4D-Var system for SSH, SST, and subsurface temperature (Fig. 8a–h).

The discontinuities presented here do not exactly correspond to the analysis increments. We have presented the differences in the ocean state between day 4 of the previous (5 d) forecast and the beginning of the subsequent forecast (which correspond to concurrent times). For 4D-Var, the ocean state at the beginning of the forecast is taken from the previous cycle analysis, and so the difference presented here represents the difference between the forecast (or the background) at day 4 and analysis at day 4 (once data assimilation has been performed on that assimilation cycle). This is essentially the “analysis increment at day 4”. However for a 4D-Var system the analysis increments typically refer to the adjustments to the initial conditions, boundary forcing, and surface forcing that are made to generate the analysis. For EnOI, the analysis increments refer to the difference between the background model and the analysis (both centred on a single time and computed daily in this case). However, here we take the analyses every 4 d and perform 5 d forecasts, and the differences presented here refer to the difference between day 4 of the forecast and the analysis that provides initial conditions for the subsequent forecast.

With 4D-Var we are able to represent the entirety of the observations collected over a time window (in this case 5 d), placing them in dynamical context using the (linearised) model equations. In contrast, EnOI performs discrete minimisations with observations centred on a single time (in this case every day). The estimate of the ocean over the observation window that is created with the 4D-Var assimilation system results is smaller discontinuities between forecast cycles, on average, compared to the EnOI system, as a continuous field evolves by the nonlinear primitive equations as opposed to starting a forecast from a discrete estimate, which can “shock” the system. Our results of the improved predictability achieved by the 4D-Var system support the understanding that a continual and dynamically balanced analysis field is advantageous to the quality of future predictions.

Figure 8 Root-mean-squared difference between the initial conditions (from the analysis) and the previous forecast field at that time for (a, b) SSH, (c, d) SST, (e, f) temperature at 400 m, and (g, h) temperature at 1000 m, for 4D-Var system (top row) and EnOI system (bottom row). (i) Time series over an example period to illustrate the differences between the end of the forecast window and the analysis conditions in EnOI compared to 4D-Var, for surface temperature at 34° S (location shown in Fig. 11). SST observations within two grid cells and temperature observations from SYD140 in the upper 25 m are also shown for comparison.

4.2 Energetics

The modelled velocities are used to compute eddy kinetic energy (EKE) and mean kinetic energy (MKE) over the 2012–2013 simulation period. MKE is given by $\mathrm{MKE}=\frac{\mathrm{1}}{\mathrm{2}}({\stackrel{\mathrm{‾}}{U}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{V}}^{\mathrm{2}})$, where $\stackrel{\mathrm{‾}}{U}$ and $\stackrel{\mathrm{‾}}{V}$ are the time mean velocity components, and the EKE is given by $\mathrm{EKE}=\frac{\mathrm{1}}{\mathrm{2}}({U^{\prime }}^{\mathrm{2}}+{V^{\prime }}^{\mathrm{2}})$, where U^′ and V^′ are the velocity anomalies. The MKE describes the energy associated with the mean currents, and the EKE describes the energy associated with the perturbations from the mean. Figure 9 shows the MKE and EKE averaged over the upper 400 m and from 400–1200 m.

Comparisons of MKE above 400 m show that the EAC core is narrower and more confined to the slope in the 4D-Var system, while MKE for the EnOI system is more spread out and with higher MKE directly over the continental shelf (Fig. 9a, e, i). This difference is despite the identical SSH observations being assimilated, noting that SSH observations in water depth <1000 m are not assimilated, and the identical forward numerical model. In the 4D-Var simulation, the MKE is greater below 400 m than the EnOI simulation downstream of 27.5° S to the typical EAC separation zone (Fig. 9b, f, j). This is consistent with Kerry and Roughan (2020a), who use a long-term integration of the free-running simulation to describe a downstream deepening of the EAC before separation.

The spatial structure of the EKE is similar across the two systems. Above 400 m, the EnOI system has elevated EKE over the EAC jet (Fig. 9k, blue regions), while the 4D-Var system has elevated EKE in the eddy-dominated regions (Fig. 9k, red regions). The elevated EKE for the EnOI system (in the more coherent region) relates to the greater discontinuities between the subsequent forecasts, which manifests itself as greater low-frequency >1 d variability over the 5 d forecasts as the 5 d model run adjusts to the “shocks” to the system. In contrast, the elevated EKE in the 4D-Var system outside of the coherent jet relates to the greater near-inertial variability. This is explored in Sect. 4.3 and Fig. 12. At depth (400–1200 m), EKE is elevated for EnOI compared to 4D-Var in the EAC southern extension.

Figure 9 The 4D-Var simulation (a) 0–400 m, (b) 400–1200 m MKE, and (c) 0–400 m and (d) 400–1200 m time-averaged EKE are shown. The EnOI (e) 0–400 m, (f) 400–1200 m MKE, and (g) 0–400 m and (h) 400–1200 m time-averaged EKE are shown. In (i)–(l), the difference in each respective field between the 4D-Var and EnOI simulations is shown, where a positive difference indicates more energy in the 4D-Var simulation.

Eddies can form through barotropic instability in the mean flow or baroclinic instability in the vertical density structure. It is important for a model to correctly represent these instabilities, as they represent the pathways by which eddies are generated. Following Kang and Curchitser (2015), we calculate the barotropic conversion rate (KmKe) as

$$\begin{array}{}\text{(9)}& \mathrm{KmKe}={\mathit{\rho }}_{\mathrm{0}}\left[\stackrel{\mathrm{‾}}{U^{\prime }{U}^{\prime }}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{U}}{\partial x}}+\stackrel{\mathrm{‾}}{U^{\prime }{V}^{\prime }}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{U}}{\partial y}}+\stackrel{\mathrm{‾}}{V^{\prime }{U}^{\prime }}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{V}}{\partial x}}+\stackrel{\mathrm{‾}}{V^{\prime }{V}^{\prime }}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{V}}{\partial y}}\right],\end{array}$$

where ρ₀=1025 kg m⁻³. The baroclinic conversion rate (PeKe), from eddy potential energy to EKE, is calculated as

$$\begin{array}{}\text{(10)}& \mathrm{PeKe}=-g\stackrel{\mathrm{‾}}{{\mathit{\rho }}^{\prime }{W}^{\prime }},\end{array}$$

where the acceleration due to gravity is g=9.81 m s⁻², and ρ^′ and W^′ are the density and vertical velocity anomalies. KmKe and PeKe have been previously used to explore eddy generation rates in the EAC (e.g. Li et al., 2021, 2022; Gwyther et al., 2023).

Barotropic and baroclinic energy conversions are computed from the model forecast fields and averaged over the 2-year period (Fig. 10). Both the 4D-Var and EnOI systems show similar magnitude and overall spatial structure of the barotropic and baroclinic energy conversions, as well as similar partitioning between barotropic and baroclinic instabilities. The similarities are likely due to the common model and atmospheric forcing. The barotropic conversion (compare Fig. 10a, c) represents instabilities in the depth-mean flow, which 4D-Var and EnOI represent similarly. The baroclinic conversion (compare Fig. 10b–d) is also similar between the DA configurations in overall spatial structure and the zonally integrated magnitudes (Fig. 10e), although the EnOI baroclinic conversion rate contains more high-wavenumber spatial patterns, which likely relate to unbalanced adjustments upon assimilation. This is further explored in Fig. 14 and the associated discussion.

Figure 10 For the 4D-Var simulation, the (a) barotropic (KmKe) and (b) baroclinic (PeKe) conversions are shown. The EnOI simulation (c) barotropic and (d) baroclinic conversion rates are shown. Conversion rates are calculated as the depth-mean conversion for each model column from the surface to 450 m. In (e), the zonally averaged conversions are shown for both simulations. Averaging is performed in the across-shelf direction in a band extending approximately from the coast to ∼3° offshore, as indicated by dashed lines in panel (a).

4.3 Temporal and spatial scales of variability

When observations are assimilated the goal is to provide an improved fit to the observations while retaining a dynamically consistent ocean state that can be used as initial conditions for the subsequent forecast. The background numerical model produces an estimate of the ocean state whose frequency and wavenumber spectra are limited by the resolution of the model and the processes resolved. If the observations sample time and space scales that cannot be resolved by the model, it is standard DA practice to either remove these scales of variability from the observations or account for them in the observation error terms (e.g. Kerry and Powell, 2022). If the model background is deficient at some space scale and/or timescale (which it is able to resolve), then these may be corrected by DA so that the analyses and forecasts are better. However, if the assimilation process introduces energy at different, non-physical scales, this may negatively impact the forecast skill. By presenting the temporal and spatial scales of variability of the forecast ocean state, we can understand how the assimilation has changed the ocean's energy distribution and understand the differences in error growth across the two DA systems.

The subsurface structure of the model fields and their variability is shown in Fig. 11. The EnOI system has more temperature variability near the surface (upper 200–500 m) compared to 4D-Var. The greater near-surface temperature variability in EnOI compared to 4D-Var is greater in the eddy-dominated region (34° S), where adjustments are greater (Fig. 8d, f) compared to more coherent, upstream region (28° S). For velocity variability, 4D-Var shows elevated variability almost everywhere except in the upper 250 m near the shelf at 34° S. Both data-assimilating configurations show elevated variability in temperature and velocity over the upper ∼ 1000 m compared to the non-data assimilating simulation (hereafter referred to as the Free-run). The differences are further illustrated in Figs. 12 and 13, where the frequencies of the variability are revealed.

Figure 11 Column 1: mean temperature across all 5 d forecasts for sections 28° S (top panels) and 34° S (bottom panels) for the 4D-Var system, the EnOI system, the difference in mean temperature between the systems (4D-Var – EnOI), the difference in mean temperature between 4D-Var and the Free-run, and the difference in mean temperature between EnOI and the Free-run. Column 2: temperature variability for the 4D-Var system, the EnOI system, the difference in variability between the systems (4D-Var – EnOI), the difference in variability between 4D-Var and the Free-run, and the difference in variability between EnOI and the Free-run. Temperature variance is computed for every 5 d forecast, averaged over all forecast windows, and the square root taken. Column 3 shows the same as column 1 but for alongshore velocity. Column 4 shows the same as column 2 but for alongshore velocity. The 4D-Var – EnOI variability panels show points chosen to present frequency spectra (Figs. 12 and 13). Sections are shown in Fig. 1a.

Frequency spectral analysis is first performed for all 5 d forecast windows and then averaged (Fig. 12). A 5 d window with the model output 4-hourly gives a frequency range from $\mathrm{1}/\mathrm{5}$ d to $\mathrm{1}/\mathrm{8}$ h, with 15 points in frequency space due to the short time series (31 points). Surface velocity (but not temperature) and subsurface temperature and velocity display elevated energy in the 16–24 h band for the 4D-Var system compared to EnOI and the Free-run (Fig. 12), corresponding to the near-inertial band. This inertial energy is introduced through the assimilation adjustments which, due to the nature of 4D-Var, must satisfy the model equations. Increased near-inertial variability upon 4D-Var data assimilation was also shown in Matthews et al. (2012) and Kerry and Powell (2022). Matthews et al. (2012) found that the increased inertial energy had minimal impact on the mesoscale circulation. Using observing system simulation experiments, Kerry and Powell (2022) showed that, while the 4D-Var system displayed elevated near-inertial variability (compared to their free running truth simulation), near-inertial frequencies did not influence energy at other frequencies, and predictability at both higher frequencies (in their case internal tides) and lower frequencies (associated with the mesoscale circulation) was good.

The differences between EnOI and the Free-run and EnOI and 4D-Var (as revealed in Fig. 11) are difficult to decipher from Fig. 12 as they exist at low frequencies (periods greater than 1 d). In order to resolve the low frequencies, we concatenate the forecast cycles in time to produce a full 2-year time series. As a longer time series allows a higher resolution in frequency space, Fig. 13 show a higher frequency resolution compared to Fig. 12, for which the spectra are computed for all 5 d periods and averaged. Concatenation of the time series requires removal of the 1 d overlap (the last day of each cycle is excluded) such that time is monotonously increasing. Because of the assimilation updates, discontinuities exist between the cycles every 4 d (this is not the case for the Free-run). These discontinuities (displayed in Fig. 8 and discussed in Sect. 4.1) manifest as harmonics of the $\mathrm{1}/\mathrm{4}$ d frequency and are most pronounced for EnOI in temperature at 34° S. For example, the surface temperature spectra for the EnOI system at 34° S show spikes at harmonics of 0.25 (0.5, 0.75 ,1.0 ,1.25, 1.5, etc. cycles per day). Nevertheless, the spectra are useful in showing the differences in variability across the DA systems and the Free-run particularly for low frequencies. The Free-run displays less energy than both DA systems in both temperature and velocity in the eddy-dominated region, consistent with the reduced variability shown in Fig. 11. This relates to less variability at low frequencies (periods greater than 1 d) in the Free-run compared to the DA systems, and for the 4D-Var system, less variability at inertial frequencies also.

The elevated energy in the EnOI system compared to 4D-Var and the Free-run relates to periods greater than 1 d for both temperature and velocity (Fig. 13). Greater variability at low frequencies in the EnOI system compared to 4D-Var exists in both temperature and velocity and is most pronounced in the upper 500 m and in the eddy-dominated region (34° S compared to the more coherent region at 28° S, Fig. 13). This increased low-frequency variability in EnOI compared to 4D-Var dominates the total variability (displayed in Fig. 11) for near-surface temperature, but it is masked by greater inertial-period variability in the 4D-Var system for velocity. That is, despite the low-frequency velocity variability being greater in EnOI (Fig. 13), the total velocity variability is greater for 4D-Var (Fig. 11). We find that the greater low-frequency variability for EnOI compared to 4D-Var is associated with greater discontinuities between the subsequent forecasts. The discontinuities also exist in 4D-Var, but they are less pronounced (Fig. 8).

Figure 12 Frequency spectra in model space for temperature and alongshore velocity at the surface, 400 m, and 1000 m at 28 and 34° S. Spectra are computed for each 5 d forecast window and then averaged. Points are chosen in the core of the EAC based on the long-term alongshore velocity mean (from Kerry and Roughan, 2020a) where the shelf slope depth is 1500 m at 28° S and just offshore of the shelf slope where the water depth is 3500 m at 34° S. The daily period is shown by the vertical dashed lines.

Figure 13 Frequency spectra for the same variables and points shown in Fig. 12. However rather than averaging the spectra for all 5 d periods, the forecast cycles are concatenated to make a full 2-year time series (with the 1 d overlap removed). The discontinuities between the cycles every 4 d manifest as harmonics of the $\mathrm{1}/\mathrm{4}$ d frequency and are most pronounced for EnOI in temperature at 34° S. The daily and 12-hourly periods are shown by the vertical dashed lines. In computing the spectra, four ensembles and four bands are used to increase the statistical significance.

The spatial scales of the forecast ocean state can be represented by wavenumber spectra. Here we present cross-shore wavenumber kinetic energy spectra through sections at 28 and 34° S (Fig. 1a) for days 1 and 5 of the forecasts, for the Free-run and for AVISO gridded geostrophic velocities (Fig. 14). The observational data product used is the AVISO gridded velocities from altimetry and drifters using multiscale interpolation, version 0100 (Ballarotta et al., 2022), with a $\mathrm{1}/\mathrm{10}$° spatial resolution and temporal coverage from 1 July 2016–30 June 2020. Note that the spatial resolution is the same as that of the assimilated SSH observations (Sect. 2.2.1). For the model forecasts, wavenumber kinetic energy spectra are computed for days 1 and 5 of all (186) cycles, and the averages are plotted. For the AVISO observations, wavenumber kinetic energy spectra are computed for every day of the available time period and the average plotted. Model spectra are shown at the surface and at depths of 400 and 1000 m; AVISO data provide geostrophic velocities, and the corresponding spectra are plotted on the surface velocity panels.

At the surface all systems, except the EnOI at day 1, display consistent kinetic energy spectra at 28° S. The AVISO velocities show less energy at spatial scales between 15–80 km compared to the Free-run, the 4D-Var system across all forecast days, and the EnOI system at day 5. At 34° S, where eddy variability is high, the Free-run underrepresents the kinetic energy across all spatial scales at all depths. At the surface, the 4D-Var system across all forecast days and the EnOI system at day 5 represent the AVISO spectrum well, with the AVISO velocities again showing slightly lower energy at spatial scales between 15–80 km.

For the first day of the EnOI forecasts (representative of the analyses), there is elevated kinetic energy at finer length scales and this energy dissipates by day 5 of the forecast. This elevated energy is most pronounced at the surface and near-surface (upper 200 m, not shown). Specifically, elevated kinetic energy exists in the EnOI initial states at length scales less than 100 km at 28° S and between 20–80 km at 34° S. For the 4D-Var system the wavenumber kinetic energy spectra remain relatively unchanged over the forecast window, with the day 1 and day 5 wavenumber spectra tracking closely. Compared to the Free-run, both the 4D-Var and EnOI assimilation systems introduce more kinetic energy across all spatial scales throughout the water column in the eddy-dominated region (illustrated by the sections through 34° S in Fig. 14).

We include the idealised spectral slopes of $k^{-\mathrm{5}/\mathrm{3}}$ and k⁻³ in Fig. 14 for reference. The wavenumber kinetic energy spectra approximately match the $k^{-\mathrm{5}/\mathrm{3}}$ slope for the mesoscale range, the k⁻³ slope for the submesoscale range for the 4D-Var ocean state on day 1 and day 5, and the EnOI forecasts on day 5. However, we note that the submesoscale range is only partially resolved by the 2.5–6 km resolution model and even less so by the AVISO observations. The $k^{-\mathrm{5}/\mathrm{3}}$ and k⁻³ slopes have been shown to represent surface quasi-geostrophic and quasi-geostrophic dynamics, respectively; however realistic simulations show that other slopes are possible (Xu and Fu, 2011).

Figure 14 Cross-shore wavenumber kinetic energy spectra for the models at the surface, 400 and 1000 m and for AVISO geostrophic velocities at the surface, and at 28 and 34° S. The length scales 200, 100, and 20 km are shown by the vertical dashed lines. The $-\mathrm{5}/\mathrm{3}$ and −3 spectral slopes are shown on the first panel for comparison. In computing the spectra, two ensembles and two bands are used to increase the statistical significance.

We have shown that energy is elevated for shorter (less than 100 km) length scales in the EnOI analyses, and upon integration of the forecast model this energy dissipates to match the energy associated with the 4D-Var system. Wavenumber kinetic energy analysis of the atmosphere by Skamarock (2004) showed the contrary: the increase of energy at small scales upon integration of the forecast model. They showed that the initial states of high-resolution NWP model forecasts lacked the fine-scale (mesoscale in the case of the atmosphere) energy because “observations to initialise the fine scales are not generally available and data assimilation methods that can use high-resolution observations are not yet mature”. The fine-scale portion of the kinetic energy spectrum was spun up in the forecasts in 6–12 h, providing increased value to the NWP forecasts. In our study we observe the introduction of energy at small spatial scales upon DA with EnOI, and this elevated small-scale energy is lost by day 5. This implies that the small-scale energy dissipates over the 5 d forecast. The elevated kinetic energy at scales less than 100 km is not a physical space scale that is resolved by the observations, as shown by the AVISO kinetic energy spectra, and does not exist in the 4D-Var system (Fig. 14). Rather it comes about due to EnOI's adjustments upon assimilation. It is likely that the increased error growth (hence poorer forecast skill) for the EnOI system (compared to 4D-Var) in this study relates to these adjustments. The consistency of the wavenumber spectra over the 4D-Var 5 d forecast windows likely relates to the constraint that the analysis is a complete solution of the model nonlinear equations, requiring dynamically balanced adjustments.

5 Conclusions

This study shows in a quantified manner that the smoother and more dynamically balanced fit between the observations and the model's time-evolving flow achieved by the 4D-Var system results in improved predictability against both assimilated and non-assimilated observations. The EnOI system does not produce as tight as fit to the SSH data as the 4D-Var system (although this may be related to tuneable parameters in the DA formulation); however, the SSH error grows at the same rate in the EnOI and 4D-Var forecasts (Fig. 2). The surface expression of the EAC and its associated eddies is associated with the barotropic mode, and our results show that the barotropic energy conversion rates are generally consistent across the two systems (Fig. 10a, c). However, the baroclinic conversion rate has small spatial scale variability in the EnOI forecasts compared to the 4D-Var (Fig. 10b, d), and the EnOI analyses (the forecast initial conditions) display elevated energy at fine (<100 km) spatial scales (Fig. 14). This is accompanied by reduced predictive skill for both surface and in situ temperature, in situ salinity, and surface velocities (Figs. 3, 4, 5, 6, 7). For SST (Fig. 3) and temperature in the upper 600 m (Fig. 4c, d), the analyses have errors of similar magnitude for the EnOI and 4D-Var systems, but error growth is considerably greater in the EnOI forecasts. Note that the upper 600 m is the region of greatest variability in both temperature and salinity (Fig. 4c, d, g, h, blue lines). The improved forecasts of SST and in situ temperature in the upper 600 m for 4D-Var after 5 d (Figs. 3, 4d) are a demonstration of improved dynamical balance of the model initial conditions. This is evident by the smaller magnitude of the increments for 4D-Var (Fig. 8a, c, e, g) compared to EnOI (Fig. 8b, d, f, h). The bias-corrected salinity errors also show similar errors at forecast day 1 for both systems, with greater error growth in the EnOI system compared to 4D-Var by day 5 (Fig. 4g, h).

Independent surface velocity observations as measured by the high-frequency radar array at 30° S are less well represented by the EnOI system compared to the 4D-Var system from day 1 through to day 5 of the forecasts (Fig. 5). Independent in situ temperature observations from gliders show only slightly lower analysis errors for 4D-Var compared to EnOI, and the subsurface temperature forecasts degrade faster over the 5 d window for EnOI compared to 4D-Var (Fig. 6), consistent with the forecast errors associated with assimilated in situ temperature observations (Fig. 4). For salinity, EnOI and 4D-Var perform equally well on the shelf (observations above 200 m in Fig. 6g, h are dominated by shelf gliders), but EnOI displays higher errors below 200 m by day 5. The 4D-Var system displays improved velocity forecasts compared to the EnOI system for the upstream moorings (EAC2 and SEQ400, Fig. 7), while downstream and on the shelf the forecasts are comparable. This indicates the benefit of 4D-Var including the northern boundary conditions in the cost function. Generally, we show that the benefits of 4D-Var over EnOI are most pronounced in the (5 d) forecasts, rather than the fit of the analyses to the observations, consistent with the paper “Why does 4D‐Var beat 3D‐Var?” by Lorenc and Rawlins (2005).

The EnOI system displays greater discontinuities between the end of the forecast and the subsequent analysis, particularly for near-surface temperature (about the thermocline), and the discontinuities have greater magnitude in the downstream eddy-dominated region (Fig. 8). These assimilation “shocks” manifest as increased low-frequency variability (periods greater than 1 d, Figs. 11 and 13). The 4D-Var system displays elevated energy in the near-inertial frequency band for both temperature and velocity (Figs. 12 and 13). Consistent with Kerry and Powell (2022) and Matthews et al. (2012), the energy at near-inertial frequencies does not appear to affect the mean low-frequency energetics associated with the mesoscale circulation. While the EnOI DA system introduces elevated energy at fine (<100 km) spatial scales, 4D-Var maintains the kinetic energy distribution in wavenumber space upon assimilation (Fig. 14).

This study chose to compare two DA methods across a common modelling framework and observational network. The two methods were chosen as EnOI has been widely used by the Australian ocean forecasting community (Oke et al., 2008b, 2010; Chamberlain et al., 2021b), and 4D-Var has been implemented to study predictability and observation impact in the EAC (Kerry et al., 2016, 2018; Siripatana et al., 2020; Gwyther et al., 2022, 2023). It made sense for the two user groups (operational and research) to come together to objectively compare the two methods. Each system was tuned by its developers (Australian Bureau of Meteorology for EnOI and UNSW for 4D-Var). We note that the degree of fit between an analysis and the assimilated observations of a specific DA system is sensitive to the prior choice of various parameters, such as the observation and background error covariances, and that the system performance is influenced by the DA system configuration, such as size of the ensemble for ensemble methods and the assimilation window length for 4D-Var (Moore et al., 2020; Santana et al., 2023). For example, the EnOI system presented here could be further tuned to provide an improved fit to SSH observations (Fig. 2), and different ensemble sizes could be tested. For the 4D-Var system, different window lengths could be tested and the sensitivity to changes in B could be studied. However, the goal of this study was not to compare various versions of each DA method. Rather we compare a single version of the two methods, carefully tuned by each user group, and set a baseline for future comparisons. The focus of this paper is not the fit in the analyses but the rate of forecast error growth and the response of the ocean state to the assimilation methodology. As such, the study's utility and relevance is significant without a large number of comparisons with different prior specified parameters or DA system configurations.

The EnOI system is ∼ 25 times cheaper than the 4D-Var system presented here. It is noted that EnOI has been effective for long-term reanalysis products where analyses were created every day (Oke et al., 2008b; Chamberlain et al., 2021b) and forecasts were not required. With increasing computational capacity and the pursuit of more accurate ocean forecasts, this study's comparison motivates the use of 4D-Var over EnOI for ocean forecasts of the EAC region. This result is likely to be applicable over similar, highly variable, oceanic regions such as WBCs. More generally, the comparison advocates for the use of advanced time-dependent DA schemes over time-independent methods. We illustrate how a DA scheme can influence forecast skill which motivates future development of DA methods. It is noted that Australia's operational ocean model (OceanMAPS) recently transitioned to an EnKF DA method (from EnOI). The new system achieves lower mean error and error variance in WBC extensions regions (Chamberlain et al., 2021a; Brassington et al., 2023), with lower increments to SSH and subsurface velocities, and less kinetic energy at depth in the analyses, due to more dynamically balanced adjustments, compared to the EnOI system.

Our future work specifically aims to directly address the need to improve predictive skill in WBC regions. Time-independent schemes (e.g. 3D-Var and EnOI) are useful for intermittent cycling DA at synoptic scales and are capable of resolving slowly evolving flows governed by simple balance relationships. Time-dependent DA methods (e.g. 4D-Var and EnKF) are greatly beneficial for highly intermittent flows with irregularly sampled observations as the time-variable dynamics of the model are used to evolve the error covariances. Furthermore, these methods allow the entirety of observations over a time interval to be minimised rather than discrete minimisations. The time-evolving state is required to truly exploit many novel observation types that are nonlinearly or indirectly related to the model state. Indeed, the two techniques that are the most promising in NWP and ocean DA are 4D-Var and EnKF (Moore et al., 2019). In recent years it has been recognised that a marriage of 4D-Var and EnKF perhaps represents a more optimal approach since it capitalises on the advantages of both approaches (i.e. the dynamical interpolation properties of the adjoint and the explicit flow-dependent error covariances that capture the “errors of the day”). The relative performance of 4D-Var and EnKF methods in regional ocean models has been assessed by Moore et al. (2020), and the differences are due primarily to the properties of the background error covariances, so it is anticipated that the performance of a system using a hybrid covariance will be superior to either 4D-Var or the EnKF alone. Such ensemble-variational methods have been studied extensively for atmospheric DA (e.g. Lorenc et al., 2015) with improvements in forecast skill achieved particularly in dynamically active systems (Raynaud et al., 2011; Lorenc and Jardak, 2018).