the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Comparison of 4dimensional variational and ensemble optimal interpolation data assimilation systems using a Regional Ocean Modeling System (v3.4) configuration of the eddydominated East Australian Current system
Colette Gabrielle Kerry
Moninya Roughan
Shane Keating
David Gwyther
Gary Brassington
Adil Siripatana
Joao Marcos A. C. Souza
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 shortterm 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 2year period. We compare the timedependent 4dimensional variational (4DVar) data assimilation system with the more computationally efficient, timeindependent ensemble optimal interpolation (EnOI) system, across a common modelling and observational framework. Both systems assimilate the same observations: satellitederived 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 4DVar 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, rootmeansquared forecast errors are 20 %–60 % higher for the EnOI system compared to the 4DVar system. The 4DVar 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 lowfrequency (>1 d) surfaceintensified variability in temperature and elevated kinetic energy at length scales less than 100 km at the beginning of the forecast windows. The 4DVar system displays elevated energy in the nearinertial range throughout the water column, with the wavenumber kinetic energy spectra remaining unchanged upon assimilation. Overall, this comparison shows quantitatively that the 4DVar system results in improved predictability as the analysis provides a smoother and more dynamically balanced fit between the observations and the model's timeevolving flow. This advocates the use of advanced, timedependent data assimilation methods, particularly for highly variable oceanic regions, and motivates future work into further improving data assimilation schemes.
 Article
(6542 KB)  Fulltext XML
 BibTeX
 EndNote

The predictive performances of two ocean data assimilation systems (EnOI and 4DVar) are assessed in a Regional Ocean Modeling System (ROMS) configuration of the East Australian Current over 5 d forecast horizons.

The forecast skill of the 4DVar system surpasses the EnOI system against both assimilated and independent observations at the surface and subsurface.

The EnOI system has greater analysis increments, elevated lowfrequency (>1 d) surfaceintensified variability in temperature, and elevated kinetic energy at length scales less than 100 km at the beginning of the forecast windows.

The dynamically balanced 4DVar system displays elevated energy in the nearinertial range throughout the water column, with the wavenumber kinetic energy spectra remaining unchanged upon assimilation.
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 southwardflowing 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 warmcore 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, timeindependent methods such as 3dimensional variational data assimilation (3DVar) 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 MediumRange Weather Forecasts uses 3DVar 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 flowdependent 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, timedependent DA techniques (Edwards et al., 2015; Moore et al., 2019). Advanced DA methods make use of the timevariable 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, timeindependent 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 4dimensional variational data assimilation (4DVar) and the ensemble Kalman filter (EnKF), and ocean DA is following suit (Moore et al., 2019).
In 4DVar 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 4DVar, 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 4DVar (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 timeindependent DA technique (EnOI) with the timedependent technique (4DVar) 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 shortrange (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.
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 (SEACOFS, Roughan and Kerry, 2023). ROMS is a widely used freesurface, hydrostatic, terrainfollowing, 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 crossshore 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 alongshore direction. Higher resolution over the shelf allows the steep topography to be maintained while minimising pressure gradient errors that emerge in terrainfollowing coordinate schemes, which otherwise may result in artificial alongslope 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 nearconstantdepth 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 flowrelaxation 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 EarthSystem 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 freerunning 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 eddyrelated 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 4DVar) for comparison in this study. These include satellitederived 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 Satellitederived 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 alongtrack satellite altimetry data, computed with respect to a 7year 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 delayedtime 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 spatialscale 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 alongtrack altimetry, for this comparison study. Current work including the development of a highresolution coastal ocean forecast system (Roughan and Kerry, 2023) is now making use of alongtrack SSH data successfully with 4DVar.
2.2.2 Satellitederived sea surface temperature
SST data from the US Naval Oceanographic Office's Global Area Coverage Advanced Very High Resolution Radiometer level2 product (NAVOCEANO's GAC AVHRR L2P SST) are used for this study. Data are available 2–3 times per day. We remove daytime SST observations and any nighttime observations when wind speed <2 m s^{−1} (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 (AndreuBurillo et al., 2010), which is 0.4–0.5 °C.
2.2.3 Satellitederived sea surface salinity
We use the Level3 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 (freedrifting 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 highfrequency coastal radar (HF radar); temperature, salinity, and velocity observations from continentalshelf moorings off the coast of New South Wales (NSW) and South East Queensland (SEQ); temperature, salinity, and velocity observations from five deepwater 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 4DVar (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 (4DVar and EnOI systems that assimilate TRAD observations).
2.3 Data assimilation experiments
In this paper, we refer to three different configurations of the SEACOFS model which differ in DA type and/or the observations assimilated. Each case is performed over the 2year period from January 2012 and December 2013 and is described below.

4DVar TRAD refers to the 4DVar 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.

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

4DVar FULL refers to the 4DVar 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 4DVar TRAD and the FULL systems was presented in Siripatana et al. (2020). The purpose of this paper is to compare the 4DVar 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 4DVar 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
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 4DVar, K can be expressed as
where B is the background covariance, R is the observation error covariance, and G performs the mapping from model space to observation space.
For timedependent methods (4DVar 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 4DVar, 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 3DVar and EnOI, observations are all centred at a single time and, rather than using the model physics to constrain the model versus observation error, timeinvariant covariances are prescribed.
2.4.1 EnOI
Ensemble methods (which include the timedependent EnKF and the timeindependent EnOI) use an ensemble of model anomalies to estimate the background error covariances. The EnKF allows for the timevarying 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 statedependent “errors of the day”. For EnOI, the ensemble of model anomalies is generated from a long nonassimilating model run. This makes the assumption that the background error covariances are not statedependent and are well represented by a stationary or seasonally varying ensemble. This method is considerably less expensive than the timedependent EnKF or 4DVar 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 timevarying 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 timeindependent; all observations are colocated 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
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 2week 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 4DVar 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 socalled K factor (Sandery and Sakov, 2017) with the value of K=2.
For comparison with the 4DVar 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 4DVar and EnOI systems (not the analysis skill).
2.4.2 4DVar
The 4DVar 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 leastsquares 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
representing the increments to the initial conditions (time t_{0}) and the surface forcing and boundary conditions for model times t_{1} to t_{n}. The cost function can then be written as
where $\mathbf{G}={H}_{i}\mathbf{M}({t}_{i},{t}_{\mathrm{0}})$, and M(t_{i},t_{0}) represents the tangent linear version of the nonlinear model equations ℳ, integrated from t_{0} 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
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 ∇_{δz}J=0 and is given by
for the dual form (in observation space).
In practice, with 4DVar, 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
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).
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 4DVar 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 4DVar 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; ZavalaGaray 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 4DVar 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 4DVar and EnOI methods. Another significant difference is the computational expense. For the 15 inner loops and single outer loop used in this study, the 4DVar 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 4DVar) 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 4DVar 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 4DVar methods and explore hybrid ensemble–4DVar methods that capitalise on the advantages of both (i.e. the dynamical interpolation properties of the adjoint used in 4DVar and the explicit flowdependent 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 4DVar method, across a common modelling framework and observational network.
3.1 Assimilated observations
We begin by assessing the performance of the EnOI and 4DVar 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 2year period. Both systems have an identical number of observations. Compared to the EnOI, the 4DVar 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 4DVar 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 4DVar 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).
The performance of the two systems relative to SSH, SST, and Argo observations is presented in more detail using the rootmeansquare difference (RMSD) between the model forecasts at the observation locations and the observation values. Figure 2a and b show the RMSD between the forecasts (4DVar and EnOI, respectively) and observations for SSH across the model domain, averaged over the 2year period. The EnOI forecasts display higher SSH errors across the model domain, with both systems showing higher errors in the eddydominated 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 2year 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 eddydominated region (Fig. 2e) for each day of the 5 d forecast horizon. With SSH, the forecast errors are consistently lower for the 4DVar system due to lower errors in the initial conditions, while the rate of error increase is similar between the 4DVar and EnOI systems. At day 5, the domainaveraged (eddydominated region averaged) rootmeansquared (rms) SSH forecast errors are 61 % (64 %) higher for the EnOI system compared to the 4DVar system.
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 eddydominated region. The EnOI forecasts display higher SST errors across the model domain, with the most pronounced difference in the eddydominated region (Fig. 3a, b). The time series of RMSD for EnOI and 4DVar (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 4DVar system and the 4DVar system outperforms the EnOI system for SST forecasts after 1 d. At day 5, the domainaveraged (eddydominated region averaged) rms SST forecast errors are 21 % (29 %) higher for the EnOI system compared to the 4DVar system.
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 rootmeansquared (rms) observation anomaly, which describes the variability of the observations within each depth bin. For in situ temperature, both the 4DVar and EnOI forecasts display similar skill on the first day of the forecasts (Fig. 4c); however by day 5 the 4DVar forecasts display lower errors compared to the EnOI forecasts over the upper 600 m, with a maximum difference in RMSD (biascorrected 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 biascorrected RMSD values are less that the total RMSD (Fig. 4g, h).
3.2 Independent observations
As described in Sect. 2.2, a number of observations were withheld from the 4DVar and EnOI DA systems presented in this paper, allowing the system performances to be assessed against independent observations. In this section, forecasts from the 4DVar 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 4DVar TRAD displays similar predictive skill to the FULL forecasts. The EnOI forecasts are worse than the 4DVar TRAD across the 5 d, showing that the 4DVar system provides better representation of the circulation under the HF radar footprint in the analyses and forecasts.
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 4DVar 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 (4DVar TRAD, EnOI TRAD, and 4DVar FULL). The 4DVar 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 4DVar 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 biascorrected RMSD for EnOI still exceeds the biascorrected RMSD for 4DVar TRAD at both day 1 and day 5 (Fig. 6c, d).
For salinity, the 4DVar 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 offshelf missions into the Tasman Sea), forecast errors peak at 300 m reaching 0.30 for EnOI at day 5, compared to 0.23 for 4DVar. Similar to the Argoobserved salinity (Fig. 4f, g, h), salinity bias dominates the errors associated with glider observed salinity below 500 m for 4DVar TRAD and below 200 m for EnOI.
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 4DVar 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 4DVar including the northern boundary conditions in the cost function. On the shelf at 30° S (CH100) and 34° S (SYD100), the EnOI and 4DVar systems show similar predictive skill.
As shown in both Figs. 5 and 7, the 4DVar 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.
We have shown that the 4DVar 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 4DVar 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.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 4DVar 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 nearsurface temperature data) are shown for reference. The discontinuities between the forecasts are less pronounced for the 4DVar system compared to the EnOI system. Over the entire 2year 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 4DVar 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 4DVar, 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 4DVar 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 4DVar 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 4DVar 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 4DVar system support the understanding that a continual and dynamically balanced analysis field is advantageous to the quality of future predictions.
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{\u203e}}{U}}^{\mathrm{2}}+{\stackrel{\mathrm{\u203e}}{V}}^{\mathrm{2}})$, where $\stackrel{\mathrm{\u203e}}{U}$ and $\stackrel{\mathrm{\u203e}}{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 4DVar 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 4DVar 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 longterm integration of the freerunning 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 4DVar system has elevated EKE in the eddydominated 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 lowfrequency >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 4DVar system outside of the coherent jet relates to the greater nearinertial variability. This is explored in Sect. 4.3 and Fig. 12. At depth (400–1200 m), EKE is elevated for EnOI compared to 4DVar in the EAC southern extension.
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
where ρ_{0}=1025 kg m^{−3}. The baroclinic conversion rate (PeKe), from eddy potential energy to EKE, is calculated as
where the acceleration due to gravity is g=9.81 m s^{−2}, 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 2year period (Fig. 10). Both the 4DVar 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 depthmean flow, which 4DVar 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 highwavenumber spatial patterns, which likely relate to unbalanced adjustments upon assimilation. This is further explored in Fig. 14 and the associated discussion.
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, nonphysical 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 4DVar. The greater nearsurface temperature variability in EnOI compared to 4DVar is greater in the eddydominated region (34° S), where adjustments are greater (Fig. 8d, f) compared to more coherent, upstream region (28° S). For velocity variability, 4DVar shows elevated variability almost everywhere except in the upper 250 m near the shelf at 34° S. Both dataassimilating configurations show elevated variability in temperature and velocity over the upper ∼ 1000 m compared to the nondata assimilating simulation (hereafter referred to as the Freerun). The differences are further illustrated in Figs. 12 and 13, where the frequencies of the variability are revealed.
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 4hourly 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 4DVar system compared to EnOI and the Freerun (Fig. 12), corresponding to the nearinertial band. This inertial energy is introduced through the assimilation adjustments which, due to the nature of 4DVar, must satisfy the model equations. Increased nearinertial variability upon 4DVar 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 4DVar system displayed elevated nearinertial variability (compared to their free running truth simulation), nearinertial 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 Freerun and EnOI and 4DVar (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 2year 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 Freerun). 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 Freerun particularly for low frequencies. The Freerun displays less energy than both DA systems in both temperature and velocity in the eddydominated 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 Freerun compared to the DA systems, and for the 4DVar system, less variability at inertial frequencies also.
The elevated energy in the EnOI system compared to 4DVar and the Freerun 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 4DVar exists in both temperature and velocity and is most pronounced in the upper 500 m and in the eddydominated region (34° S compared to the more coherent region at 28° S, Fig. 13). This increased lowfrequency variability in EnOI compared to 4DVar dominates the total variability (displayed in Fig. 11) for nearsurface temperature, but it is masked by greater inertialperiod variability in the 4DVar system for velocity. That is, despite the lowfrequency velocity variability being greater in EnOI (Fig. 13), the total velocity variability is greater for 4DVar (Fig. 11). We find that the greater lowfrequency variability for EnOI compared to 4DVar is associated with greater discontinuities between the subsequent forecasts. The discontinuities also exist in 4DVar, but they are less pronounced (Fig. 8).
The spatial scales of the forecast ocean state can be represented by wavenumber spectra. Here we present crossshore wavenumber kinetic energy spectra through sections at 28 and 34° S (Fig. 1a) for days 1 and 5 of the forecasts, for the Freerun 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 Freerun, the 4DVar system across all forecast days, and the EnOI system at day 5. At 34° S, where eddy variability is high, the Freerun underrepresents the kinetic energy across all spatial scales at all depths. At the surface, the 4DVar 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 nearsurface (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 4DVar 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 Freerun, both the 4DVar and EnOI assimilation systems introduce more kinetic energy across all spatial scales throughout the water column in the eddydominated 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^{−3} 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^{−3} slope for the submesoscale range for the 4DVar 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^{−3} slopes have been shown to represent surface quasigeostrophic and quasigeostrophic dynamics, respectively; however realistic simulations show that other slopes are possible (Xu and Fu, 2011).
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 4DVar 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 highresolution NWP model forecasts lacked the finescale (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 highresolution observations are not yet mature”. The finescale 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 smallscale energy is lost by day 5. This implies that the smallscale 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 4DVar 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 4DVar) in this study relates to these adjustments. The consistency of the wavenumber spectra over the 4DVar 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.
This study shows in a quantified manner that the smoother and more dynamically balanced fit between the observations and the model's timeevolving flow achieved by the 4DVar system results in improved predictability against both assimilated and nonassimilated observations. The EnOI system does not produce as tight as fit to the SSH data as the 4DVar 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 4DVar 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 4DVar (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 4DVar 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 4DVar 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 4DVar (Fig. 8a, c, e, g) compared to EnOI (Fig. 8b, d, f, h). The biascorrected salinity errors also show similar errors at forecast day 1 for both systems, with greater error growth in the EnOI system compared to 4DVar by day 5 (Fig. 4g, h).
Independent surface velocity observations as measured by the highfrequency radar array at 30° S are less well represented by the EnOI system compared to the 4DVar 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 4DVar compared to EnOI, and the subsurface temperature forecasts degrade faster over the 5 d window for EnOI compared to 4DVar (Fig. 6), consistent with the forecast errors associated with assimilated in situ temperature observations (Fig. 4). For salinity, EnOI and 4DVar 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 4DVar 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 4DVar including the northern boundary conditions in the cost function. Generally, we show that the benefits of 4DVar 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 nearsurface temperature (about the thermocline), and the discontinuities have greater magnitude in the downstream eddydominated region (Fig. 8). These assimilation “shocks” manifest as increased lowfrequency variability (periods greater than 1 d, Figs. 11 and 13). The 4DVar system displays elevated energy in the nearinertial 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 nearinertial frequencies does not appear to affect the mean lowfrequency energetics associated with the mesoscale circulation. While the EnOI DA system introduces elevated energy at fine (<100 km) spatial scales, 4DVar 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 4DVar 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 4DVar). 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 4DVar (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 4DVar 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 4DVar system presented here. It is noted that EnOI has been effective for longterm 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 4DVar 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 timedependent DA schemes over timeindependent 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. Timeindependent schemes (e.g. 3DVar and EnOI) are useful for intermittent cycling DA at synoptic scales and are capable of resolving slowly evolving flows governed by simple balance relationships. Timedependent DA methods (e.g. 4DVar and EnKF) are greatly beneficial for highly intermittent flows with irregularly sampled observations as the timevariable 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 timeevolving 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 4DVar and EnKF (Moore et al., 2019). In recent years it has been recognised that a marriage of 4DVar 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 flowdependent error covariances that capture the “errors of the day”). The relative performance of 4DVar 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 4DVar or the EnKF alone. Such ensemblevariational 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).
The ROMS model code is available from https://www.myroms.org (last access: 6 March 2024; DOI: https://doi.org/10.5281/zenodo.8294716, Roughan and Kerry, 2023). SEACOFS model configuration is accessible at https://doi.org/10.26190/5e683944e1369 (Kerry and Roughan, 2020b), https://doi.org/10.26190/5ebe1f389dd87 (Kerry et al., 2020b), and https://doi.org/10.5281/zenodo.8294716 (Roughan and Kerry, 2023).
The observations were sourced from the Integrated Marine Observing System (IMOS). IMOS is a national collaborative research infrastructure, supported by the Australian Government (https://www.imos.org.au, last access: 6 March 2024). Observations are available at https://portal.aodn.org.au/ (last access: 6 March 2024).
Argo data were collected and made freely available by the international Argo programme and the national programmes that contribute to it. (http://www.argo.ucsd.edu, last access: 6 March 2024). The Argo programme is part of the Global Ocean Observing System (https://doi.org/10.17882/42182, Notarstefano, 2020).
We acknowledge AVISO for the delayedtime SLA data. The Ssalto/Duacs altimeter products were produced and distributed by the Copernicus Marine and Environment Monitoring Service (CMEMS) (https://marine.copernicus.eu, last access: 6 March 2024).
CGK developed the ROMS model configuration of the EAC system, processed the observations, and developed the 4DVar DA configuration. CGK performed the 5 d forecasts, given the EnOI analyses. CGK analysed the results to produce Figs. 1–8 and 11–14. DG produced Figs. 9–10. CGK wrote the manuscript with some original input from AS. AS generated the results in Table 1. We acknowledge Pavel Sakov, who generated the EnOI analyses, given the ROMS model configuration and the processed observations from CGK. MR, SK, GB, and JMACS provided useful guidance and input into the scope of the project and interpretation of results.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
For this research, David Gwyther and Adil Siripatana were partially supported by the Australian Research Council Industry Linkage grant no. LP170100498 to Moninya Roughan, Colette Gabrielle Kerry, and Shane Keating. Prior model development was supported by the Australian Research Council grant nos. DP140102337 and LP160100162. CSIRO Marine and Atmospheric Research and Wealth from Oceans Flagship Program, Hobart, Tasmania, Australia, provided BRAN2020 output for boundary conditions.
This research has been supported by the Australian Research Council (grant nos. LP170100498, DP140102337 and LP160100162).
This paper was edited by Deepak Subramani and reviewed by two anonymous referees.
AndreuBurillo, I., Brassington, G., Oke, P., and Beggs, H.: Including a new data stream in the BLUElink Ocean Data Assimilation System, Aust. Meteorol. Ocean., 59, 77–86, 2010. a
Archer, M., Roughan, M., Keating, S., and Schaeffer, A.: On the variability of the East Australian Current: Jet structure, meandering, and influence on shelf circulation, J. Geophys. Res.Oceans, 122, 8464–8481, 2017. a
AVISO: SSALTO/DUACS user handbook: (M)SLA and (M)ADT nearreal time and delayed time products. CLSDOSNT06034, SALPMUPEA21065CLS, CNES, 66 pp., http://www.aviso.oceanobs.com/fileadmin/documents/data/tools/hdbk_duacs.pdf (last access: 6 March 2024), 2015. a
Ballarotta, M., Ubelmann, C., Veillard, P., Prandi, P., Etienne, H., Mulet, S., Faugère, Y., Dibarboure, G., Morrow, R., and Picot, N.: Improved global sea surface height and current maps from remote sensing and in situ observations, Earth Syst. Sci. Data, 15, 295–315, https://doi.org/10.5194/essd152952023, 2023. a
Brassington, G. B., Sakov, P., Divakaran, P., Aijaz, S., SweeneyVan Kinderen, J., Huang, X., and Allen, S.: OceanMAPS v4. 0i: a global eddy resolving EnKF ocean forecasting system, in: OCEANS 2023Limerick, IEEE, 1–8, https://doi.org/10.1109/OCEANSLimerick52467.2023.10244383, 2023. a, b, c
Broquet, G., Edwards, C. A., Moore, A., Powell, B. S., Veneziani, M., and Doyle, J. D.: Application of 4DVariational data assimilation to the California Current System, Dynam. Atmos. Oceans, 48, 69–92, 2009. a
Brousseau, P., Berre, L., Bouttier, F., and Desroziers, G.: Flowdependent backgrounderror covariances for a convectivescale data assimilation system, Q. J. Roy. Meteor. Soc., 138, 310–322, 2012. a
Bull, C. Y. S., Kiss, A. E., Jourdain, N. C., England, M. H., and Van Sebille, E.: Wind forced variability in eddy formation, eddy shedding, and the separation of the East Australian Current., J. Geophys. Res.Oceans, 122, 9980–9998, 2017. a
Cetina Heredia, P., Roughan, M., Van Sebille, E., and Coleman, M.: Longterm trends in the East Australian Current separation latitude and eddy driven transport, J. Geophys. Res., 119, 4351–4366, https://doi.org/10.1002/2014JC010071, 2014. a, b
Chamberlain, M., Oke, P., Brassington, G., Sandery, P., Divakaran, P., and Fiedler, R.: Multiscale data assimilation in the Bluelink ocean reanalysis (BRAN), Ocean Model., 166, 101849, https://doi.org/10.1016/j.ocemod.2021.101849, 2021a. a
Chamberlain, M. A., Oke, P. R., Fiedler, R. A. S., Beggs, H. M., Brassington, G. B., and Divakaran, P.: Next generation of Bluelink ocean reanalysis with multiscale data assimilation: BRAN2020, Earth Syst. Sci. Data, 13, 5663–5688, https://doi.org/10.5194/essd1356632021, 2021b. a, b, c
De Souza, J. M. A. C., Powell, B., CastilloTrujillo, A. C., and Flament, P.: The Vorticity Balance of the Ocean Surface in Hawaii from a Regional Reanalysis, J. Phys. Oceanogr., 45, 424–440, 2015. a
Di Lorenzo, E., Moore, A. M., Arango, H. G., Cornuelle, B. D., Miller, A. J., Powell, B. S., Chua, B. S., and Bennett, A. F.: Weak and Strong Constraint Data Assimilation in the inverse Regional Ocean Modelling System (ROMS): development and application for a baroclinic coastal upwelling system, Ocean Model., 16, 160–187, 2007. a
Donlon, C., Minnett, P., Gentemann, C., Nightingale, T., Barton, I., Ward, B., and Murray, M.: Toward improved validation of satellite sea surface skin temperature measurements for climate research, J. Climate, 15, 353–369, 2002. a
Edwards, C. A., Moore, A. M., Hoteit, I., and Cornuelle, B. D.: Regional ocean data assimilation, Annu. Rev. Mar. Sci., 7, 21–42, 2015. a
Evensen, G.: Sequential data assimilation for nonlinear dynamics: the ensemble Kalman filter, in: Ocean Forecasting: Conceptual basis and applications, 97–116, Springer, https://doi.org/10.1007/9783662226483_6, 2002. a
Fairall, C. W., Bradley, E. F., Rogers, D. P., Edson, J. B., and Young, G. S.: Bulk parameterization of airsea fluxes for tropical oceanglobal atmosphere CoupledOcean Atmosphere Response Experiment, J. Geophys. Res., 101, 3747–3764, 1996. a
Feron, R. C. V.: The southern ocean Western Boundary Currents: Comparison of fine resolution Antarctic model results with Geosat altimeter data, J. Geophys. Res., 100, 4959–4975, 1995. a
Gaspari, G. and Cohn, S. E.: Construction of correlation functions in two and three dimensions, Q. J. Roy. Meteor. Soc., 125, 723–757, 1999. a, b
Gwyther, D. E., Kerry, C., Roughan, M., and Keating, S. R.: Observing system simulation experiments reveal that subsurface temperature observations improve estimates of circulation and heat content in a dynamic western boundary current, Geosci. Model Dev., 15, 6541–6565, https://doi.org/10.5194/gmd1565412022, 2022. a, b
Gwyther, D. E., Keating, S. R., Kerry, C., and Roughan, M.: How does 4DVar data assimilation affect the vertical representation of mesoscale eddies? A case study with observing system simulation experiments (OSSEs) using ROMS v3.9, Geosci. Model Dev., 16, 157–178, https://doi.org/10.5194/gmd161572023, 2023. a, b
Haidvogel, D. B., Arango, H. G., Hedstrom, K., Beckmann, A., MalanotteRizzoli, P., and Shchepetkin, A. F.: Model evaluation experiments in the North Atlantic Basin: simulations in nonlinear terrainfollowing coordinates, Dynam. Atmos. Oceans, 32, 239–281, 2000. a
Haney, R. L.: On the Pressure Gradient Force over Steep Topography in Sigma Coordinate Ocean Models, J. Phys. Oceanogr., 21, 610–619, 1991. a
Houtekamer, P. L. and Zhang, F.: Review of the ensemble Kalman filter for atmospheric data assimilation, Mon. Weather Rev., 144, 4489–4532, 2016. a, b
Imawaki, S., Bower, A. S., Beal, L., and Qiu, B.: Chapter 13 – Western Boundary Currents, in: Ocean Circulation and Climate, in: International Geophysics, edited by: Siedler, G., Griffies, S. M., Gould, J., and Church, J. A., vol. 103, 305–338, https://doi.org/10.1016/B9780123918512.000131, 2013. a
Janeković, I., Powell, B. S., Matthews, D., McManus, M. A., and Sevadjian, J.: 4DVar Data Assimilation in a Nested, Coastal Ocean Model: A Hawaiian Case Study, J. Geophys. Res., 118, 5022–5035, https://doi.org/10.1002/jgrc.20389, 2013. a
Kang, D. and Curchitser, E. N.: Energetics of Eddy–Mean Flow Interactions in the Gulf Stream Region, J. Phys. Oceanogr., 45, 1103–1120, https://doi.org/10.1175/jpod140200.1, 2015. a
Kerry, C. and Roughan, M.: Downstream Evolution of the East Australian Current System: Mean Flow, Seasonal, and Intraannual Variability, J. Geophys. Res.Oceans, 125, e2019JC015227, https://doi.org/10.1029/2019JC015227, 2020a. a, b, c, d
Kerry, C. and Roughan, M.: A highresolution, 22year, freerunning, hydrodynamic simulation of the EAC System using the ROMS, unsworks [data set], https://doi.org/10.26190/5e683944e1369, 2020b. a
Kerry, C., Powell, B., Roughan, M., and Oke, P.: Development and evaluation of a highresolution reanalysis of the East Australian Current region using the Regional Ocean Modelling System (ROMS 3.4) and Incremental StrongConstraint 4Dimensional Variational (IS4DVar) data assimilation, Geosci. Model Dev., 9, 3779–3801, https://doi.org/10.5194/gmd937792016, 2016. a, b, c, d, e, f, g, h, i, j, k, l, m
Kerry, C. G., Roughan, M., and Powell, B. S.: Observation Impact in a Regional Reanalysis of the East Australian Current System, J. Geophys. Res.Oceans, 123, 7511–7528, https://doi.org/10.1029/2017JC013685, 2018. a, b, c, d, e, f
Kerry, C., Roughan, M., and Powell, B.: Predicting the submesoscale circulation inshore of the East Australian Current, J. Marine Syst., 204, 103286, https://doi.org/10.1016/j.jmarsys.2019.103286, 2020a. a, b
Kerry, C., Roughan, M., Powell, B., and Oke, P.: A highresolution reanalysis of the East Australian Current System assimilating an unprecedented observational data set using 4DVar data assimilation over a twoyear period (2012–2013), Version 2017, UNSW [data set], https://doi.org/10.26190/5ebe1f389dd87, 2020b. a, b
Kerry, C., Roughan, M., and Azevedo Correia de Souza, J. M.: Drivers of upper ocean heat content extremes around New Zealand revealed by Adjoint Sensitivity Analysis, Frontiers in Climate, 205, https://doi.org/10.3389/fclim.2022.980990, 2022. a
Kerry, C. G. and Powell, B. S.: Including tides improves subtidal prediction in a region of strong surface and internal tides and energetic mesoscale circulation, J. Geophys. Res.Oceans, 127, e2021JC018314, https://doi.org/10.1029/2021JC018314, 2022. a, b, c, d
Li, J. and Roughan, M.: Energetics of Eddy–Mean Flow Interactions in the East Australian Current System, J. Phys. Oceanogr., 53, 595–612, 2023. a
Li, J., Roughan, M., and Kerry, C.: Dynamics of Interannual Eddy Kinetic Energy Modulations in a Western Boundary Current, Geophys. Res. Lett., 48, e2021GL094115, https://doi.org/10.1029/2021gl094115, 2021. a
Li, J., Roughan, M., and Kerry, C.: Variability and Drivers of Ocean Temperature Extremes in a Warming Western Boundary Current, J. Climate, 35, 1097–1111, https://doi.org/10.1175/JCLID210622.1, 2022. a
Liu, J. and Kalnay, E.: Estimating observation impact without adjoint model in an ensemble Kalman filter, Q. J. Roy. Meteor. Soc., 134, 1327–1335, 2008. a
Lorenc, A. C.: Analysis methods for numerical weather prediction, Q. J. Roy. Meteor. Soc., 112, 1177–1194, 1986. a
Lorenc, A. C. and Jardak, M.: A comparison of hybrid variational data assimilation methods for global NWP, Q. J. Roy. Meteor. Soc., 144, 2748–2760, 2018. a, b
Lorenc, A. C. and Rawlins, F.: Why does 4DVar beat 3DVar?, Q. J. Roy. Meteor. Soc., 131, 3247–3257, 2005. a, b
Lorenc, A. C., Bowler, N. E., Clayton, A. M., Pring, S. R., and Fairbairn, D.: Comparison of hybrid4DEnVar and hybrid4DVar data assimilation methods for global NWP, Mon. Weather Rev., 143, 212–229, 2015. a, b
Macdonald, H. S., Roughan, M., Baird, M. E., and Wilkin, J.: A numerical modeling study of the East Australian Current encircling and overwashing a warmcore eddy, J. Geophys. Res., 118, 301–315, https://doi.org/10.1029/2012JC008386, 2013. a
Malan, N., Roughan, M., Hemming, M., and Schaeffer, A.: Mesoscale Circulation Controls Chlorophyll Concentrations in the East Australian Current Separation Zone, J. Geophys. Res.Oceans, 128, e2022JC019361, https://doi.org/10.1029/2022JC019361, 2023. a
Marchesiello, P. and Middleton, J. H.: Modeling the East Australian Current in the Western Tasman Sea, J. Phys. Oceanogr., 30, 2956–2971, 2000. a
Mata, M. M., Tomczak, M., Wijffels, S. E., and Church, J. A.: East Australian Current volume transports at 30° S: Estimates from the World Ocean Circulation Experiment hydrographic sections PR11/P6 and the PCM3 current meter array, J. Geophys. Res., 105, 28509–28526, 2000. a
Mata, M. M., Wijffels, S. E., Church, J. A., and Tomczak, M.: Eddy shedding and energy conversions in the East Australian Current, J. Geophys. Res., 111, C09034, https://doi.org/10.1029/2006JC003592, 2006. a
Matthews, D., Powell, B. S., and Janeković, I.: Analysis of Fourdimensional Variational State Estimation of the Hawaiian Waters, J. Geophys. Res., 117, C03013, https://doi.org/10.1029/2011JC007575, 2012. a, b, c, d
Mellor, G. L., Ezer, T., and Oey, L. Y.: The pressure gradient error conundrum of sigma coordinate ocean models, J. Atmos. Ocean. Tech., 11, 1126–1134, 1994. a
Mogensen, K., Balmaseda, M., and Weaver, A.: The NEMOVAR ocean data assimilation system as implemented in the ECMWF ocean analysis for System 4, European Centre for MediumRange Weather Forecasts, Reading, UK, https://doi.org/10.21957/x5y9yrtm, 2012. a
Moore, A., Martin, M., Akella, S., Arango, H., Balmaseda, M., Bertino, L., Ciavatta, S., Cornuelle, B., Cummings, J., Frolov, S., Lermusiaux, P., Oddo, P., Oke, P., Storto, A., Teruzzi, A., Vidard, A., and Weaver, A.: Synthesis of Ocean Observations Using Data Assimilation for Operational, RealTime and Reanalysis Systems: A More Complete Picture of the State of the Ocean, Frontiers in Marine Science, 6, 90, https://doi.org/10.3389/fmars.2019.00090, 2019. a, b, c
Moore, A., ZavalaGaray, J., Arango, H. G., Edwards, C. A., Anderson, J., and Hoar, T.: Regional and basin scale applications of ensemble adjustment Kalman filter and 4DVar ocean data assimilation systems, Prog. Oceanogr., 189, 102450, https://doi.org/10.1016/j.pocean.2020.102450, 2020. a, b, c
Moore, A. M., Arango, H. G., Di Lorenzo, E., Cornuelle, B. D., Miller, A. J., and Neilson, D. J.: A comprehensive ocean prediction and analysis system based on the tangent linear and adjoint of a regional ocean model, Ocean Model., 7, 227–258, 2004. a
Moore, A. M., Arango, H. G., Broquet, G., Edwards, C., Veneziani, M., Powell, B. S., Foley, D., Doyle, J., Costa, D., and Robinson, P.: The Regional Ocean Modeling System (ROMS) 4dimensional variational data assimilation systems: Part II – Performance and application to the California Current System, Prog. Oceanog., 91, 50–73, https://doi.org/10.1016/j.pocean.2011.05.003, 2011a. a
Moore, A. M., Arango, H. G., Broquet, G., Edwards, C., Veneziani, M., Powell, B. S., Foley, D., Doyle, J., Costa, D., and Robinson, P.: The Regional Ocean Modeling System (ROMS) 4dimensional variational data assimilation systems: Part III – Observation impact and observation sensitivity in the California Current System, Prog. Oceanog., 91, 74–94, https://doi.org/10.1016/j.pocean.2011.05.005, 2011b. a
Moore, A. M., Arango, H. G., Broquet, G., Powell, B. S., Weaver, A. T., and ZavalaGaray, J.: The Regional Ocean Modelling System (ROMS) 4dimensional variational data assimilation systems: Part 1 – System overview and formulation, Prog. Oceanogr., 91, 34–49, 2011c. a
Moore, A. M., Arango, H. G., Broquet, G., Powell, B. S., ZavalaGaray, J., and Weaver, A. T.: The Regional Ocean Modeling System (ROMS) 4dimensional variational data assimilation systems: Part I – System overview and formulation, Prog. Oceanog., 91, 34–49, https://doi.org/10.1016/j.pocean.2011.05.004, 2011d. a
Notarstefano, G.: Argo Float Data and Metadata from Global Data Assembly Centre (Argo GDAC), SEANOE [data set], https://doi.org/10.17882/42182, 2020. a
Oke, P., Sakov, P., Cahill, M. L., Dunn, J. R., Fiedler, R., Griffin, D. A., Mansbridge, J. V., Ridgway, K. R., and Schiller, A.: Towards a dynamically balanced eddyresolving ocean reanalysis: BRAN3, Ocean Model., 67, 52–70, 2013. a
Oke, P. R. and Griffin, D. A.: The coldcore eddy and strong upwelling off the coast of New South Wales in early 2007, DeepSea Res. Pt. II, 58, 574–591, https://doi.org/10.1016/j.dsr2.2010.06.006, 2011. a
Oke, P. R. and Middleton, J. H.: Topographically Induced Upwelling off Eastern Australia, J. Phys. Oceanogr., 30, 512–530, 2000. a
Oke, P. R., Brassington, G. B., Griffin, D. A., and Schiller, A.: The Bluelink ocean data assimilation system (BODAS), Ocean Model., 21, 46–70, 2008a. a
Oke, P. R., Brassington, G. B., Griffin, D. A., and Schiller, A.: The Bluelink ocean data assimilation system (BODAS), Ocean Model., 21, 46–70, 2008b. a, b
Oke, P. R., Brassington, G. B., Griffin, D. A., and Schiller, A.: Ocean data assimilation: a case for ensemble optimal interpolation, Aust. Meteorol. Ocean., 59, 67–76, 2010. a, b
Oke, P. R., Roughan, M., CetinaHeredia, P., Pilo, G. S., Ridgway, K. R., Rykova, T., Archer, M. R., Coleman, R. C., Kerry, C. G., Rocha, C., Schaeffer, A., and Vitarelli, E.: Revisiting the circulation of the East Australian Current: its path, separation, and eddy field, Prog. Oceanogr., 176, 102139, https://doi.org/10.1016/j.pocean.2019.102139, 2019. a
Ott, E., Hunt, B. R., Szunyogh, I., Zimin, A. V., Kostelich, E. J., Corazza, M., Kalnay, E., Patil, D., and Yorke, J. A.: A local ensemble Kalman filter for atmospheric data assimilation, Tellus, 56A, 415–428, 2004. a
Pilo, G. S., Mata, M. M., and Azevedo, J. L. L.: Eddy surface properties and propagation at Southern Hemisphere western boundary current systems, Ocean Sci., 11, 629–641, https://doi.org/10.5194/os116292015, 2015a. a
Pilo, G. S., Oke, P. R., Rykova, T., Coleman, R., and Ridgway, K.: Do East Australian Current anticyclonic eddies leave the Tasman Sea?, J. Geophys. Res.Oceans, 120, 8099–8114, 2015b. a
Pilo, G. S., Oke, P. R., Coleman, R., Rykova, T., and Ridgway, K.: Patterns of vertical velocity induced by eddy distortion in an ocean model, J. Geophys. Res.Oceans, 123, 2274–2292, 2018. a
Powell, B. S.: Quantifying How Observations Inform a Numerical Reanalysis of Hawaii, J. Geophys. Res.Oceans, 122, 8427–8444, https://doi.org/10.1002/2017JC012854, 2017. a
Powell, B. S. and Moore, A. M.: Estimating the 4DVAR analysis error of GODAE products, Ocean Dynam., 59, 121–138, 2008. a
Powell, B. S., Arango, H. G., Moore, A. M., Di Lorenzo, E., Milliff, R. F., and Foley, D.: 4DVAR data assimilation in the IntraAmericas Sea with the Regional Ocean Modeling System (ROMS), Ocean Model., 25, 173–188, 2008. a
Powell, B. S., Kerry, C. G., and Cornuelle, B. D.: Using a numerical model to understand the connection between the ocean and acoustic traveltime measurements, J. Acoust. Soc. Am., 134, 3211–3222, 2013. a
Puri, K., Dietachmayer, G., Steinle, P., Dix, M., Rikus, L., Logan, L., Naughton, M., Tingwell, C., Xiao, Y., Barras, V., Bermous, I., Bowen, R., Deschamps, L., Franklin, C., Fraser, J., Glowacki, T., Harris, B., Lee, J., Le, T., Roff, G., Sulaiman, A., Sims, H., Sun, X., Sun, Z., Zhu, H., Chattopadhyay, M., and Engel, C.: Operational implementation of the ACCESS Numerical Weather Prediction system, Aust. Meteorol. Ocean., 63, 265–284, 2013. a
Raynaud, L., Berre, L., and Desroziers, G.: An extended specification of flowdependent background error variances in the MétéoFrance global 4DVar system, Q. J. Roy. Meteor. Soc., 137, 607–619, 2011. a
Roughan, M. and Kerry, C.: South East Australian Coastal Ocean Forecast System (SEACOFS), Zenodo [data set], https://doi.org/10.5281/zenodo.8294716, 2023. a, b, c, d, e
Roughan, M., Keating, S., Schaeffer, A., Cetina Heredia, P., Rocha, C., Griffin, D., Robertson, R., and Suthers, I.: A tale of two eddies: The biophysical characteristics of two contrasting cyclonic eddies in the e ast a ustralian currents ystem, J. Geophys. Res.Oceans, 122, 2494–2518, 2017. a, b
Sakov, P. and Oke, P. R.: A deterministic formulation of the ensemble Kalman filter: an alternative to ensemble square root filters, Tellus A, 60, 361–371, 2008. a
Sandery, P. and Sakov, P.: Ocean forecasting of mesoscale features can deteriorate by increasing model resolution towards the submesoscale, Nat. Commun., 8, 1566, https://doi.org/10.1038/s41467017015950, 2017. a
Santana, R., Macdonald, H., O'Callaghan, J., Powell, B., Wakes, S., and H. Suanda, S.: Data assimilation sensitivity experiments in the East Auckland Current system using 4DVar, Geosci. Model Dev., 16, 3675–3698, https://doi.org/10.5194/gmd1636752023, 2023. a
Schaeffer, A. and Roughan, M.: Influence of a western boundary current on shelf dynamics and upwelling from repeat glider deployments, Geophys. Res. Lett., 42, 121–128, 2015. a
Schaeffer, A., Roughan, M., and Wood, J. E.: Observed bottom boundary layer transport and uplift on the continental shelf adjacent to a western boundary current, J. Geophys. Res.Oceans, 119, 4922–4939, https://doi.org/10.1002/2013JC009735, 2014. a
Schaeffer, A., Gramoulle, A., Roughan, M., and Mantovanelli, A.: Characterizing frontal eddies along the E ast A ustralian C urrent from HF radar observations, J. Geophys. Res.Oceans, 122, 3964–3980, 2017. a
Shchepetkin, A. F. and McWilliams, J. C.: The regional oceanic modeling system (ROMS): a splitexplicit, freesurface, topographyfollowingcoordinate oceanic model, Ocean Model., 9, 347–404, 2005. a
Siripatana, A., Kerry, C., Roughan, M., Souza, J. M. A., and Keating, S.: Assessing the impact of nontraditional ocean observations for prediction of the East Australian Current, J. Geophys. Res.Oceans, 125, e2020JC016580, https://doi.org/10.1029/2020JC016580, 2020. a, b, c, d, e
Skamarock, W. C.: Evaluating mesoscale NWP models using kinetic energy spectra, Mon. Weather Rev., 132, 3019–3032, 2004. a
Sloyan, B. M., Ridgway, K. R., and Cowley, R.: The East Australian Current and Property Transport at 27° S from 2012 to 2013, J. Phys. Oceanogr., 46, 993–1008, https://doi.org/10.1175/JPOD150052.1, 2016. a
Souza, J., Powell, B. S., CastilloTrujillo, A. C., and Flament, P.: The Vorticity Balance of the Ocean Surface in Hawaii from a Regional Reanalysis, J. Phys. Oceanogr., 45, 424–440, 2014. a
Stammer, D.: Global Characteristics of Ocean Variability Estimated from Regional TOPEX/ POSEIDON Altimeter Measurements, J. Phys. Oceanogr., 27, 1743–1769, 1997. a
Weaver, A. and Courtier, P.: Correlation modelling on the sphere using generalized diffusion equation., Q. J. Roy. Meteor. Soc., 127, 1815–1846, 2001. a, b
Wilkin, J., Levin, J., Moore, A., Arango, H., López, A., and Hunter, E.: A dataassimilative model reanalysis of the US Mid Atlantic Bight and Gulf of Maine: Configuration and comparison to observations and global ocean models, Prog. Oceanogr., 209, 102919, https://doi.org/10.1016/j.pocean.2022.102919, 2022. a
Xu, L.: 4DVar Data Assimilation for Navy Mesoscale NWP, https://apps.dtic.mil/sti/citations/ADA598257 (last access: 6 March 2024), 2013. a
Xu, Y. and Fu, L.L.: Global variability of the wavenumber spectrum of oceanic mesoscale turbulence, J. Phys. Oceanogr., 41, 802–809, 2011. a
ZavalaGaray, J., Wilkin, J. L., and Arango, H. G.: Predictability of mesoscale variability in the East Australian Current given strongconstraint data assimilation, J. Phys. Oceanogr., 42, 1402–1420, 2012. a
 Abstract
 Highlights
 Introduction
 Model and data assimilation system configuration
 System performance: assessing predictive skill
 Comparisons in model space
 Conclusions
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Highlights
 Introduction
 Model and data assimilation system configuration
 System performance: assessing predictive skill
 Comparisons in model space
 Conclusions
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References