An ensemble Kalman filter system with the Stony Brook Parallel Ocean Model v1.0

This study develops an ensemble Kalman filter (EnKF)-based regional ocean data assimilation system in which the local ensemble transform Kalman filter (LETKF) is implemented with version 1.0 of the Stony Brook Parallel Ocean Model (sbPOM) to assimilate satellite and in situ observations at a daily frequency. A series of sensitivity experiments are performed with various settings of the incremental analysis update (IAU) and covariance inflation methods, for which the relaxation-to-prior perturbations and spread (RTPP and RTPS, respectively) and multiplicative inflation (MULT) are considered. We evaluate the geostrophic balance and the analysis accuracy compared with the control experiment in which the IAU and covariance inflation are not applied. The results show that the IAU improves the geostrophic balance, degrades the accuracy, and reduces the ensemble spread, and that the RTPP and RTPS have the opposite effect. The experiment using a combination of the IAU and RTPP results in a significant improvement for both balance and analysis accuracy when the RTPP parameter is 0.8–0.9. The combination of the IAU and RTPS improves the balance when the RTPS parameter is ≤0.8 and increases the analysis accuracy for parameter values between 1.0 and 1.1, but the balance and analysis accuracy are not improved significantly at the same time. The experiments with MULT inflating the forecast ensemble spread by 5 % do not demonstrate sufficient skill in maintaining the balance and reproducing the surface flow field regardless of whether the IAU is applied or not. The 11 d ensemble forecast experiments show consistent results. Therefore, the combination of the IAU and RTPP with a parameter value of 0.8–0.9 is found to be the best setting for the EnKF-based ocean data assimilation system.

1 Introduction

The ensemble Kalman filter (EnKF; Evensen, 1994, 2003) estimates optimal analyses using model forecasts and observations with their error covariance. The EnKF is advantageous in that it includes flow-dependent forecast errors from an ensemble of model forecasts and is relatively easy to implement with various models. Therefore, various EnKF-based ocean data assimilation systems have been developed thus far (see Table 1).

The number of available observations has increased dramatically with enhanced observations of temperature and salinity in the ocean interior by Argo profiling floats as well as measurements of sea surface temperature, salinity, and height (SST, SSS, and SSH, respectively) by satellites. The Himawari-8 geostationary satellite (Bessho et al., 2016; Kurihara et al., 2016) has an infrared sensor that has been observing SSTs in the Pacific region since July 2015, although there are missing values where cloud obscures the sea surface. Its geostationary orbit and short observation interval allow Himawari-8 to provide better daily coverage within the observation area than a polar-orbiting satellite with a microwave sensor, such as the Global Change Observation Mission-Water (GCOM-W; https://gportal.jaxa.jp/gpr, last access: 11 November 2022), which can capture SSTs even in cloudy regions. Satellite SSS observations by the Soil Moisture and Ocean Salinity (SMOS) mission started in June 2010, and previous studies have demonstrated the positive impacts of these observations, using their ocean data assimilation systems, to better represent the ocean interior structure, such as mixed and barrier layers, low-salinity water caused by river discharge, and prediction of the El Niño–Southern Oscillation (ENSO; Chakraborty et al., 2014; Hackert et al., 2014; Toyoda et al., 2015). The Surface Water and Ocean Topography (SWOT; https://swot.jpl.nasa.gov/, last access: 11 November 2022) satellite that has a new type of altimeter which can observe SSH anomalies (SSHAs) in two dimensions over a 120 km wide swath is scheduled for launch in 2022.

Table 1 Overview of the EnKF-based ocean data assimilation systems developed after 2010. The abbreviations used in the table are as follows: PEODAS – Predictive Ocean Atmosphere Model for Australia (PAOMA) Ensemble Ocean Data Assimilation System, DEnKF – deterministic EnKF (Sakov and Oke, 2008), LETKF – local ensemble transform Kalman filter (Hunt et al., 2007), EAKF – ensemble adjustment Kalman filter (Anderson, 2001), LESKTF – local error subspace Kalman transform filter (Nerger et al., 2012), T – temperature, S – salinity, SST – sea surface temperature, SSH – sea surface height, and MULT – multiplicative inflation. Adaptive MULT was proposed by Miyoshi (2011). Dashes are used to indicate no application. “Inflated obs. error” in TOPAZ4 indicates that observation errors are inflated when ensemble analyses are calculated.

To take advantage of such enhanced observations, frequent data assimilation is important. Here, dynamical imbalances in the analysis field may cause an initial shock with high-frequency gravity waves and may degrade the analysis accuracy. He et al. (2020) described the relationship between the assimilation interval and accuracy using an atmospheric data assimilation system. As seen in Table 1, most of the recent ocean data assimilation systems have an assimilation interval longer than 5 d; in particular, 5 d and 7 d assimilation intervals are employed in the existing ocean reanalysis datasets of the Predictive Ocean Atmosphere Model for Australia Ensemble Ocean Data Assimilation System (PEODAS; Yin et al., 2011) and TOPAZ4 (Sakov et al., 2012), respectively. PEODAS assimilates only in situ temperature and salinity data, whereas TOPAZ4 uses all types of observations but with inflation of observation errors. Although the ocean data assimilation systems constructed by Karspeck et al. (2013) and Miyazawa et al. (2012) have short assimilation intervals of 1 and 2 d, respectively, the former assimilates only in situ temperature and salinity data, and the latter conducts a data assimilation experiment for a short period of 20 d because unrealistic fields are detected if the experiment is performed over several months (Yasumasa Miyazawa, 2022, personal communication). Although Brüning et al. (2021) recently established regional data assimilation systems for the North Sea and Baltic Sea at a frequent interval of 12 h, only satellite SSTs are assimilated. Therefore, the existing systems might mitigate the effects of initial shocks by using the longer assimilation interval, inflating observation errors, and reducing the number of assimilated observations. This is also the case for atmosphere–ocean coupled data assimilation systems (e.g., Brune et al., 2015; Chang et al., 2013; Counillon et al., 2016; Tang et al., 2020). To provide accurate analyses in an EnKF-based ocean data assimilation system in which satellite and in situ observations are assimilated at a frequent interval, it is necessary to investigate an optimal setting for both dynamical balance and accuracy.

The incremental analysis update (IAU; Bloom et al., 1996; see Sect. 2.1) has been proposed to reduce noise from high-frequency gravity waves associated with initial shocks. Covariance relaxation methods such as relaxation-to-prior perturbations (RTPP; Zhang et al., 2004) and relaxation-to-prior spread (RTPS; Whitaker and Hamill, 2012) (see Sect. 2.3), in which the analysis ensemble perturbations are relaxed towards the forecast ensemble perturbations, would also mitigate the initial shock (Houtekamer and Zhang, 2016; Ying and Zhang, 2015). In EnKF-based ocean data assimilation systems, the method used to apply the analysis update to the model evolution and the technique used to inflate the ensemble spread could make significant differences for the dynamical balance and accuracy. However, the IAU and RTPP/RTPS have not been widely used in EnKF-based ocean data assimilation systems (Table 1). Therefore, this study aims to develop an EnKF-based ocean data assimilation system with a frequent assimilation interval of 1 d in order to take advantage of frequent satellite observations and to explore the optimal settings by performing sensitivity experiments with various settings of the IAU and covariance inflation methods.

This paper is organized as follows: Sect. 2 describes the data and methods of IAU, RTPP, and other schemes as well as how to evaluate geostrophic balance and accuracy relative to observations; details of the EnKF-based ocean data assimilation system and sensitivity experiment are described in Sect. 3; Sect. 4 presents the results for geostrophic balance and accuracy in the sensitivity experiments; Sect. 5 compares the prescribed multiplicative inflation (MULT) parameter with the sensitivity experiment with RTPP and IAU; and Sect. 6 provides a summary.

2 Data and methods

In this section, we provide details of the methods used to alleviate some of the problems associated with high-frequency assimilation. Section 2.1 presents the IAU designed to cut off noise from high-frequency gravity waves, and Sect. 2.2 describes perturbed boundary conditions. Covariance inflation methods to prevent the underestimation of ensemble-based forecast error covariance by various factors, such as the limited ensemble size and model imperfections, are introduced in Sect. 2.3, and the methods used to evaluate geostrophic balance and accuracy relative to observations are given in Sect. 2.4.

2.1 IAU

In this study, we implement the IAU (Bloom et al., 1996) based on existing ocean data assimilation systems (Balmaseda et al., 2015; Martin et al., 2015). The procedure for one assimilation cycle is as follows: (i) conduct model integration up to the middle of an assimilation window; (ii) assimilate observations within the window and save the analysis increments in temperature, salinity, and horizontal velocity; and (iii) conduct model integration over the assimilation window adding the increments equally distributed to each time step. The IAU reduces noise from high-frequency gravity waves associated with initial shocks, but the computational cost of the model integration is 1.5 times that of the standard method in which the analyses performed at the beginning of the window are used for the model initial conditions. Following Miyazawa et al. (2012), all analysis variables (SSH, temperature, salinity, and horizontal velocities) are used for initial conditions in the standard method. Although there are various IAU methods, the SSH increments are not included in most of the existing ocean data assimilation systems (Table 2 of Martin et al., 2015), mainly because the SSH increments tend to cause initial shocks. Even without the SSH increments, the SSH would be modified properly in response to the temperature and salinity increments. Therefore, we adopt the analysis increments of temperature, salinity, and horizontal velocity except for SSH.

2.2 Perturbed boundary conditions

Following previous studies (Kunii and Miyoshi, 2012; Penny et al., 2013; Torn et al., 2006), atmospheric and lateral boundary conditions are artificially perturbed for each ensemble member. Atmospheric forcing of the ith ensemble member at a time t, w⁽ⁱ⁾ (t), is given by

$$\begin{array}{}\text{(1)}& {\mathbf{w}}^{\left(i\right)}\left(t\right)=\mathbf{w}\left(t\right)+\mathit{\alpha }\left[\mathbf{w}\left(t+\mathit{\delta }{t}_i\right)-{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n\mathbf{w}\left(t+\mathit{\delta }{t}_i\right)\right],\end{array}$$

where w(t) is atmospheric forcing at a time t, α (=0.2) is an arbitrary constant, w(t+δt_i) is atmospheric forcing at the same time as w(t) but in a different year, and n (=100) is the ensemble size. Here, the year in w(t+δt_i) is changed every month. As is clear from Eq. (1), the ensemble mean of the atmospheric forcing w⁽ⁱ⁾(t) is equivalent to w(t).

Lateral boundary conditions for each ensemble member are obtained from a monthly mean global ocean reanalysis dataset for different years. Namely, the ensemble mean of the lateral boundary condition corresponds to a monthly climatology. These perturbed atmospheric and lateral boundary conditions play a role equivalent to additive inflation (Houtekamer and Zhang, 2016).

2.3 Covariance inflation methods

Three covariance inflation methods (MULT, RTPP, and RTPS) are adopted in this study. MULT inflates forecast error covariance P^f by a factor of ρ (>1):

$$\begin{array}{}\text{(2)}& P_{\inf }^{\mathrm{f}}=\mathit{\rho }{P}_{\mathrm{orig}}^{\mathrm{f}},\end{array}$$

where the subscripts inf and orig denote inflated and original (i.e., before inflation), respectively. Both RTPP and RTPS restore the analysis ensemble perturbation towards the forecast ensemble perturbations maintaining the analysis ensemble mean, as represented by

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathbf{X}}_{\inf }^{\mathrm{a}}={\mathit{\alpha }}_{\mathrm{RTPP}}{\mathbf{X}}^{\mathrm{f}}+\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{RTPP}}\right){\mathbf{X}}_{\mathrm{orig}}^{\mathrm{a}}\mathrm{and}\text{(4)}& {\displaystyle }{\mathbf{X}}_{\inf }^{\mathrm{a}\left(i\right)}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{RTPS}}{\mathit{\sigma }}^{\mathrm{f}\left(i\right)}+(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{RTPS}}){\mathit{\sigma }}_{\mathrm{orig}}^{\mathrm{a}\left(i\right)}}{{\mathit{\sigma }}_{\mathrm{orig}}^{\mathrm{a}\left(i\right)}}}{\mathbf{X}}_{\mathrm{orig}}^{\mathrm{a}\left(i\right)}.\end{array}$$

Here, $\mathbf{X}[=({\mathit{x}}^{\left(\mathrm{1}\right)}-\stackrel{\mathrm{‾}}{x},\mathrm{\dots },{\mathit{x}}^{\left(i\right)}-\stackrel{\mathrm{‾}}{x},\mathrm{\dots },{\mathit{x}}^{\left(\mathbf{n}\right)}-\stackrel{\mathrm{‾}}{x}\left)\right]$ is the ensemble perturbation matrix whose ith column consists of the perturbations of the ith ensemble member, where x⁽ⁱ⁾ and $\stackrel{\mathrm{‾}}{x}$ are the state vector of the ith ensemble member and ensemble mean; the superscripts a and f denote analysis and forecast; and α_RTPP and α_RTPS are the relaxation parameters in the RTPP and RTPS, respectively. σ⁽ⁱ⁾ is the ensemble spread of the ith variable of state vector x, as represented by

$$\begin{array}{}\text{(5)}& {\mathit{\sigma }}^{\left(i\right)}=\sqrt{{\left(n-\mathrm{1}\right)}^{-\mathrm{1}}{\mathbf{X}}^{\left(i\right)}{\left({\mathbf{X}}^{\left(i\right)}\right)}^T}.\end{array}$$

In the RTPP and RTPS, the relaxation parameters are generally defined between 0 and 1, where α_RTPP=0 and α_RTPS=0 correspond to no inflation, and α_RTPP=1 and α_RTPS=1 correspond to the inflated analysis ensemble spread being equivalent to the forecast ensemble spread. RTPP and RTPS are thought to have side effects in maintaining the dynamic balance (Houtekamer and Zhang, 2016; Ying and Zhang, 2015).

2.4 Validation

2.4.1 Nonlinear balance equation (NBE)

Surface horizontal velocity can be represented as the sum of surface geostrophic and ageostrophic velocities under the geostrophic approximation. Here, the ageostrophic velocity is defined as being caused by the surface wind stress curl except for the vertical geostrophic shear, according to the classical Ekman theory (Cronin and Tozuka, 2016). In this study, the atmospheric field is not included in the model state vector; therefore, there are no differences between the forecast and analysis ageostrophic velocities. Consequently, writing the geostrophic balance equation in terms of analysis increments, we obtain

$$\begin{array}{}\text{(6)}& f\mathit{k}\times \mathit{\delta }\mathbf{u}=-g{\mathbf{\nabla }}_{\mathrm{h}}\mathit{\delta }\mathit{\eta },\end{array}$$

where f is the vertical component of the Coriolis parameter, k is a unit vector in the vertical upward direction, δ is the analysis increment, u is the horizontal velocity at the sea surface, g (=9.8 m s⁻²) is the gravitational acceleration, ${\mathrm{\nabla }}_{\mathrm{h}}=(\partial /\partial x,\partial /\partial y)$ is the horizontal gradient operator, and η is the SSH. By taking $\partial /\partial x$ of the x component of Eq. (6) plus $\partial /\partial y$ of the y component, Eq. (6) can be reduced to the nonlinear balance equation (NBE; Shibuya et al., 2015; Zhang et al., 2001):

$$\begin{array}{}\text{(7)}& -f\mathit{\delta }\mathit{\zeta }+\mathit{\beta }\mathit{\delta }u+g{\mathbf{\nabla }}_{\mathrm{h}}^{\mathbf{2}}\mathit{\delta }\mathit{\eta }=\mathrm{0},\end{array}$$

where $\mathit{\zeta }(=\partial v/\partial x-\partial u/\partial y)$ is the relative vorticity at the sea surface, and $\mathit{\beta }(=\partial f/\partial y)$ is the planetary vorticity gradient. If geostrophic balance is not satisfied in the analysis field, there is an absolute residual of the NBE, Δ NBE:

$$\begin{array}{}\text{(8)}& \mathrm{\Delta }\mathrm{NBE}\equiv \left|-f\mathit{\delta }\mathit{\zeta }+\mathit{\beta }\mathit{\delta }u+g{\mathbf{\nabla }}_{\mathrm{h}}^{\mathrm{2}}\mathit{\delta }\mathit{\eta }\right|,\end{array}$$

where | ⋅ | denotes taking the absolute value. A smaller (larger) Δ NBE indicates more (less) geostrophic balance in the analysis field. Few initial shocks would occur if the analysis increments of SSH and surface horizontal velocity satisfy the geostrophic balance.

2.4.2 Improvement ratio (IR)

To compare the geostrophic balance and accuracy among sensitivity experiments using a statistical method, we calculate improvement ratios (IRs) of area-averaged Δ NBE and root-mean-square deviations (RMSDs) relative to observations as represented by

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathrm{IR}}_{\mathrm{N}}={\displaystyle \frac{{\left(\mathrm{\Delta }\mathrm{NBE}\right)}_{\mathrm{CTL}}-(\mathrm{\Delta }\mathrm{NBE}{)}_{\mathrm{EXP}}}{(\mathrm{\Delta }\mathrm{NBE}{)}_{\mathrm{CTL}}}}\times \mathrm{100},\mathrm{and}\text{(10)}& {\displaystyle }{\mathrm{IR}}_{\mathrm{R}}={\displaystyle \frac{{\left(\mathrm{RMSD}\right)}_{\mathrm{CTL}}-(\mathrm{RMSD}{)}_{\mathrm{EXP}}}{(\mathrm{RMSD}{)}_{\mathrm{CTL}}}}\times \mathrm{100},\end{array}$$

respectively. The subscripts CTL and EXP indicate control and sensitivity experiments, respectively. Significant improvement and degradation of the dynamical balance and accuracy are detected by applying the bootstrap method, where the IRs of the area-averaged Δ NBE and RMSDs are resampled for 10 000 cycles, and a 99 % confidence level is used to detect the significance in all sensitivity assimilation experiments.

2.4.3 Observations

To validate the accuracy of the sensitivity experiments, we use observational gridded SSH and SSHA datasets from Archiving Validation and Interpretation of Satellite Oceanographic data (AVISO; Ducet et al., 2000) with a horizontal resolution of 0.25^∘, in situ surface horizontal velocity from surface drifting buoys of the Global Drifter Program (Elipot et al., 2016), in situ temperature and salinity in the depth range 1–525 m, and horizontal velocity in the 8–36 m depth range at 32.3^∘ N, 144.6^∘ E, south of the Kuroshio Extension (KE) from the Kuroshio Extension Observatory (KEO) buoy (https://www.pmel.noaa.gov/ocs/, last access: 11 November 2022; see Fig. 5a). The mean dynamical ocean topography (MDOT) of the AVISO is estimated from a geoid model, satellite altimetry, and in situ drifter buoy data. The AVISO dataset is not an independent observational dataset because satellite SSHAs are used for the assimilation in this study, whereas the surface drifter and KEO buoys are independent. Although validation in the ocean interior might not be sufficient, this is due to the limitation of available independent observations.

Table 2 Overview of the regional ocean model in the ocean data assimilation system.

3 EnKF-based ocean data assimilation system

3.1 Ocean model

The σ-coordinate regional ocean model used in this study is based on version 1.0 of the Stony Brook Parallel Ocean Model (sbPOM; Jordi and Wang, 2012) and constructed for the northwestern Pacific region (15–50^∘ N, 117–180^∘ E) with a horizontal resolution of 0.25^∘ and 50 σ-layers (Table 2). The bottom topography is derived from ETOPO1, a 1 arcmin global relief model of Earth's surface (Amante and Eakins, 2009). We apply a Gaussian filter with e-folding scales of 200 km to the topography to reduce the pressure gradient errors in σ-coordinate models caused by steep bottom slopes (Mellor et al., 1994) and fulfill the condition $\left|H_{i+\mathrm{1}}-H_i\right|/\left|H_{i+\mathrm{1}}+H_i\right|<\mathrm{0.2}$, where H_i and H_i+1 are bottom topographies at adjacent grids. Monthly (seasonal) temperature and salinity climatologies from the World Ocean Atlas 2018 (WOA18; Locarnini et al., 2019; Zweng et al., 2019) with a horizontal resolution of 1^∘ and 57 (102) layers are used for an initial condition over depths shallower (deeper) than 1500 m. Lateral boundary conditions for temperature, salinity, and horizontal velocity are obtained from version 3.7.2 of Simple Ocean Data Assimilation (SODA; Carton et al., 2018) with a horizontal resolution of 0.5^∘ and 50 layers. Here, to satisfy volume conservation, flow relaxation (Guo et al., 2003) is applied to the horizontal velocity at the lateral boundary. The Japanese 55-year Reanalysis (JRA-55; Kobayashi et al., 2015) with horizontal and temporal resolutions of 1.25^∘ and 6 h, respectively, is adopted for the atmospheric boundary conditions, including air temperature and specific humidity at 2 m, wind velocity at 10 m, shortwave radiation, total cloud fraction, sea level pressure, and precipitation. We also use river discharge from the Japan Aerospace Exploration Agency (JAXA)'s land surface and river simulation system, Today's Earth Global (TE-Global; https://www.eorc.jaxa.jp/water/, last access: 11 November 2022), with horizontal and temporal resolutions of 0.25^∘ and 3 h, respectively. The atmospheric and lateral boundary conditions are perturbed as described in Sect. 2.2, except for the rainfall and river discharge.

The model is driven by wind stresses as well as heat and freshwater fluxes using bulk formulae in which bulk coefficients are estimated from the Coupled Ocean–Atmosphere Response Experiment (COARE), version 3.5, bulk algorithm (Brodeau et al., 2017; Edson et al., 2013). The horizontal diffusivity coefficient is calculated by a Smagorinsky type formulation with a coefficient of 0.1 (Smagorinsky et al., 1965) and is assumed to be one-fifth of the horizontal viscosity coefficient. The vertical diffusivity coefficient is estimated by the Level 2.5 version of Nakanishi and Niino (2009). The model is spun up from 1 January 2011 to 6 July 2015 using the initial condition with no motion. During the spin-up period, simulated temperatures and salinity are nudged towards the monthly and seasonal climatologies from WOA18 with a 90 d timescale to damp northward overshooting of the Kuroshio. We have confirmed that the perturbed boundary conditions substantially increase the ensemble spread even with the nudging (not shown).

Table 3 Overview of data assimilation in the ocean data assimilation system.

3.2 Data assimilation

We implement the three-dimensional local ensemble transform Kalman filter (3D-LETKF; Hunt et al., 2007; Miyoshi and Yamane, 2007) with 100 ensemble members to assimilate the following observations on a 1 d assimilation interval (Table 3): satellite SSTs from Himawari-8 and GCOM-W, SSS from the SMOS (http://www.esa.int/Applications/Observing_the_Earth/SMOS, last access: 11 November 2022) and Soil Moisture Active Passive (SMAP) version 4.3 (Meissner et al., 2018), SSH consisting of satellite SSH anomalies from the Copernicus Marine Environment Monitoring Service (CMEMS; http://marine.copernicus.eu/, last access: 11 November 2022) and MDOT estimated from simulated SSH averaged in 2012–2014, and in situ temperature and salinity from the Global Temperature and Salinity Profile Programme (GTSPP; Sun et al., 2010) and Advanced automatic QC (AQC) Argo Data version 1.2a (https://www.jamstec.go.jp/argo_research/dataset/aqc/index_dataset.html, last access: 11 November 2022). We exclude satellite SSS within 100 km of the coasts, SSH for bottom topography shallower than 200 m, in situ temperature and salinity duplicated between the GTSPP and AQC Argo datasets, and observations without the best quality flags or whose differences from the forecasts are larger than the values in the gross error check in Table 3. Following Miyazawa et al. (2012) and Penny et al. (2013), the localization scales based on a Gaussian function are chosen to be 300 km and 100 m in the horizontal and vertical directions, respectively. An observational error covariance matrix is assumed to be diagonal using the observation errors in Table 3.

3.3 Sensitivity experiments

We conduct sensitivity experiments combining the IAU and covariance inflation methods (no inflation (NO INFL), RTPP, RTPS, and MULT) to investigate their impacts on the geostrophic balance and accuracy. We set the relaxation parameters in the RTPP and RTPS experiments to ${\mathit{\alpha }}_{\mathrm{RTPP}}={\mathit{\alpha }}_{\mathrm{RTPS}}=\mathrm{0.9}$ without the IAU and to ${\mathit{\alpha }}_{\mathrm{RTPP}}={\mathit{\alpha }}_{\mathrm{RTPS}}=\mathrm{0.5},\mathrm{0.6},\mathrm{\dots },\mathrm{1.2}$ with the IAU, and we set the inflation parameter to ρ=1.05² (inflating the forecast ensemble spread by 5 %) in the MULT experiments, regardless of the application of the IAU. In this study, we do not explore all values of the relaxation and inflation parameters because of the limitations of computational resources. Hereafter, we refer to the RTPP experiments implemented with and without the IAU as the RTPP+IAU and RTPP experiments, respectively, and we refer to the RTPP+IAU experiment with a relaxation parameter of 0.5 as the RTPP05+IAU experiment. Kotsuki et al. (2017) indicated that the RTPP and RTPS do not consider the model error explicitly and that the optimal relaxation parameter may be larger than 1.0. Therefore, we perform experiments with a relaxation parameter value >1. To clarify the effects of the IAU and covariance inflation methods, the NO INFL experiment is defined as a control experiment in this study (see Eqs. 9 and 10).

We integrate the LETKF-based ocean data assimilation system from 7 July 2015 at the start date of the Himawari-8 observations to 31 December 2016, applying the SSS nudging with 90 d timescale to damp a surface freshening drift, as in the model spin-up described in Sect. 3.1. Furthermore, we conduct 11 d ensemble forecast experiments initialized on the first day of each month in 2016 by the forecasts from the NO INFL, NO INFL+IAU, RTPP09, RTPP09+IAU, RTPS09, and RTPS09+IAU experiments, with the SSS nudging applied with a 90 d timescale. We estimate Δ NBE from the ensemble analysis increments on days 1 and 16 of each month, the RMSDs from the daily averaged ensemble mean analyses and forecasts, and the ensemble spread from the daily mean ensemble analyses. As described in Sect. 2.4.2, the statistical analyses are applied to IRs of area-averaged Δ NBE and analysis RMSDs in all analysis experiments. The results of the RTPP11+IAU and RTPP12+IAU experiments are not shown because numerical instability developed.

4 Results

4.1 Geostrophic balance

We first compare the geostrophic balance for the various sensitivity experiments using spatiotemporally averaged Δ NBE over the whole system domain for 2016 (Fig. 1). The NO INFL+IAU experiment has the best geostrophic balance with significant improvement relative to the NO INFL experiment; thus, the IAU plays a role in enhancing the balance, probably because the IAU reduces noise of the high-frequency gravity waves associated with initial shocks. This result is consistent with Yan et al. (2014), who demonstrated that the IAU reduces the spurious oscillation of vertical velocity in twin experiments using a relatively idealized EnKF-based ocean data assimilation system. In contrast, as the RTPP09 and RTPS09 experiments show significantly larger Δ NBE than the NO INFL experiment, RTPP and RTPS contribute to breaking the balance. The MULT+IAU and MULT experiments give such a large Δ NBE ($\mathrm{2.11}\times {\mathrm{10}}^{-\mathrm{10}}$ and $\mathrm{5.22}\times {\mathrm{10}}^{-\mathrm{10}}$ s⁻², respectively) that the MULT breaks the balance considerably even if the IAU is applied.

Figure 1 Spatiotemporally averaged Δ NBE over the whole domain for 2016 in the NO INFL (black star), RTPP (red), RTPS (blue), NO INFL+IAU (gray), RTPP+IAU (orange), and RTPS+IAU (cyan) experiments as a function of the relaxation parameters. Open circles and triangles indicate significant improvement and degradation relative to the NO INFL experiment at a 99 % confidence level, respectively, and closed circles and triangles denote improvement and degradation relative to the NO INFL experiment with no significant differences. The RTPS12+IAU, MULT+IAU, and MULT experiments show significant degradation with an averaged Δ NBE of $\mathrm{2.94}\times {\mathrm{10}}^{-\mathrm{10}}$, $\mathrm{2.11}\times {\mathrm{10}}^{-\mathrm{10}}$, and $\mathrm{5.22}\times {\mathrm{10}}^{-\mathrm{10}}$ s⁻², respectively (not shown). The RTPP+IAU experiments for the relaxation parameters of α_RTPP≥1.1 are not shown because numerical instability developed.

The RTPP+IAU and RTPS+IAU experiments provide significant improvement when the relaxation parameters are α_RTPP≤0.9 and α_RTPS≤0.8, respectively. The RTPP11+IAU experiment becomes numerically unstable in December 2015, and the RTPS11+IAU experiment also significantly degrades the balance; thus, relaxation parameters larger than 1.0 do not appear to be appropriate for the EnKF-based ocean data assimilation system. The combinations of the IAU and RTPP/RTPS, in which the relaxation parameters are set to α_RTPP≤0.9 and α_RTPS≤0.8, appear to maintain geostrophic balance, likely because the IAU counteracts the RTPP/RTPS by improving the balance.

To investigate spatial characteristics of the geostrophic balance, Δ NBE is temporally averaged over the whole year 2016 (Fig. 2). Here, the RTPP09+IAU and RTPS11+IAU experiments are shown from the RTPP+IAU and RTPS+IAU experiments because they have the best accuracy, as seen in Sect. 4.2. The NO INFL, RTPP09, and RTPS09 experiments produce less-balanced fields in the midlatitude region, especially around the KE (Fig. 2a, c, e). In the RTPS11+IAU experiment, the balance is also lost in higher-latitude regions (Fig. 2f). In the MULT and MULT+IAU experiments, there are almost no balanced regions with Δ NBE smaller than $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{10}}$ s⁻² (not shown). The NO INFL+IAU and RTPP09+IAU experiments show substantial improvement around the KE region, although a relatively large Δ NBE remains along the KE in the RTPP09+IAU experiment (Fig. 2b, d). Thus, in general, although the balance in the analysis field is not maintained around the KE region, this imbalance is substantially reduced in the NO INFL+IAU and RTPP09+IAU experiments.

Figure 2 Δ NBE (colors) and SSH (white contours) averaged over 2016 in the (a) NO INFL, (b) NO INFL+IAU, (c) RTPP09, (d) RTPP09+IAU, (e) RTPS09, and (f) RTPS11+IAU experiments. Thin (thick) contour intervals are 0.2 m (1.0 m).

4.2 Accuracy

4.2.1 Surface flow field

We evaluate the accuracy of the surface flow field in the sensitivity experiments, calculating the analysis RMSDs relative to the AVISO observational SSH and SSHA gridded datasets as well as surface zonal and meridional velocity from the drifter buoys. We also estimate the ensemble spread in observational space. As described in Sect. 2.4.3, the AVISO dataset is not independent, as it uses satellite SSHAs assimilated in our system, whereas the drifter buoys are independent. The different results from the SSH and SSHA RMSDs are caused by the different MDOT between the AVISO dataset and the system, as described in Sects. 2.4.3 and 3.2, respectively. The analysis RMSDs and ensemble spreads are averaged over the whole domain for 2016 in the SSH and SSHA fields (Fig. 3) and the surface zonal and meridional velocity fields (Fig. 4). Compared with the NO INFL experiment, the NO INFL+IAU experiment has significantly larger RMSDs and smaller ensemble spreads in most of the variables, whereas the RTPP09 and RTPS09 experiments show significantly smaller RMSDs and larger spreads (Figs. 3, 4). This indicates that the IAU has a significant effect on reducing the accuracy in the surface flow field because the relatively small ensemble spread leads to small analysis increments and because the IAU does not use the SSH analysis increments. However, the small analysis increments result in a better dynamical balance, as shown in Sect. 4.1. The result is consistent with Yan et al. (2014), who demonstrated that the IAU degrades the accuracy of SSH, temperature, and horizontal velocities using twin experiments. In contrast, the RTPP and RTPS lead to significant improvement by inflating the ensemble spread. The large analysis increments caused by the large ensemble spread might reduce the dynamical balance. The MULT and MULT+IAU experiments yield poor accuracy and very large ensemble spreads in the flow fields; for example, the averaged SSH RMSDs of 0.22 and 0.24 m and the averaged SSH and SSHA ensemble spreads of 0.41 and 0.74 m, respectively. Thus, the MULT does not have sufficient skill in reproducing the flow field.

Figure 3 As in Fig. 1 but for the analysis RMSDs of (a) SSH and (b) SSHA relative to the AVISO dataset. Panel (c) shows the spatiotemporally averaged ensemble spreads of SSH and SSHA over the whole domain for 2016 in observational space (circles). The RMSDs of SSH and SSHA in the RTPS12+IAU experiment are 0.164 and 0.137 m, respectively (not shown).

In both RTPP+IAU and RTPS+IAU experiments, the ensemble spreads are increased in all of the variables for the larger relaxation parameters (Figs. 3c, 4c, d). It appears that the larger relaxation parameters maintain the large ensemble spread induced by the perturbed boundary conditions. In the RTPP+IAU experiments, the accuracy of SSH and SSHA is the highest for α_RTPP=0.8, although there is no significant improvement relative to the NO INFL experiment (Fig. 3a, b). The accuracy in both zonal and meridional velocity improves with larger relaxation parameter, and significantly improves for α_RTPP=0.8–1.0 (Fig. 4a, b). Consequently, α_RTPP=0.8–0.9 in the RTPP+IAU experiment may be appropriate to represent the flow field more accurately.

In the RTPS+IAU experiment, the accuracy of the SSH, SSHA, and horizontal velocity tends to improve as the relaxation parameter increases, and then significant degradation suddenly occurs for α_RTPS=1.2 (Figs. 3a, b; 4a, b). For α_RTPS=1.1, the RTPS+IAU experiment has the best accuracy for the SSH and SSHA but significant improvement only in the SSHA (Fig. 3a, b). For α_RTPS=1.0, the accuracy in both zonal and meridional velocity is significantly higher (Fig. 4a, b). Therefore, α_RTPS=1.0–1.1 seems to be the best among the RTPS+IAU experiments. We note that the accuracy of the SSH and SSHA in the RTPP+IAU and RTPS+IAU experiments does not surpass the RTPP09 and RTPS09 experiments, probably because the IAU method does not use the SSH analyses. Furthermore, the comparison between the RTPP+IAU and RTPS+IAU experiments suggests that the combination of the IAU and RTPP has higher skill in reproducing the flow field.

Figure 4 As in Fig. 1 but for the analysis RMSDs of the surface (a) zonal and (b) meridional velocity relative to the drifter buoys as well as the ensemble spreads of the surface (c) zonal and (d) meridional velocity. The RMSDs of surface zonal and meridional velocity in the RTPS12+IAU experiment are 0.293 and 0.277 m s⁻¹, respectively (not shown). The RMSD in panel (b) and ensemble spreads in panels (c) and (d) in the RTPP09 experiment are slightly offset for visualization.

To examine the spatial features of the analysis accuracy and ensemble spread, the analysis RMSDs and ensemble spreads in the SSHA are also averaged over 2016 (Figs. 5 and 6, respectively). In most experiments, large RMSDs and ensemble spreads are distributed around the KE region, where there are abundant fronts and eddies. Compared with the NO INFL experiment, the ensemble spreads become smaller in the midlatitude region in the NO INFL+IAU experiment; thus, the accuracy around the KE region is degraded. The RTPP09, RTPS09, RTPP09+IAU, and RTPS11+IAU experiments show larger ensemble spreads, leading to improvement of the accuracy around the KE region. However, the larger ensemble spread is also seen in the subtropical region in the RTPS11+IAU experiment. This does not seem reasonable because a free ensemble experiment does not demonstrate such spread even if the perturbed atmospheric and lateral boundary conditions are applied (not shown).

Figure 5 As in Fig. 3 but for the analysis RMSDs relative to the SSHA from the AVISO dataset (color). The black star in panel (a) indicates the KEO buoy location (32.3^∘ N, 144.6^∘ E).

Figure 6 As in Fig. 3 but for the SSHA ensemble spreads.

To investigate the forecast accuracy, we calculate the spatiotemporally averaged forecast RMSDs of the 11 d ensemble forecast experiments for each month in 2016 (i.e., a total of 12 cases) relative to the AVISO and drifter buoys, and the 12 cases are averaged to obtain the forecast RMSDs over 2016 (Fig. 7). As shown in Figs. 3, 4, and 7, the results of the forecast RMSDs generally agree with those of the analysis RMSDs, except for the RTPP09+IAU and RTPS09+IAU experiments showing smaller forecast SSHA RMSDs than the NO INFL experiment. Overall, the combination of the IAU and RTPP09 seems to be the most suitable for not only constructing analysis products but also conducting ensemble forecasts.

Figure 7 Spatiotemporally averaged RMSDs of the 11 d ensemble forecast in the NO INFL (black), RTPP09 (red), RTPS09 (blue), NO INFL+IAU (gray), RTPP09+IAU (orange), and RTPS09+IAU (cyan) experiments relative to (a) SSH and (b) SSHA from the AVISO dataset and surface (c) zonal and (d) meridional velocities from the drifter buoys. The RMSDs in the RTPS09 in panels (b), (c), and (d) are slightly offset for visualization.

4.2.2 The KEO buoy

We also calculate the analysis and forecast RMSDs relative to independent observations of temperature, salinity, and horizontal velocity from the KEO buoy located south of the KE (Fig. 5a). Here, only the temperature and salinity results are shown because there is basically no improvement in the horizontal velocity. There is almost no difference between the NO INFL+IAU and NO INFL experiments in the temperature analysis accuracy, whereas the salinity analysis accuracy is significantly degraded around 0–200 m depth in the NO INFL+IAU experiment (Fig. 8). Therefore, the IAU may reduce the analysis accuracy, although this is not as obvious as for the flow field shown in Sect. 4.2.1. The RTPP and RTPS experiments give significantly better analysis accuracy than the NO INFL experiment for both temperature and salinity; thus, the RTPP and RTPS play a role in enhancing the analysis accuracy. These results are qualitatively the same as the forecast accuracy (Fig. 9).

Figure 8 Analysis RMSDs of (a) temperature and (b) salinity relative to the KEO buoy averaged over 2016 in the NO INFL (black star), RTPP09 (red), RTPS09 (blue), NO INFL+IAU (gray), RTPP09+IAU (orange), and RTPS11+IAU (cyan) experiments. Open circles and triangles denote significant improvement and degradation relative to the NO INFL experiment at a 99 % confidence level, respectively. Closed circles and triangles indicate improvement and degradation with no significant differences, respectively.

When the relaxation parameter α_RTPP is 0.7–1.0 in the RTPP+IAU experiment, the temperature analysis accuracy is significantly enhanced around 200–500 m depth, although there is slight degradation around 50–150 m depth (Fig. 10a). For the parameter values in that range, the salinity analysis accuracy is also significantly improved at almost all depths (Fig. 10c). As the temperature analysis accuracy is the best at α_RTPP=0.8–0.9 and the salinity analysis accuracy improves as the relaxation parameter increases, the appropriate relaxation parameter would be α_RTPP=0.8–0.9 in the RTPP+IAU experiment.

Figure 9 As in Fig. 8 but for the forecast RMSDs of the 11 d ensemble forecast experiments. We note that the results of the RTPS09+IAU experiment are shown here, whereas the analysis RMSDs of the RTPS11+IAU experiment are shown in Fig. 8; thus, the relaxation parameters are different in the RTPS+IAU experiments.

Figure 10 Temperature analysis RMSDs (black contours) and IRs (color shading and white contours) between the KEO buoy and the (a) RTPP+IAU and (b) RTPS+IAU experiments averaged over 2016. Panels (c) and (d) are the same as panels (a) and (b) but for salinity. Open circles and triangles indicate significant improvement and degradation relative to the NO INFL experiment at a 99 % confidence level, respectively. Thin (thick) black contour intervals are 0.2 (1.0) ^∘C in panels (a) and (b), whereas they are 0.1 (0.2) in panels (c) and (d); thin (thick) white contour intervals are 10 % (100 %).

In the RTPS+IAU experiments, the temperature analysis accuracy below 200 m depth is significantly improved for α_RTPS=0.8–1.1, whereas that above 200 m depth is significantly degraded for α_RTPS=0.5–0.8 and α_RTPS=1.2 (Fig. 10b). The salinity analysis accuracy improves over almost the whole depth when the relaxation parameter is α_RTPS=1.0–1.2, whereas there is significant degradation around 0–200 m depth when the relaxation parameter is α_RTPS=0.5–0.9 (Fig. 10d). Therefore, a suitable relaxation parameter is α_RTPS=1.0–1.1 in the RTPS+IAU experiment.

The RTPP09+IAU and RTPS11+IAU experiments have higher analysis accuracy (Fig. 8), and the RTPP09 and RTPP09+IAU experiments show higher forecast accuracy (Fig. 9) than the other experiments. Therefore, the combination of the IAU and RTPP09 is the most appropriate for the analysis and ensemble forecasts.

Table 4 Schematic summarizing the evaluation of the geostrophic balance and analysis accuracy of the AVISO SSH and SSHA, the surface zonal and meridional velocity from the drifter buoys, and the temperature and salinity at the KEO buoy in the sensitivity experiments. Open circles and crosses indicate improvement and degradation relative to the NO INFL experiment, respectively, and asterisks denote significant improvement and degradation. In rows two to four, symbols and asterisks are used only if both variables have the same results; otherwise, dashes are used to indicate no significant difference from the NO INFL experiment. Parentheses in the RTPP+IAU and RTPS+IAU experiments denote the best relaxation parameter in the second row and the range of the relaxation parameter with significant improvement in the other rows.

5 Comparison of the prescribed MULT parameter with the RTPP09+IAU experiment

To investigate how much the inflation in the RTPP09+IAU experiment corresponds to the MULT parameter, we estimate the MULT parameter ρ_est corresponding to the RTPP09+IAU experiment using the following equation:

$$\begin{array}{}\text{(11)}& {\mathit{\alpha }}_{\mathrm{RTPP}}{\mathbf{X}}^{\mathrm{f}}+\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{RTPP}}\right){\mathbf{X}}_{\mathrm{orig}}^{\mathrm{a}}=\sqrt{{\mathit{\rho }}_{\mathrm{est}}}{\mathbf{X}}^f.\end{array}$$

By multiplying $({\mathbf{X}}^{\mathrm{f}}{)}^T{\left[{\mathbf{X}}^{\mathrm{f}}({\mathbf{X}}^{\mathrm{f}}{)}^T\right]}^{-\mathrm{1}}$ from the right-hand side of Eq. (11),

$$\begin{array}{}\text{(12)}& {\mathit{\rho }}_{\mathrm{est}}\mathbf{I}={\left\{{\mathit{\alpha }}_{\mathrm{RTPP}}I+\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{RTPP}}\right){\mathbf{X}}_{\mathrm{orig}}^{\mathrm{a}}({\mathbf{X}}^{\mathrm{f}}{)}^T{\left[{\mathbf{X}}^{\mathrm{f}}({\mathbf{X}}^{\mathrm{f}}{)}^T\right]}^{-\mathrm{1}}\right\}}^{\mathrm{2}},\end{array}$$

where I denotes the identity matrix. In scalar format, the estimated parameter ${\mathit{\rho }}_{\mathrm{est}}^{\left(i\right)}$ at the ith variable might be represented as

$$\begin{array}{}\text{(13)}& {\mathit{\rho }}_{\mathrm{est}}^{\left(i\right)}={\left\{{\mathit{\alpha }}_{\mathrm{RTPP}}+{\displaystyle \frac{\mathrm{1}-{\mathit{\alpha }}_{\mathrm{RTPP}}}{n-\mathrm{1}}}{\displaystyle \frac{X_{\mathrm{orig}}^{\mathrm{a}\left(i\right)}(X^{\mathrm{f}\left(i\right)}{)}^T}{({\mathit{\sigma }}^{f\left(i\right)}{)}^{\mathrm{2}}}}\right\}}^{\mathrm{2}}.\end{array}$$

Using the outputs from the RTPP09+IAU experiment, we calculate the estimated MULT parameter for the SST, SSS, and SSH fields (Fig. 11). The estimated MULT parameter is large around the midlatitude region, especially around the KE region. The estimated MULT parameters averaged over the whole domain and analysis period are 1.08 (1.11) for the SST and SSS (SSH) fields, and these values correspond well to the prescribed MULT parameter $\mathit{\rho }={\mathrm{1.05}}^{\mathrm{2}}\approx \mathrm{1.10}$.

Figure 11 Estimated MULT parameters (Eq. 11) averaged over 2016 for (a) SST, (b) SSS, and (c) SSH fields using the outputs from the RTPP09+IAU experiment. Right bottom values indicate spatiotemporally averaged estimated MULT parameters. Thin (Thick) counter intervals are 0.02 (0.1).

As shown in Fig. 11, the MULT parameter might have spatial dependency; therefore, adaptive MULT (Miyoshi, 2011) may be useful. However, Ohishi et al. (2022) demonstrated that adaptive observation error inflation (AOEI; Minamide and Zhang, 2017; Zhang et al., 2016), with opposite effects to the adaptive MULT, significantly improves the dynamical balance and accuracy of the temperature, salinity, and surface horizontal velocities. This is because the AOEI suppresses the erroneous temperature and salinity analysis increments associated with the representation errors around the KE region, which result in strong vertical salinity diffusion through weakening density stratification and therefore degrade the low-salinity structure in the intermediate layer. This implies that the adaptive MULT would increase the analysis increments and degrade the dynamical balance and accuracy. Therefore, it is difficult to find an appropriate MULT parameter.

6 Summary

In this study, we have developed an EnKF-based ocean data assimilation system with an assimilation interval of 1 d to take advantage of frequent satellite observations; moreover, we have conducted sensitivity experiments to explore the best combination of the IAU and covariance inflation methods by evaluating the geostrophic balance and analysis accuracy. Table 4 summarizes the overall evaluation in this study. The IAU and RTPP/RTPS have opposite effects to each other; namely, the IAU improves the balance but degrades the accuracy, reducing the ensemble spread, whereas the RTPP and RTPS degrade the balance and improve the accuracy by inflating the ensemble spread. Large RTPP and RTPS parameters maintain large ensemble spread inflated by the perturbed boundary conditions, and the resulting large analysis increments degrade the balance but improve the accuracy. The RTPP+IAU experiment provides significantly better balance for relaxation parameters of α_RTPP≤0.9 as well as better accuracy when the relaxation parameter is α_RTPP=0.8–1.0. Therefore, this study demonstrates that the appropriate parameter is α_RTPP=0.8–0.9 when the IAU and RTPP are combined. In contrast, the RTPS+IAU experiment does not significantly improve the balance and accuracy at the same time, as the balance is significantly better for a relaxation parameter of α_RTPS≤0.8, whereas the accuracy is significantly higher when the relaxation parameter is α_RTPS=1.0–1.1. Therefore, this study demonstrates that the combination of the IAU and RTPP with a relaxation parameter of α_RTPP=0.8–0.9 is the most suitable for the EnKF-based ocean data assimilation system. The 11 d ensemble forecast experiments show consistent results of forecast accuracy with the analysis accuracy.

In the combination of the IAU and RTPP, the large relaxation parameter of α_RTPP=0.8–0.9 maintains the ensemble spread induced by perturbed boundary conditions and leads to the improvement of the analysis accuracy but the degradation of the dynamical balance; concurrently, the IAU improves the degradation of the dynamical balance by the RTPP. As a result, this would lead to further improvement of the forecast and analysis accuracy by reducing the initial shocks in frequent data assimilation. Compared with the RTPS (RTPS+IAU) experiments, the RTPP (RTPP+IAU) experiments show better balance and result in smaller initial shocks. As a result, the combination of the IAU and RTPP leads to better accuracy than that of the IAU and RTPS.

The MULT with a 5 % inflation of the forecast ensemble spread does not have sufficient skill in maintaining the balance and accurately reproducing the flow field, regardless of whether or not the IAU is applied. Although it is difficult to find an appropriate MULT parameter, as described in Sect. 5, it might be possible that MULT produces analyses with good balance and accuracy by tuning the inflation parameter. However, as the computational cost of tuning the parameters in all covariance inflation methods is high, this study focuses on the combination of the RTPP/RTPS and IAU with good balance and accuracy. This system still contains other tuning parameters in the perturbed atmospheric forcing, ensemble size, localization scale, and observation errors. We note that the suitable RTPP parameter in the RTPP+IAU experiment would be different depending on those parameter settings. Further experiments are required to determine the best settings for a given computational resource, and we will address this issue in future studies.

The results of this study would also be useful for constructing EnKF-based data assimilation systems in other fields in which gravity waves have substantial impacts. Furthermore, this study may help improve the accuracy of existing EnKF-based data assimilation systems. Table 1 shows that there are no eddy-resolving EnKF-based ocean reanalysis datasets in the Pacific region. We are now planning to construct such analysis datasets and real-time ensemble prediction systems.