Mitigation of model bias influences on wave data assimilation with multiple assimilation systems using WaveWatch III v5.16 and SWAN v41.20

High-quality wave prediction with a numerical wave model is of societal value. To initialize the wave model, wave data assimilation (WDA) is necessary to combine the model and observations. Due to imperfect numerical schemes and approximated physical processes, a wave model is always biased in relation to the real world. In this study, two assimilation systems are first developed using two nearly independent wave models; then, “perfect” and “biased” assimilation frameworks based on the two assimilation systems are designed to reveal the uncertainties of WDA. A series of biased assimilation experiments is conducted to systematically examine the adverse impact of model bias on WDA. A statistical approach based on the results from multiple assimilation systems is explored to carry out bias correction, by which the final wave analysis is significantly improved with the merits of individual assimilation systems. The framework with multiple assimilation systems provides an effective platform to improve wave analyses and predictions and help identify model deficits, thereby improving the model.

1 Introduction

Ocean waves, referring to the ocean surface gravity waves driven by wind, are important physical processes in the study of multiscale coupled systems. Many studies show that ocean waves are necessary for upper-ocean mixing processes, whether in small-scale coastal simulations or large-scale global climate simulations (e.g., Babanin et al., 2009; Huang and Qiao, 2010; Qiao et al., 2004, 2010). The existence of ocean waves can modify the structures of both atmospheric and marine boundary layers by providing sea surface roughness, wave-induced bottom stress, breaking-wave-induced mixing and so on, which ultimately influence air–sea momentum and heat exchange. Therefore, ocean waves are an important component in atmosphere–ocean interaction flux processes (e.g., Chen et al., 2007; Doyle, 2002; Liu et al., 2011; Warner et al., 2010). In addition, the study of ocean waves can reduce and prevent marine disasters and provide guidance for the development of the social economy (e.g., Folley and Whittaker, 2009; Rusu, 2015; Wei et al., 2017). Thus, studying ocean waves is of great scientific and social significance.

At present, ocean wave observational techniques are constantly being improved (e.g., Daniel et al., 2011; Hisaki, 2005). Except for traditional buoy observations (e.g., Mitsuyasu et al., 1980; Rapizo et al., 2015; Walsh et al., 1989), satellites can provide much near-real-time wave observational information, which is beneficial for understanding the state of ocean waves (e.g., Gommenginger et al., 2003; Lzaguirre et al., 2011; Queffeulou, 2004). However, observations always represent scattered samples in time and space in the real world and therefore do not represent the complete three-dimensional structure and temporal evolution of real-world waves.

Numerical wave models are a powerful tool for studying the physical processes of ocean waves and predicting future wave states. Following the development of the previous two generations, third-generation wave models, such as WAve Modeling (WAM) (WAMDI Group, 1988), WaveWatch III (WW3) (Tolman, 1991), Simulating Waves Nearshore (SWAN) (Booij et al., 1999), and MArine Science and Numerical Modeling (MASNUM) (Yang et al., 2005), integrate the spectral action balance equation describing the two-dimensional ocean wave spectrum evolution without additional ad hoc assumptions regarding the spectral shape, and these third-generation models are more robust for arbitrary wind fields than previous models. However, there are generally three error sources in wave models. One error source is from an incomplete understanding of the physical processes, approximate expressions of the numerical discretization schemes and so on, which causes systematic errors that are usually referred to as wave model bias. The second error source is due to inaccurate wind forcings of wave models. The third error source is from the initial-condition uncertainties, which can grow due to the nonlinearity of the model equations during model forwarding. In this sense, the model-simulated waves do not represent the real world either.

Given the scattering nature of observational information and the approximate characteristics of wave modeling, wave model data assimilation (WDA) is necessary to combine the advantages of both the model and observations. WDA optimizes the model initial conditions to produce more accurate wave forecasts and produces a more accurate evolution of three-dimensional wave states to elucidate the underlying mechanisms; this approach dates back to the 1980s (e.g., Esteva, 1988; Janssen et al., 1989). Since then, many advanced WDA methods have been developed (e.g., Abdalla et al., 2013; Bauer et al., 1996; Greenslade and Young, 2004; Jesus and Cavaleri, 2015; Lionello et al., 1992; Sun et al., 2017; Vorrips et al., 1999), and their applications have been assessed (e.g., Francis and Stratton, 1990; Heras et al., 1994; Stopa and Cheung, 2014). Furthermore, various observation types, such as buoy, radar and satellite, have been applied to WDA (e.g., Bhatt et al., 2005; Breivik et al., 1998; Feng et al., 2006; Greenslade, 2001; Hasselmann et al., 1997; Qi and Cao, 2016; Voorrips, 1999; Waters et al., 2013), and the wave forecasts have also been directly addressed (e.g., Almeida et al., 2016; Emmanouil et al., 2012; Lionello et al., 1995; Qi and Fan, 2013; Sannasiraj et al., 2006; Voorrips, 1999; Wang and Yu, 2009; Zhang et al., 2003).

Due to the approximate nature of the numerical discretization and physical processes, a systematic difference between a model and the real world (i.e., model bias) exists. As noted by Zhang et al. (2012), since model bias is not well defined in observational space, the influence of model bias on data assimilation is a challenging research topic. Alternatively, one can simulate model bias using a pair of models and study the adverse impacts on data assimilation. Inspired by previous work (e.g., Dee, 2005; Zhang et al., 2012), here, we use a simple data assimilation scheme with two wave models (WW3 and SWAN) to explore the influences of different error sources on WDA. The adverse impacts of wind forcing errors and initial-condition uncertainties as well as wave model bias on WDA are studied first, and then two simple statistical methods for bias correction are developed to mitigate assimilation errors and improve wave analysis.

This paper is organized as follows. After the introduction, the methodology is presented in Sect. 2, including a brief description of the employed models and observations, the development of the two WDA systems using the WW3 and SWAN models, as well as the design of experiments throughout the study. Section 3 presents the model bias analysis and the adverse impacts of model bias on WDA. In Sect. 4, the method used to mitigate model bias influences on wave assimilation is explored. Finally, the discussion and conclusion are given in Sect. 5.

2 Methodology

2.1 Models and data

2.1.1 Three models

In the wave models, the variance spectrum or energy density E(σ,θ) is a quantity that represents the wave energy distribution in the radian frequency (σ) and propagation direction (θ). Without ambient ocean currents, the variance or energy of a wave package is conserved. However, if the current is involved, due to the work done by the current on the mean momentum transfer of waves (Longuet-Higgins and Stewart, 1961, 1962), the energy of a spectral component is no longer conserved. In general, an action density spectrum defined as $N=E/\mathit{\sigma }$ is considered within the models. Then, the governing equation of the wave model can be written as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial N}{\partial t}}+{\mathrm{\nabla }}_{\mathit{x}}\cdot \left({\mathit{c}}_{\text{g}}N\right)+{\displaystyle \frac{\partial c_{\mathit{\sigma }}N}{\partial \mathit{\sigma }}}+{\displaystyle \frac{\partial c_{\mathit{\theta }}N}{\partial \mathit{\theta }}}={\displaystyle \frac{S_{\text{tot}}}{\mathit{\sigma }}}.\end{array}$$

The left-hand side of this equation is the kinematic process during wave propagation. The second term describes the wave energy propagation in two-dimensional geographical space denoted by x. The c_g is the group velocity that follows the dispersion relation. The third term represents the effect of shifting the radian frequency due to variation in depth. The fourth term represents the depth-induced refraction. c_σ and c_θ are wave velocities in the frequency σ and direction θ, respectively. On the right-hand side, S_tot is the nonconservative source–sink term representing all physical processes that generate, dissipate or redistribute wave energy. Typically, there are three important physical processes that contribute to S_tot, which include the atmosphere–wave interaction, nonlinear wave–wave interaction and wave–ocean interaction. In a shallow-water case, additional processes must be considered, such as wave–bottom interaction, depth-induced breaking and triad wave–wave interaction.

In this study, we use three advanced third-generation spectrum models – WAM (Cycle 4.5.4) (WAMDI Group, 1988), WaveWatch III (WW3, version 5.16) (Tolman, 1991) and SWAN (version 41.20) (Booij et al., 1999) – which have improved a lot from their predecessors, such as numerical and physical approaches. These wave models have the same form of governing equation; however their numerical implement processes are different due to different considerations. For example, SWAN is more focused on wave propagation processes in shallow water, while WAM and WW3 pay more attention to deep water. For more comparisons, please see Sect. 2.2.

2.1.2 Model configurations

Three wave models use two-dimensional spectral space containing 29 frequencies that cover from 0.035 to 0.555 Hz with a logarithmic distribution and 24 equidistant directions. The geographic space is from 180^∘ W to 180^∘ E in the zonal direction and 75^∘ S to 75^∘ N in the meridional direction with a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid resolution. The topography in this study is taken from the high-resolution ETOPO1 dataset provided by NOAA (website: https://www.ngdc.noaa.gov/mgg/global/, last access: 3 March 2020). The wind forcing has two sources, both of which have 6 h time intervals. The first dataset is the ERA-Interim reanalysis from the European Centre for Medium-Range Weather Forecasts (ECMWF), with a resolution of $\mathrm{0.75}{}^{\circ }\times \mathrm{0.75}{}^{\circ }$ (http://apps.ecmwf.int/datasets/data/interim-full-daily/, last access: 3 March 2020). The second dataset is the CFSRv2 dataset from NCEP, with a resolution of $\mathrm{0.205}{}^{\circ }\left(\text{longitude}\right)\times \mathrm{0.204}{}^{\circ }\left(\text{latitude}\right)$ (https://rda.ucar.edu/datasets/ds094.1/, last access: 3 March 2020). The time step of all three models is 15 min. All relevant parameters above are set to be identical for every wave model (WW3, SWAN and WAM) in this study.

2.1.3 Data

The AVISO (Archiving, Validation and Interpolation of Satellite Oceanographic) data (https://www.aviso.altimetry.fr/en/data/products/, last access: 3 March 2020) are the satellite observational products used in this study. For ocean waves, the AVISO has two satellite altimetry products: along-track data and gridded data. The along-track data are used as the observational data source for the wave data in the simulation (sampled from “truth” in the “twin” experiments). The gridded data are used to validate the wave simulation and assimilation ($\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ resolution with 1 d time intervals). During wave simulation, the significant wave height (SWH) is used as a basic observational variable for data assimilation, which is provided from three ongoing satellites: Jason-2, Jason-3, and Satellite for Argos and ALtiKa (SARAL). Figure 1 shows a one-cycle ground orbit by taking Jason-2 and SARAL as examples. Jason-3 is the successor of Jason-2, and both satellites share the same orbit.

Figure 1 The ground projection tracks of satellite Jason-2 (blue line, with an inclination of 66^∘ N(S)) and SARAL (red line, with an inclination of 88^∘ N(S)) in one cycle over approximately 10 d and 35 d, respectively, in the (a) global and (b) East Asia domains (zoomed out of green box in a).

2.2 Different modeling strategies in WW3 and SWAN

Since the observations are only a sample of real-world information, the model bias (i.e., systematic difference between a numerical wave model and the real world) is not well defined against the real world. In this study, we use the systematic difference between the WW3 and SWAN models to simulate the model bias and study the influences on wave data assimilation (WDA).

First, let us identity the difference in physical and numerical aspects to comprehend the causes of “bias” between these two models. In general, WW3 addresses global scales, and SWAN is more applicable in shallow water. Although the two models have most of the same physical processes, such as the wind input and nonlinear wave–wave interactions, each can provide multiple parameterization schemes to choose. For example, the nonlinear wave–wave interactions in SWAN include the discrete interaction approximation (DIA) (Hasselmann et al., 1985) and the Webb–Resio–Tracy (Resio and Perrie, 1991; Van Vledder, 2006; Webb, 1978), while there are more choices in WW3, such as the Generalized Multiple DIA (Toman, 2004, 2013), the two-scale approximation and full Boltzmann integral (Perrie et al., 2013; Perrie and Resio, 2009; Resio et al., 2011; Resio and Perrie, 2008), as well as the nonlinear filter scheme (Tolman, 2011). In order to reduce more uncertainties, the same scheme is chosen for the same physical process. In numerical aspects, there exist different implementation strategies such as the differencing method, which also contributes to bias.

2.3 Two data assimilation systems using WW3 and SWAN

To explore the model bias influences on WDA, we develop two data assimilation systems based on WW3 and SWAN in this study.

Generally, based on the program structure of wave models, we insert the assimilation module between calculations of the two-dimensional wave spectrum and outputs of wave parameters so that at the assimilation time, we call on the assimilation module to update the spectrum and SWH. When building the data assimilation systems, we need to consider the different structures of parallelism method, data storage and information exchange in WW3 and SWAN models as noted in Sect. 2.2.

To clearly demonstrate the influences of model bias on WDA and minimize its adverse impact, the analysis scheme in both assimilation systems is optimal interpolation (OI), which also is low cost and easy to operate. We implement the OI analysis in three to four steps. The first step uses two Gaussian convolutions of the background and observed SWHs to compute the observational increment of SWH at the observational location. The second step projects the SWH observational increment onto the model grids centered at the observational location but within an impact radius using linear regression. The third step transforms the analyzed SWH to update the spectrum of model waves. The fourth step corrects wind forcing using the observational SWH.

Step 1: computing observational increment by convolution of two Gaussians

Starting from the idea of the ensemble adjustment Kalman filter (Anderson, 2001), an observational increment at the observational location k, $\mathrm{\Delta }{H}_k^{\text{O}}$ (H represents SWH), is computed by the convolution of two Gaussians of the model background and observation, which can usually be obtained from model ensemble members and observational samples. $\mathrm{\Delta }{H}_k^{\text{O}}$ is formulated as follows (Zhang et al., 2007):

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{H}_k^{\text{O}}={\displaystyle \frac{\frac{\mathrm{1}}{{\left({\mathit{\sigma }}^{\text{M}}\right)}^{\mathrm{2}}}{\stackrel{\mathrm{‾}}{H}}^{\text{M}}+\frac{\mathrm{1}}{{\left({\mathit{\sigma }}^o\right)}^{\mathrm{2}}}{\stackrel{\mathrm{‾}}{H}}^{\text{O}}}{\frac{\mathrm{1}}{{\left({\mathit{\sigma }}^{\text{M}}\right)}^{\mathrm{2}}}+\frac{\mathrm{1}}{{\left({\mathit{\sigma }}^{\text{O}}\right)}^{\mathrm{2}}}}}+{\displaystyle \frac{\mathrm{\Delta }{H}_k^{\text{M}}}{\sqrt{\mathrm{1}+{\left(\frac{{\mathit{\sigma }}^{\text{M}}}{{\mathit{\sigma }}^{\text{O}}}\right)}^{\mathrm{2}}}}}-H_k^{\text{M}}.\end{array}$$

Here, the first and second terms on the right-hand side of the equation adjust the ensemble mean and ensemble spread, respectively, and $\mathrm{\Delta }{H}_k^{\text{M}}$ represents the prior model spread. Superscripts “O” and “M” denote the observation and prior quantity estimated by the model, respectively. σ is the corresponding error standard deviation of SWH and varies in time and space and differs in every wave model. The overbar denotes the ensemble mean. In this simplified case, we specify σ^M=0.6 m, σ^O=0.25 m, assuming that σ is constant in all conditions similarly to previous studies (e.g., Qi and Fan, 2013), and use single-model and observational values as the ensemble mean.

Step 2: regressing the observational increment onto model grids

The second step projects the observational increment $\mathrm{\Delta }{H}_k^{\text{O}}$ onto the related model grids using background error covariance, which is a key step in the analysis. To simplify the problem and improve the computational efficiency, many studies use a flow-independent distance function to sample the background error covariance for computing the analysis increment at the model grid i, $\mathrm{\Delta }{H}_i^A$, as $\mathrm{\Delta }{H}_i^A=({\mathit{\sigma }}_i^{\text{s}}{)}^{\mathrm{2}}\exp \left(-\left(\frac{d_{i,k}}{L}\right)\right)\times \mathrm{\Delta }{H}_k^{\text{O}}$. Usually, such an expression is only a symmetrical approximation of the correlation function and cannot represent the spatial structure and propagation characteristics of waves. Here, with the assumption that wave covariance has the same structure as wave error covariance, we modify the covariance formula to increase its representativeness for wave structure by superimposing a statistical correlation coefficient obtained from the outputted SWH of model simulation onto the formula. After analysis, the equation becomes

$$\begin{array}{ll}\text{(3)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{H}_i^A={\displaystyle \frac{{\mathit{\sigma }}_i^{\text{s}}}{{\mathit{\sigma }}_k^{\text{s}}}}{r}_{i,k}^{\text{s}}\exp \left(-\left({\displaystyle \frac{d_{i,k}}{L}}\right)\right)\times \mathrm{\Delta }{H}_k^{\text{O}},\text{if}{d}_{i,k}\le R{\displaystyle }& {\displaystyle }\mathrm{\Delta }{H}_i^A=\mathrm{0},\text{if}{d}_{i,k}>R,\end{array}$$

where L is the characteristic length and d_i,k is the distance between the model grid i and observational point k. When d_i,k is larger than the impact radius R, there is no observational impact on this model point from observation k. All variables with superscript “s” represent the model statistics from free model control results. For example, $r_{i,k}^{\text{s}}$ is the SWH covariance between the model grid i and observation k, which is evaluated from the model data time series in corresponding experiments. To ensure the local characteristics of ocean waves, in this study, the characteristics length L and impact radius R (or the largest d_i,k) are the same, causing this incremental projection to reach to the e-folding scale. Referring to previous studies (e.g., Lionello et al., 1992; Qi and Fan, 2013), we tested different values of L and R as 300, 800 and 1000 km and found no essential improvement with larger L and R values. Trading-off with computational efficiency, we set L and R to 300 km throughout this study. As shown in Fig. 2, the new covariance represents more wave physics; i.e., the correlation has more asymmetrical and wave-dependent characteristics.

Figure 2 Spatial distributions of (a) background correlation coefficients by the empirical correlation model (blue), model data statistics (green) as well as their combination (red) of Eq. (3), (b) background adjustment increments of SWH projected from the observational increment with the empirical correction model (blue line) and their combination model (red filled), and (c) the difference in background adjustment increments of SWH from (b) by an analysis process given the single observation obtained at 18.90^∘ N, 114.09^∘ E, denoted by the asterisk (unit: m).

Step 3: transforming the SWH to wave spectrum

The assimilation SWH $H_i^A$ is a sum of the prior $H_i^{\text{M}}$ and the analysis increment from step 2 ($H_i^A=H_i^{\text{M}}+\mathrm{\Delta }{H}_i^A$). In the wave model, the form of ocean waves is a two-dimensional wave spectrum that is distributed over frequency and phase. Thus, transforming the assimilation SWH to a wave spectrum is necessary to update other wave parameters. Following the previous study (Qi and Fan, 2013), we assume that the change in wave spectrum is proportional to the energy change that is expressed by the square of SWH. Then, the analyzed spectrum $S_i^A(f,\mathit{\theta })$ can be written as follows:

$$\begin{array}{}\text{(4)}& S_i^A\left(f,\mathit{\theta }\right)={\left({\displaystyle \frac{H_i^A}{H_i^{\text{M}}}}\right)}^{\mathrm{2}}{S}_i^{\text{M}}\left(f\mathit{\theta }\right),\end{array}$$

where f is the wave frequency and θ is the phase direction.

Step 4: correcting wind forcing using SWH data

If the assimilation only adjusts the wave spectrum as described in Step 3, the updated spectral structure may be quickly overwritten by erroneous wind. In this step, we describe a simple scheme using the observed SWH data to correct the wind forcing. Starting from a first guess of wind (the ERA-Interim reanalysis, for instance), the analyzed wind $W_{i,j}^A$ at model grid (i,j) can be written as follows:

$$\begin{array}{}\text{(5)}& W_{i,j}^A=W_{i,j}^{\text{M}}+\mathrm{\Delta }{W}_{i,j},\end{array}$$

where W represents either the u or the v component of wind. ΔW_i,j is the corrected wind increment transformed from the updated SWH. While the details of the transformation scheme can be found in Lionello et al. (1992, 1995), we comment on certain aspects relevant to our study. Regardless of boundaries, in general, the energy of ocean waves is determined by the wind speed and duration, which can also be expressed by SWH. In that sense, a function equation can be built, in which the left-hand side is an expression of wind speed and duration, while the right-hand side is an expression of SWH, and they are balanced through wave energy. Then, the analyzed wind speed can be resolved under the assumption that the duration is the same in both the prior and the analyzed fields.

With respect to the configuration of wave model data assimilation, the model time step is 15 min and the assimilation interval is 1 h. At the assimilation time, we assimilate the along-track observations within a 1 h time window centered at the time. After 10 d, all the observations will cover the global area. The wind data from the reanalysis products (ERA-Interim and NCEP-CFSR in this case) are available every 6 h. To incorporate the wind correction into the wind forcing of the model, we distribute the wind correction to the adjacent two time levels of wind data. As the process is looped forward as the wave model state is updated, the wind forcing is adjusted through the SWH assimilation.

2.4 Experimental design

Throughout this study, we use the symbol MA${}_{O\left(s\right)}^{\text{WF}}$ as the name for the assimilation experiment. Here, “MA” stands for the “assimilation model” and the subscript “O(s)” (or superscript “WF”) represents the observing system (or wind forcing) in the assimilation. The wind forcing is either the ECMWF ERA-Interim (hereafter known as ERAI) or NCEP-CFSR wind (hereafter known as CFSR). Wind forcing can also be corrected by observations of SWH (under this circumstance, the superscript “WF” is replaced by “ASSW”). The observations used in the assimilation could be the model data but are projected on the along-track points of satellite(s) if being used for the twin experiments. Under this circumstance, “O” represents the “model that produces observations” and “(s)” represents the satellite tracks used (J2-Jason-2, J3-Jason-3 and SA-SARAL, for instance). Otherwise, in the real-data assimilation experiments, the subscript “O(s)” directly lists the satellites that measure the SWH. For example, a symbol named ${\text{WW3}}_{\text{SWAN(J2)}}^{\mathrm{ERAI}}$ means that the assimilation model is WW3 (here, “MA” = “WW3”) forced by ERA-Interim wind (“WF” = “ERAI”), and the model producing observation is SWAN with the Jason-2 satellite track (“O” = “SWAN” and “(s)” = “(J2)”).

2.4.1 Twin experiments

Twin experiments refer to a type of observing system simulation experiment (OSSE), in which a model simulation is used to define the “true” solution of a data assimilation problem and the other model simulation is used to start the assimilation. The “observations” are samples of the truth with some white noise to simulate the observational errors. When the truth and assimilation are conducted by different (or identical) models, the framework is a “biased” (or “perfect”) model twin experiment. Within a twin experiment framework, any aspect of assimilation skills can be measured as the degree to which the truth is recovered through the assimilation.

(a) Perfect twin experiment

In a perfect twin experiment, we assume that the assimilation model and the observation are unbiased; i.e., both the instrument measuring and numerical modeling processes are sampling the same stochastic dynamical system. Such sampling only has random sampling errors without any systematic difference (bias). We can build this perfect model framework by using the same model to produce the truth as the assimilation model but with different initial conditions and wind forcings.

The observations from the observational time window (1 h) centered at the assimilating time can be created by sampling the truth SWH with the tracks of Jason-2, Jason-3 and SARAL, which will cover the global area in 10 d. In this circumstance, if WW3 (or SWAN) is used as the assimilation model, the truth is produced by the same WW3 (or SWAN) model. In the assimilation, we may start the model with different initial conditions and/or wind forcings to examine the influences of initial errors and wind forcing errors on the wave assimilation. Such a perfect twin experiment can be named WW3${}_{\text{WW3(s)}}^{\text{WF}}$ or SWAN${}_{\text{SWAN(s)}}^{\text{WF}}$.

(b) Biased twin experiment

To study the impact of model errors on wave assimilation, we use two models to design a biased twin experiment. Again, due to the scattering nature of the observations, it is difficult to obtain a complete picture of the model bias against the real world. Given the difference between the WW3 and SWAN models described in Sect. 2.2, we use these two models and their assimilation systems here to simulate the model bias and examine its influences on the WDA. We use the ERA-Interim reanalysis wind to force the WW3 (or SWAN) to produce the truth and observations but use the SWAN (or WW3) assimilation system to assimilate the observations. The degree to which the truth produced by different model-based assimilation systems is recovered by assimilating the observations is an assessment of the model bias influences on the WDA. Such a biased twin experiment can be named WW3${}_{\text{SWAN(s)}}^{\text{WF}}$ or SWAN${}_{\text{WW3(s)}}^{\text{WF}}$.

Under the biased twin experiment framework, we also conduct experiments to examine the impacts of observing systems on wave assimilations by increasing the observational information based on multiple satellite tracks. For example, we can examine the assimilation results of WW3${}_{\text{SWAN(J2)}}^{\text{WF}}$, WW3${}_{\mathrm{SWAN}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3})}^{\text{WF}}$ and WW3${}_{\mathrm{SWAN}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{WF}}$ (or SWAN${}_{\text{WW3(J2)}}^{\text{WF}}$, SWAN${}_{\mathrm{WW}\mathrm{3}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3})}^{\text{WF}}$ and SWAN${}_{\mathrm{WW}\mathrm{3}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{WF}}$) to understand the impacts of observing systems on different model-based assimilations.

2.4.2 Real-data assimilation experiments

In this study, we also conduct real-data assimilation experiments using WW3 and SWAN assimilation systems with real track data from Jason-2, Jason-3 and SARAL. Through real-data assimilation experiments with different model-based assimilation systems, we can (1) increase our understanding of the influences of model errors on the WDA and (2) study the method to reduce the model error influences on the assimilation results. The real-data assimilation experiments can be directly named, i.e., WW3${}_{\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA}}^{\text{WF}}$ or SWAN${}_{\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA}}^{\text{WF}}$.

Table 1 List of perfect model twin experiments.

3 Error sources in wave models and WDA

3.1 Influences of initial and wind forcing errors

Usually, wave numerical simulation can be improved by three methods: (1) reducing the errors in the initial conditions, (2) enhancing the accuracy of the wind forcing, and (3) improving the representation of the wave model and its parameterization.

In this section, we use perfect model twin experiments (as described in Sect. 2.4.1) to exclude model errors and explore the impact of wind forcings and initial conditions on the wave simulations. To compare the performances of the WW3 and SWAN models, we conduct separate experiments with these two models. The truth and model control runs are two basic experiments of the perfect twin experiment framework. We use the ERA-Interim wind to drive WW3 (or SWAN) and generate a long time series of model states as the “truth,” which is called WW3^ERAI (or SWAN^ERAI) for the WW3 (or SWAN) perfect model twin experiment. The observations are created by interpolating the corresponding truth SWH onto the along-track points of satellite orbits. Then, we use the NCEP-CFSR wind to force WW3 (or SWAN), called the model control WW3^CFSR (or SWAN^CFSR), and the data assimilation is named WW3${}_{\text{WW3(s)}}^{\text{CFSR}}$ (or SWAN${}_{\text{SWAN(s)}}^{\text{CFSR}}$). Starting from an independent initial condition produced by the model control, we can conduct the assimilation with the ERA-Interim or NCEP-CFSR wind forcing. The error verification of the assimilation results against the truth simulation compared to the error of the model control is an evaluation of the initial error and/or wind forcing error influences on the WDA. All perfect model twin experiments are listed in Table 1.

First, we conduct two sets of model control experiments WW3^CFSR and SWAN^CFSRfor 80 d (from December 2017 to February 2018). To explore the effect of the initial conditions, we perform the model spin-up for a long time to adequately reach a steady state. Then using the 45th-day model states as the initial conditions, we conduct one more model simulation and data assimilation experiments for each model system as WW3^ERAI and WW3${}_{\text{WW3(J2)}}^{\text{ERAI}}$ as well as SWAN^ERAI and SWAN${}_{\text{SWAN(J2)}}^{\text{ERAI}}$. The root mean square errors (RMSEs) of these experiments against the truth are shown in Fig. 3 as the red (for WW3^CFSR+WW3^ERAI and SWAN^CFSR+SWAN^ERAI) and pink (for ${\text{WW3}}^{\mathrm{CFSR}}+{\text{WW3}}_{\text{WW3(J2)}}^{\mathrm{ERAI}}$ and ${\text{SWAN}}^{\mathrm{CFSR}}+{\text{SWAN}}_{\text{SWAN(J2)}}^{\mathrm{ERAI}}$) lines.

Figure 3 The time series of RMSEs of the (a) WW3 and (b) SWAN perfect model experiments in the model control run with the NCEP CFSRv2 wind (black, denoted as WW3^CFSR in panel (a) and SWAN^CFSR in (b), assimilating the “observed” data sampled by the tracks of Jason-2 (blue, denoted as WW3${}_{\text{WW3(J2)}}^{\text{CFSR}}$ and SWAN${}_{\text{SWAN(J2)}}^{\text{CFSR}}$), Jason-2 and 3 (green, denoted as WW3${}_{\mathrm{WW}\mathrm{3}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3})}^{\text{CFSR}}$ and SWAN${}_{\mathrm{SWAN}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3})}^{\text{CFSR}}$), as well as Jason-2 and 3 and SARAL (cyan, denoted as WW3${}_{\mathrm{WW}\mathrm{3}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{CFSR}}$ and SWAN${}_{\mathrm{SWAN}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{CFSR}}$) against the truth simulation forced by the ERA-Interim wind. The red and pink are forced by the NCEP CFSRv2 wind in the first 45 d, but the next 35 d are forced using the ERA-Interim wind (same as truth) without (denoted as WW3^CFSR+WW3^ERAI and SWAN^CFSR+SWAN^ERAI) or with (denoted as WW3^CFSR + ${\text{WW3}}_{\text{WW3(J2)}}^{\mathrm{ERAI}}$ and SWAN^CFSR + ${\text{SWAN}}_{\text{SWAN(J2)}}^{\mathrm{ERAI}}$) the assimilation of Jason-2 data. The number in parentheses for each color is the corresponding RMSE averaged globally over the verification time period (30 d after the 45 d model spin-up and 5 d assimilation spin-up). The observed data are produced by projecting the truth SWH onto the satellite orbit.

The SWH RMSE is approximately 0.34 m in the WW3 or SWAN model control with the NCEP-CFSR wind. Once the wind forcing is changed to the perfect wind (the ERA-Interim) on the 45th day, the RMSE quickly drops and is close to zero after approximately 10 d, and the SWAN model takes longer to accomplish this change than the WW3. If data assimilation is added, the RMSE reduces much faster than the model controls (roughly half of the timescale of the correct wind forcing). From the analyses above, we learned (1) in wave models, the wind forcing plays an important role and an incorrect wind forcing could be a significant error source of WDA, and (2) the WDA can rapidly reduce the initial error and improve the predictability of a wave model even when it is forced by an inaccurate wind forcing.

3.2 Impact of the observational system

In this section, using the same model states (at the 45th day) in the corresponding model control as in Sect. 3.1 as the initial conditions, we conduct two sets of assimilation experiments: WW3${}_{\text{WW3(J2)}}^{\text{CFSR}}$, WW3${}_{\mathrm{WW}\mathrm{3}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3})}^{\text{CFSR}}$, WW3${}_{\mathrm{WW}\mathrm{3}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{CFSR}}$ and SWAN${}_{\text{SWAN(J2)}}^{\text{CFSR}}$, SWAN${}_{\mathrm{SWAN}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3})}^{\text{CFSR}}$, SWAN${}_{\mathrm{SWAN}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{CFSR}}$. Through examining the assimilation quality with one satellite (Jason-2), two satellites (Jason-2+Jason-3) and three satellites (Jason-2+Jason-3+SARAL), we attempt to understand the impact of improving the observing system on the WDA, considering the NCEP-CFSR wind forcing errors against the ECMWF ERA-Interim based on a perfect assimilation model. The RMSEs of all the above assimilation experiments are plotted in Fig. 3 as the blue (assimilating Jason-2 only), green (assimilating Jason-2+Jason-3) and cyan (assimilating Jason-2+Jason-3+SARAL) lines.

From Fig. 3, we can see that in both models, the assimilation errors are reduced when more observational information is used. The corresponding RMSE reductions in these three experiments from the model control run are 24 %, 32 % and 38 % for WW3 and 26 %, 35 % and 38 % for SWAN, respectively. However, when more satellite observations are assimilated into the model, the magnitude of improvement becomes small (further reduced by 8 % in WW3 and 9 % in SWAN when Jason-3 is added as well as only a 6 % in WW3 and 3 % in SWAN further reduction when SARAL is further added). These results suggest that given wind forcing errors, increasing observational information can help to improve the model behavior, but the improvement is limited.

3.3 Adverse impact of model bias

As described in Sect. 2.2, the WW3 and SWAN models discretize the wave-action-governing equation with different physical processes, parameterization schemes and differencing schemes. These differences result in each wave model having its own distinguished characteristics. To study the adverse impact of the model bias on the wave assimilation, the biased twin experiments described in Sect. 2.4.1 are used in this section, where the truth model and the assimilation model are different between WW3 and SWAN. For example, the WW3${}_{\text{SWAN(J2)}}^{\text{ERAI}}$ (or SWAN${}_{\text{WW3(J2)}}^{\text{ERAI}}$) experiment uses WW3 (or SWAN) as the assimilation model to assimilate the Jason-2 track point observations, but the observed values are produced by SWAN (or WW3), and all models are forced by the ERA-Interim wind. All related experiments for the biased model framework are described in detail in Table 2.

Table 2 List of biased model twin experiments.

The RMSEs and correlation coefficients produced by the all biased model assimilation experiments are plotted in Fig. 4. The black line in each panel represents the result of the WW3 model control forced by the ERA-Interim wind (WW3^ERAI) against the truth simulation by the SWAN model with the same wind forcing (SWAN^ERAI) (panels a and b) (vice versa in panels c and d). Both the WW3^ERAI and SWAN^ERAI experiments are initialized from a cold start by the wind and integrated for 80 d, and the results of the last 40 d are shown in Fig. 4. It is clear that the WW3 and SWAN model simulations are quite different even though both simulations use identical forcings and start from identical initial conditions. The RMSEs of the two model simulations are both 0.58 m, which is much larger than the errors produced by a perfect model but with different wind (∼0.34 m, see black lines in Fig. 3).

Figure 4 The time series of RMSEs (a, c) and spatial correlation coefficients (b, d) of the WW3 (a, b) and SWAN (c, d) biased model experiments in the model control run forced by the ERA-Interim wind (black, denoted as WW3^ERAI and SWAN^ERAI), assimilations with “observed” data from one (pink, denoted as WW3${}_{\text{SWAN(J2)}}^{\text{ERAI}}$ and SWAN${}_{\text{WW3(J2)}}^{\text{ERAI}}$), two (red, denoted as WW3${}_{\mathrm{SWAN}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3})}^{\text{ERAI}}$ and SWAN${}_{\mathrm{WW}\mathrm{3}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3})}^{\text{ERAI}}$) and three (blue, denoted as WW3${}_{\mathrm{SWAN}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{ERAI}}$ and SWAN${}_{\mathrm{WW}\mathrm{3}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{ERAI}}$) satellites, as well as the assimilation with corrected wind (dotted green, denoted as WW3${}_{\mathrm{SWAN}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{ASSW}}$ and SWAN${}_{\mathrm{WW}\mathrm{3}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{ASSW}}$) against the truth (same as in Fig. 3 but for SWAN and WW3 with the ERA-Interim wind). The numbers in the parentheses correspond to the globally averaged RMSE (in a and c) and spatial correlation coefficient (in b and d) over the last 30 d during the assimilation period. The observed data are produced by projecting the truth SWH onto the satellite orbit.

Compared with the model controls WW3^ERAI and SWAN^ERAI, the assimilation experiments WW3${}_{\text{SWAN(J2)}}^{\text{ERAI}}$ and SWAN${}_{\text{WW3(J2)}}^{\text{ERAI}}$ (pink lines in Fig. 4) can significantly reduce the SWH simulation error (by 24 % and 22 %, respectively) and enhance the correlations (by 3 % and 4 %, respectively) with the truth (SWAN^ERAI and WW3^ERAI, respectively). When the observations of Jason-3 and SARAL are added to the assimilation (i.e., WW3${}_{\mathrm{SWAN}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3})}^{\text{ERAI}}$ and WW3${}_{\mathrm{SWAN}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{ERAI}}$ as well as SWAN${}_{\mathrm{WW}\mathrm{3}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3})}^{\text{ERAI}}$ and SWAN${}_{\mathrm{WW}\mathrm{3}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{ERAI}}$) (see the red and blue lines, respectively), the model SWH error (or correlation) is further reduced (or enhanced), but the amplitude of reduction (or enhancement) gradually diminishes (10 % and 5 % for further error reduction and 1 % and 0.8 % for further correlation enhancement in the WW3 assimilation; 10 % and 7 % for further error reduction and 1.7 % and 0.7 % for further correlation enhancement in SWAN assimilation from the additions of Jason-3 and SARAL, respectively).

The results of two other sets of assimilation experiments called WW3${}_{\mathrm{SWAN}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{ASSW}}$ and SWAN${}_{\mathrm{WW}\mathrm{3}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{ASSW}}$ are also plotted by dotted green lines in Fig. 4. The superscript “ASSW” stands for the assimilation-corrected wind, meaning that the wind forcing of the assimilation model is also “corrected” by the “observed” SWH data, as described in Step 4 of Sect. 2.3. We found that in both assimilation systems, using the observed SWH data to “correct” the wind can compensate for the model errors to some degree and further reduce the assimilation errors, but the improvement is very limited. The weak improvement could be attributed to the simple wind correction method with total wave height. In the future, a more powerful correction method with wind sea wave height may have better performance.

These assimilation results clearly show that even though the wind forcing is perfect, once a biased assimilation model is used, the wave simulation has large errors. Although WDA can greatly reduce the simulation error by assimilating the observational information into the model, due to the existence of the model bias, the error remains at some significant level and cannot be eliminated entirely even by increasing the observational constraints through an improvement in the observational system and constraint of the wind forcing.

3.4 Comparison of the influence of wind forcing with model bias

Comparing the time series of SWH RMSE caused by two important error sources, model bias (0.58 m shown with black lines in panels a and c in Fig. 4) plays a stronger role than wind forcing (0.34 m shown with black lines in Fig. 3). When three satellite observations (Jason-2, Jason-3 and SARAL) are assimilated into twin experiments, the RMSE of SWH reduces by 38 % in perfect model twin experiments (cyan lines in Fig. 3) and 40 % in biased model twin experiments (blue lines of panels a and c in Fig. 4). It is obvious that the error caused by model bias is bigger than by wind forcing and their improvements are almost similar after data assimilation.

Figure 5 Distributions of SWH RMSEs caused by wind forcing (a, b) from perfect twin experiment and model bias (c, d) from biased twin experiment in the model control run (a, c) and assimilations with “observed” data with three satellites (b, d) averaged globally over the last 30 d during the assimilation period.

A spatial pattern of SWH RMSE is also displayed in Fig. 5. (Due to the similar performance inside twin experiments, here we only show the results using WW3 model as simulation/assimilation model.) In the model control run, it makes sense that the error caused by wind forcing (panel a (WW3^CFSR); truth is WW3^ERAI) has distributed in the areas where the wind is strong, such as the area of the Antarctic Circumpolar Current (ACC) and the north Pacific and Atlantic Ocean, while the error caused by model bias (panel c (WW3^ERAI); truth is SWAN^ERAI) is distributed almost globally. After assimilating with three satellite observations, the error spatial distribution of both has improved a lot (panel b (WW3${}_{\mathrm{WW}\mathrm{3}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{CFSR}}$) and panel d (WW3${}_{\mathrm{SWAN}(\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA})}^{\text{ERAI}}$)). Referring to wind distribution, we can divide the global areas roughly into three parts: the northern westerly zone (65^∘ N–30^∘ N), the trade-wind zone (30^∘ N–30^∘ S) and the circumpolar westerly zone (30^∘ S–65^∘ S). Therefore the error is caused by wind forcing (or model bias); the decreasing percentages of SWH RMSE in these three areas are 30 % (or 27 %), 45 % (or 50 %) and 46 % (or 48 %), respectively. We can find that the improvement of both error sources has a similar performance in three these areas: weak in the northern westerly zone and almost the same strength in the trade-wind zone and circumpolar westerly zone. The reason why there is a lower improvement in the north Pacific and Atlantic Ocean should be explored in the future. To sum up, the error caused by model bias is larger than wind forcing in the global area generally, especially in the equatorial ocean. However, in the north of Atlantic Ocean, wind forcing has a stronger impact. After data assimilation, both have improved greatly and have a similar spatial pattern (i.e., the bigger error is still distributed at high latitudes). However, their error gap still exists.

It is worth mentioning that there is a similar performance about the effect of the observation system on improving both error sources, whether in time series or spatial distribution. The more observation information is absorbed into assimilation systems, the better the error improvement in both twin experiments. However, due to the existence of model bias, this improvement has a limitation and stays at a certain level. If more powerful observation is absorbed (such as wave direction, wave period and two-dimensional wave spectra), the limitation maybe stay at a lower level. Next, with the results of real-data assimilation where both the model and wind forcing have errors, we will analyze and discuss how to mitigate the model bias influences on the WDA.

Figure 6 Distributions of SWH mean errors (against the merged grid altimeter data) of the WW3 (a–c) and SWAN (d–f) model simulations (a, d) and assimilations with Jason-2 (b, e), as well as all Jason-2, Jason-3 and SARAL (c, f) data forced by the ERA-Interim wind. The statistics are averaged over the last 30 d of a 70 d total assimilation period (unit: m).

4 Mitigation of model bias influences on wave assimilation

4.1 Bias characteristics of WW3 and SWAN data assimilations

From the above analyses of twin experiment results, we learned that the model bias has a strong adverse impact on WDA. To explore the method of mitigating the model bias influence on the WDA, we conduct the real-data assimilation experiments (same time range as the twin experiments) described in Sect. 2.4.2 using the WW3 and SWAN assimilation systems to assimilate the track data of Jason-2 and Jason-2+Jason-3+SARAL. To ensure the performance of the biased model WDA, a longer assimilation (more than 2 months) is conducted (a total of 70 d). The spatial distributions of the SWH errors (obtained from the difference against the merged gridded AVISO observations over the last 30 d out of 80 d) are shown in Fig. 6. Panels a, b and c (d, e, f) are for the WW3 (SWAN) simulation and assimilations: WW3^ERAI, WW3${}_{\text{J2}}^{\text{ERAI}}$ and WW3${}_{\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA}}^{\text{ERAI}}$ (or SWAN^ERAI, SWAN${}_{\text{J2}}^{\text{ERAI}}$ and SWAN${}_{\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA}}^{\text{ERAI}}$).

Comparing panel a with panel d in Fig. 6 reveals that a large difference exists in the simulations of the two models. First, the SWAN simulation errors are generally larger than the WW3 simulation errors. Second, the global error distributions are quite different: while the WW3 simulation errors appear negative (or positive) over most of the 30^∘ S north (or south) area, the SWAN simulation errors appear positive in most of the tropical oceans and negative in the middle latitudes. Both simulations show large errors in the southern ocean coastal area, but over the Antarctic Circumpolar Current area, the WW3 error is positive and the SWAN error is negative. It is interesting that under the same wind conditions, two wave models have converse performances. In a future study, it is urgently necessary to find a detailed approach to two wave models in this area.

The above systematic differences between the two model simulations have significant influences on the results of the WDA. In general, the distribution of assimilation errors shares the same patterns as the model simulation errors but with a much smaller magnitude. The net result is that both the WW3 negative (or positive) error magnitude over the 30^∘ S north (or south) area and the SWAN error magnitude as (+)(−)(+)(−) from south to north are dramatically reduced by the Jason-2 data assimilation (comparing Fig. 6b and e with Fig. 6a and d), and on this basis, incorporating more observations from Jason-3 and SARAL into the assimilation process, both model error magnitudes are further reduced to some degree (comparing Fig. 6c and f with Fig. 6b and e). From the corresponding RMSE distributions (not shown here), we learned that the large RMSEs mainly appear in places where the model bias (time mean error) is large. This finding means that the model bias has a largely adverse impact on the WDA.

Figure 7 Time series of RMSEs (a, c) and spatial correlation coefficients (b, d) of WW3 (a, b) and SWAN (c, d) produced by the model control run (black, denoted as WW3^ERAI and SWAN^ERAI) and the assimilation using the data from one (green, denoted as WW3${}_{\text{J2}}^{\text{ERAI}}$ and SWAN${}_{\text{J2}}^{\text{ERAI}}$) and three satellites (blue, denoted as WW3${}_{\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA}}^{\text{ERAI}}$ and SWAN${}_{\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA}}^{\text{ERAI}}$) with corrected wind (red, denoted as WW3${}_{\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA}}^{\text{ASSW}}$ and SWAN${}_{\mathrm{J}\mathrm{2}+\mathrm{J}\mathrm{3}+\mathrm{SA}}^{\text{ASSW}}$).

Figure 7 displays the time series of the RMSEs and spatial correlation coefficients with the global statistics in space. The RMSE (or correlation coefficient) of the SWAN model simulation is larger (or smaller) than the WW3 model simulation (0.66 m RMSE and 0.806 correlation for SWAN versus 0.61 m RMSE and 0.876 correlation for WW3). In the WW3 and SWAN assimilations with the Jason-2 data, the RMSEs are reduced by 8 % and 11 %, respectively, and the time mean correlations are enhanced by roughly 1 % and 5 %, respectively. If the data of all three satellites – Jason-2, Jason-3 and SARAL – are assimilated, the RMSEs are reduced by 11 % and 17 % in the WW3 and SWAN assimilations, respectively, and the correlations are enhanced by approximately 2 % and 8 %, respectively. In these real-data assimilation cases, for each model assimilation, both the model and wind forcing have errors. Under this circumstance, the assimilation where the SWH observations are used to adjust the model spectrum; the model wind forcing is also corrected and can further reduce the assimilation errors (red lines in Fig. 7). The red lines in Fig. 7 represent the best result of the assimilation given the WW3 and SWAN model biases, which makes full use of the observations from all three satellites to adjust both the model spectrum and wind forcing. Next, we will discuss how to use the results of two assimilation systems to mitigate the wave analysis error.

Figure 8 Time series of (a) RMSEs and (b) spatial correlation coefficients produced by two bias correction schemes (cyan and pink) through a combination of WW3 and SWAN assimilations with the data from three satellites (Jason-2, Jason-3 and SARAL) and wind correction starting from the ERA-Interim wind. The results of the individual assimilation systems are plotted as dotted and dashed red lines (taken from Fig. 7) for reference.

4.2 Mitigation of WDA errors

The mitigation of model bias is a complex issue in which improving the model is a final but long-lasting solution. From Fig. 6, we learned that the WW3 and SWAN assimilation errors have some common (or opposite) characteristics in some locations. For example, while the SWH over the southern ocean coastal area always appears to be overestimated because of the lack of adequate observations to improve in both assimilation systems, the WW3 and SWAN assimilations appear to be the opposite in the Antarctic Circumpolar Current area and the tropical oceans. The WW3 (or SWAN) assimilation errors in the Antarctic Circumpolar Current area appear positive (or negative), while the WW3 (or SWAN) assimilation errors in the tropical oceans appear negative (or positive). A question arises: as the first step of mitigating model bias influences on the WDA, can we use a pair of assimilation systems to explore a statistical approach to reduce the wave assimilation errors?

Given the opposite behaviors of two assimilation systems existing in certain places, the simplest bias correction could be conducted by a simple average. This assumes that each wave model assimilation system has its own characteristics of systemic error due to deficit physics, and these systemic errors (or model biases) from different wave models follow a Gaussian distribution with a trivial expectation. The corresponding results are shown in Fig. 8. Compared with the performance of each individual assimilation system (dashed red lines for WW3 and dotted red lines for SWAN), the results of this bias correction (cyan lines) show that the RMSE is reduced (Fig. 8a) but the spatial correlation is not greatly improved (Fig. 8b). It is reasonable that based on the opposite errors deviating from the real world in two assimilation systems, this correction method employing the mathematical average can reduce the RMSE to some extent, but it may not have a significant contribution to improving the correlation coefficient if either the sampling size of model bias is too small (only two in this case) or the bias has an asymmetric distribution.

Figure 9 Distributions of corrected SWH RMSEs of the WW3 (a–d) and SWAN (e–h) model simulations (a, e); assimilations with Jason-2 (b, f) and with Jason-2, Jason-3 and SARAL (c, g) forced by the ERA-Interim wind as well as by a corrected wind (d, h). All the results is corrected by the second bias correction method (named Bias correction2 in Fig. 8) (unit: m).

Considering the potential asymmetry of the Gaussian distribution of different model biases and small sampling size in practice, as the first step, we calculate the spatial distribution of model bias in every assimilation system (the time mean of the difference between observations and assimilation results) (Fig. 6) and extract it at each time step in the output and then calculate the expectation (average) of all assimilation systems as the results after bias correction. The corresponding results are shown with pink lines in Fig. 8. Both the RMSE and correlation are improved greatly. We also show the spatial distribution of SWH RMSEs after bias correction with the second method. From Fig. 9, we can easily find that after bias correction, all the experiments have similar error performances, the smaller at low latitudes and the larger at high latitudes. As the satellite observation increases (one in panels b and f and three in panels c and g), the error decreases gradually compared to the model control run (panels a, e). Based on the best assimilation results currently (panels c, g), the improvement with a corrected wind is more efficient in WW3 (panel d), especially at high latitudes which is the higher error area, but it is hardly to found in SWAN (panel h). If we focus on the change without (panel c) and with (panel d) the corrected wind in the WW3 model, the RMSE improvement of wind sea and swell is 27 % and 1 % averaged globally. If we divide the global area into the three parts mentioned in Sect. 3.4, the wind sea (or swell) improvement of the regional RMSE mean is 4 % (or 3 %), 3 % (or 1 %) and 6 % (or 1 %) in the northern westerly zone (65–30^∘ N), the trade-wind zone (30^∘ N–30^∘ S) and the circumpolar westerly zone (30–65^∘ S). It is obvious that the improvement of wind sea is better than swell and wind correction has a stronger impact on wind sea wave assimilation at high latitudes with strong wind. Clearly, this bias correction with physical consideration is more effective to improve the quality of WDA, but it uses observational information one more time, while the first method of bias correction processes assimilation results directly without further uses of observational information.

To verify the feasibility and applicability of the bias correction method above, three well-known wave models (WW3, SWAN and WAM) with the same data assimilation method are used to conduct longer assimilation and bias correction experiments. The calculation period lasts for 14 months (from November 2016 to December 2017) with sufficient spin-up process to reach a steady assimilation state (the first month for model spin-up and the second month for assimilation spin-up). The results of the last 12 months (for 2017) are analyzed and presented in Figs. 10 and 11 for the spatial distributions and time series of RMSEs and spatial correlation coefficients, respectively. From Fig. 10, we can see that both the RMSE and the correlation coefficient (panels d and h, respectively) have been improved by the bias correction that combines the advantages of every WDA system (panels a and e for WW3, panels b and f for SWAN, panels c and g for WAM). In Fig. 11, the bias correction of model control runs (the second bias correction method is conducted based on the result of model control run from three wave models) shows improvement but is worse than the data assimilation before bias correction (compare green with pink). Compared with the model control (blue), the assimilation results with bias correction (red) can reduce the error by 25 % and significantly enhance the correlation coefficient (from 0.88 to 0.923). This result confirms that this bias correction based on multiple assimilation systems can effectively enhance the WDA quality.

Figure 10 Distributions of SWH RMSEs (a–d) and spatial correlation coefficients (e–h) (against the merged grid altimeter data) of the WW3 (a, e), SWAN (b, f) and WAM (c, g) assimilations and the bias correction with the second method (d, h). The statistics are averaged over the entire 1-year assimilation period.

Figure 11 The time series of (a) RMSEs and (b) spatial correlation coefficients produced by the model control run (blue), data assimilation (pink) and their corresponding bias corrections (green and red) combining three wave model assimilation results (WAM, WW3 and SWAN) globally averaged over 1 year.

5 Summary and discussion

Ocean waves cause the sea surface roughness to impact the boundary conditions of the atmosphere and the wind stress of the ocean surface. Wave processes, such as wave breaking and wave-induced bottom stress, have significant effects on ocean mixing. Thus, ocean waves are important physical processes for understanding ocean mixing and air–sea interactions in coupled Earth systems. More accurately predicting ocean waves is of great societal significance. However, multiple error sources exist in wave simulations and predictions, including modeling errors, wind forcing errors and initial-condition errors.

To sort out the source of the errors of wave data assimilation (WDA), a pair of independent WDA systems is first developed using two wave models: Wave Watch III (WW3) and Simulating WAves Nearshore (SWAN). The perfect and biased model twin experiment frameworks are designed to clearly identify each error source and examine its influences on WDA. The results show that model bias is a significant error source that has a largely adverse impact. Then, two WDA systems are used to design bias correction approaches to mitigate the influences of model bias and improve the assimilation quality. Finally, long-term WDA experiments added by the third WDA system with the WAM model (WAve Modeling) (WW3, SWAN and WAM) are conducted to validate the bias correction method. Three findings are established. (1) When the model is perfect, the initial-condition error decays within 10 d, but the WDA can shorten the timescale by half. (2) When the model is biased, despite a perfect wind forcing, the wave simulation has large errors and the WDA can only reduce the error to a limited extent. (3) With the results from two assimilation systems, a statistical approach of bias correction significantly improves the quality of final wave analysis by combining the merits from individual assimilation systems.

Model bias is an obstacle to improving WDA and wave predictions. Using multiple assimilation systems to study the influences of model bias on WDA is an effective approach. As the first step, however, we have used a simple assimilation scheme and simple bias correction method. In follow-up studies, we shall consider powerful observation information (such as two-dimensional wave spectra), advanced assimilation schemes (such as ensemble Kalman Filter) and more comprehensive correction methods to help improve modeling. For example, the “online” bias correction (the “offline” bias correction is used in this paper) during the assimilation process (e.g., Dee, 2005) will be considered to improve the assimilation results and correct the instantons initial condition within individual assimilation systems; after that, a robust observation (such as a buoy) is needed to validate the quality of bias correction. In addition, improving the model is an important, inevitable and long-lasting task. In this study, under the same wind conditions, we find that three models show common bias characteristics in the Antarctic Circumpolar Current (ACC) area. If there is a similar performance forcing by other wind, this may suggest that present wave modeling may have deficits in energy spectrum expression for high-wind-speed areas. In the future, we will further examine the sensitivities of physical processes on high wind speed to mitigate such a common modeling bias. All in all, a robust bias correction method with a lower model bias and higher representation of wave physical characteristics may further improve wave analysis quality. Once a long time series of high-quality wave analyses is available, it is expected that we can improve our understanding of ocean mixing. The physical process of wave-induced mixing is linked with the structure of the ocean mixing layer (Qiao et al., 2010). This process can be expressed as a function of wave number, frequency and wave spectrum and so on, provided by wave analysis. With the framework of multiple WDA systems developed in this study, improved wave predictions can be effectively pursued. How can we further enhance the predictability of ocean waves? The first important step is to understand the physical process of ocean waves better based on a more accurate evolution of wave state from this framework. Answering these questions could be very important and interesting research topics for future studies.