Deriving hourly diagnostic surface velocity fields considering inertia and an application in the Yellow Sea

Surface currents play an important role in the transport of floating materials in the Yellow Sea, a region strongly influenced by tidal forcing and seasonal wind variability driven by the East Asian monsoon. While diagnostic models have been widely used to estimate surface currents, due to their steady-state assumption, high frequency variations such as tides and inertial oscillations cannot be resolved. To address this limitation, a time-dependent diagnostic model incorporating inertial terms into the governing equations is proposed. The performance of the proposed method is evaluated using buoy and drifter observations from 2015 to 2023. The time-dependent model in this study captures not only low frequency components (geostrophic and Ekman currents) but also high frequency variability (inertial oscillations and tides). Compared to the traditional model assuming steady-state, it shows significant improvement, achieving a correlation of 0.76 and Root-Mean-Square Error of 0.18 m s⁻¹ (compared to −0.08 and 0.43 m s⁻¹ for the steady model, that caused by wrong governing equation ignoring inertia to describe tides) because of successful consideration of high frequency variability. The decay rate of inertial oscillations is analytically derived, providing insight into the time scale for past signals in surface currents to dissipate. We expect that this study offers a practical framework for surface current estimation considering both high and low frequency signals and can be applied for quick assessments of material transport in other coastal oceans.

1 Introduction

Surface currents in oceanic environments are essential to understanding the transport of momentum and surface-floating materials. In the marginal sea such as the Yellow Sea (Fig. 1), there are many issues caused by floating materials such as jellyfish blooms (Wei et al., 2015; Choi et al., 2018), harmful algal blooms (Onitsuka et al., 2010; Kim et al., 2016), green algae blooms (Son et al., 2015; Kwon et al., 2019), and Sargassum (Shin et al., 2022; Choi et al., 2023). Thus, there is an increasing demand for estimating surface velocity fields. Numerical hydrodynamic models are widely utilized to generate velocity fields and simulate the transport of surface-floating materials (Onitsuka et al., 2010; Son et al., 2015; Wei et al., 2015; Kwon et al., 2019), with advancements in data assimilation significantly enhancing their accuracy (Storto et al., 2019). Along with numerical hydrodynamic models, the diagnostic approach, which estimates surface velocity fields from sea surface height and wind datasets, has become a popular method due to its efficiency and accuracy (Bonjean and Lagerloef, 2002; Rio et al., 2014; Dohan, 2017; Choi et al., 2023).

Figure 1 Geographic locations of the study area. Red triangle indicates location of the Southern Jeju Island buoy from the Korea Hydrographic and Oceanographic Agency (KHOA). Yellow solid lines indicate the trajectories of drifters from the National Institute of Fisheries Science (NIFS) and solid circles indicate the initial launching positions of the drifters. The orange arrow field shows the geostrophic current fields averaged over the study period (2013–2023).

The velocity fields from the previous diagnostic approach assumed that surface currents are primarily driven by the Ekman and geostrophic balance, so steady-state assumption was used. For instance, the Ocean Surface Current Analyses Real-time (OSCAR) velocity fields, provided by the National Aeronautics and Space Administration (NASA), were based on the diagnostic approach considering geostrophic and wind-driven Ekman current components (Bonjean and Lagerloef, 2002). Similarly, the Global total surface Current (GlobCurrent) dataset was produced by combining satellite-based geostrophic current and a semi-empirical model for the Ekman current component, aiming to accurately represent surface current variability (Rio et al., 2014). Recently, diagnostic velocity fields based on a simple analytical solution have been utilized and verified in the Yellow Sea to simulate algae transports and it was showed that velocity fields from the diagnostic approach can have comparable (or even locally better) accuracy to those from reanalysis datasets using numerical models and data assimilation (Choi et al., 2023). The seasonal Ekman currents driven by the East Asian Monsoon and the altimetry-based geostrophic currents, directed northeastward without a notable mean current pattern (orange arrows in Fig. 1), successfully simulated the transport of surface floating substances in the Yellow Sea (Choi et al., 2018, 2023).

Consequently, these diagnostic velocity fields developed in previous studies assumed steady states and considered only two low frequency signals (Ekman and geostrophic current components), so have limitations in capturing high frequency signals. In the Yellow Sea, there are significant high frequency variability of surface current driven by the near-inertial oscillations and tides (Teague et al., 1998; Hu et al., 2023). To consider these high frequency signals of surface current, new approach resolving the inertia is required. Time-dependent Ekman theory was advanced by Elipot and Gille (2009), Wenegrat and McPhaden (2016), and Lilly and Elipot (2021). These studies introduced and discussed analytical solutions incorporating the inertial terms, extending the classical steady Ekman theory (Ekman, 1905; Welander, 1957). Elipot and Gille (2009), in particular, assessed transfer functions of the analytical solutions based on spectral analysis of the observations.

Choi et al. (2023) showed that a simple analytical solution from the steady Ekman theory (Ekman, 1905; Welander, 1957) successfully reconstructs low frequency surface current in the Yellow Sea, that motivates this study to use a simple analytical solution from the time-dependent Ekman theory. This study builds techniques using the time-dependent Ekman theory to generate velocity fields considering both high and low frequency signals (Sect. 2). We expand the solution to consider time-dependency in the pressure gradient terms to resolve tidal velocity components. Based on the analytical solution and transfer function approach, diagnostic surface current velocity fields considering inertial term, which resolve high frequency variability, were derived and verified by buoy and drifter observations in the Yellow Sea. The velocity fields based on the time-dependent Ekman theory are compared with those based on the steady Ekman theory used by Choi et al. (2023).

2 Theory

2.1 Governing equation and the analytical solution for time-dependent Ekman theory

The governing equation and vertical boundary conditions for time-dependent Ekman theory are given by

$$\begin{array}{}\text{(1a)}& {\displaystyle \frac{\partial \mathit{u}}{\partial t}}+fi\mathit{u}=-g{\displaystyle \frac{\partial \mathit{\eta }}{\partial \mathit{n}}}+A_z{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{u}}{\partial z^{\mathrm{2}}}},\text{(1b)}& {\displaystyle }{\left.A_z{\displaystyle \frac{\partial \mathit{u}}{\partial z}}\right|}_{z=\mathrm{0}}={\displaystyle \frac{{\mathit{\tau }}^{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{0}}}},\text{(1c)}& {\displaystyle }{\left.\mathit{u}\right|}_{z=-h}=\mathrm{0}\end{array}$$

where $\mathit{u}=u+iv$ represents horizontal velocities in the complex coordinates, with i as the imaginary unit ($i=\sqrt{-\mathrm{1}}$), t is time, f is the Coriolis parameter, g is gravitational acceleration, A_z is the vertical eddy viscosity assumed to be constant, and $\partial \mathit{\eta }/\partial \mathit{n}=\partial \mathit{\eta }/\partial x+i\partial \mathit{\eta }/\partial y$ is barotropic pressure gradient in the x and y direction and presumed to be a given complex time series. In the surface and bottom boundary conditions (1b) and (1c), ${\mathit{\tau }}^{\mathrm{s}}={\mathit{\tau }}_x^{\mathrm{s}}+i{\mathit{\tau }}_y^{\mathrm{s}}$ is the surface wind stress, ρ₀ is a constant density, h is water depth. The governing Eq. (1) is identical to that of Choi et al. (2023) except for considering the inertial term. It must be noted that assuming constant eddy viscosity is a strong idealization, which can be relaxed through more complex analytical solution considering spatiotemporal variability of eddy viscosity (Elipot and Gille, 2009; Wenegrat and McPhaden, 2016; Lilly and Elipot, 2021). Nonetheless, this assumption is retained here to focus on the impact of incorporating inertial terms and to ensure consistency with the formulation in Choi et al. (2023). Using simplicity of Eq. (1), the equation will be used to inversely estimate monthly representative A_z from the velocity observations (Sect. 3.2.2). By using constant eddy viscosity consistent with Choi et al. (2023), where steady Ekman theory was used, impact of considering inertial terms can be solely assessed by the comparison between steady and time-dependent (non-steady) Ekman theory (discussed in Sect. 4).

The boundary value problem (1) is given by a linear nonhomogeneous partial differential equation. It includes an interior forcing term given by the barotropic pressure gradient $\partial \mathit{\eta }/\partial \mathit{n}$ and surface wind stress forcing τ^s via surface boundary condition (1b). This partial differential equation can be solved using typical approaches using the superposition principle. Based on this principle, the total velocity can be decomposed as $\mathit{u}={\mathit{u}}_{\mathrm{e}}+{\mathit{u}}_{\mathrm{g}}$, where u_e represents the wind-driven (Ekman) component, and u_g denotes pressure-driven (geostrophic) component. The governing Eq. (1) can be separated into equations for each component:

$$\begin{array}{}\text{(2a)}& {\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{e}}}{\partial t}}+fi{\mathit{u}}_{\mathrm{e}}=A_z{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{u}}_{\mathrm{e}}}{\partial z^{\mathrm{2}}}}\text{(2b)}& {\displaystyle }{\left.A_z{\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{e}}}{\partial z}}\right|}_{z=\mathrm{0}}={\displaystyle \frac{{\mathit{\tau }}^{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{0}}}}\text{(2c)}& {\displaystyle }{\left.{\mathit{u}}_{\mathrm{e}}\right|}_{z=-h}=\mathrm{0}\end{array}$$

$$\begin{array}{}\text{(3a)}& {\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{g}}}{\partial t}}+fi{\mathit{u}}_{\mathrm{g}}=-g{\displaystyle \frac{\partial \mathit{\eta }}{\partial \mathit{n}}}+A_z{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{u}}_{\mathrm{g}}}{\partial z^{\mathrm{2}}}}\text{(3b)}& {\displaystyle }{\left.A_z{\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{g}}}{\partial z}}\right|}_{z=\mathrm{0}}=\mathrm{0}\text{(3c)}& {\displaystyle }{\left.{\mathit{u}}_{\mathrm{g}}\right|}_{z=-h}=\mathrm{0}.\end{array}$$

Note that Eq. (2) was solved and analyzed by Elipot and Gille (2009). For the arbitrary wind stress that can be expressed as Fourier series ${\mathit{\tau }}^{\mathrm{s}}\left(t\right)=\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}{\widehat{\mathit{\tau }}}_n{e}^{-i{w}_{n}t}$, the solution at z=0 is given by

$$\begin{array}{}\text{(4)}& {\mathit{u}}_{\mathrm{e}}=\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}\underset{H_{{\mathrm{e}}_n}}{\underbrace{{\displaystyle \frac{\mathrm{2}}{{\mathit{\rho }}_{\mathrm{0}}{D}_{{\mathrm{e}}_n}(f-w_n)j}}\tanh \left({\displaystyle \frac{jh}{D_{{\mathrm{e}}_n}}}\right)}}{\widehat{\mathit{\tau }}}_n{e}^{-i{w}_{n}t}\end{array}$$

where $j=\mathrm{1}+i$ and $D_{{\mathrm{e}}_n}=\sqrt{\mathrm{2}{A}_z/\left(f-w_n\right)}$ indicating the Ekman depth modified by the presence of the inertial term. $H_{{\mathrm{e}}_n}$ represents the transfer function discussed by Elipot and Gille (2009). Similarly, solution to Eq. (3) can be obtained by the typical approach to solve the linear nonhomogeneous differential equation (e.g., method of undetermined coefficients). With $\partial \mathit{\eta }/\partial \mathit{n}=\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}{\widehat{\mathit{\eta }}}_n^{\prime }{e}^{-i{w}_{n}t}$, the solution at z=0 is given by

$$\begin{array}{}\text{(5)}& {\mathit{u}}_{\mathrm{g}}=\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}\underset{H_{{\mathrm{g}}_n}}{\underbrace{{\displaystyle \frac{ig}{f-w_n}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\cosh (jh/D_{{\mathrm{e}}_n})}}\right)}}{\widehat{\mathit{\eta }}}_n^{\prime }{e}^{-i{w}_{n}t}\end{array}$$

where $H_{{\mathrm{g}}_n}$ represents the transfer function for the velocity component. The first and second terms in Eq. (5) represent the pressure gradient-driven (geostrophic) current and its bottom frictional response for the current. The detailed method calculating the velocity components using Eqs. (4) and (5) will be illustrated in Sect. 3.2.

2.2 Dynamics of time-dependent Ekman theory and limitation of the classic steady Ekman theory

The solution to Eq. (1), given by the summation of Eqs. (4) and (5), aligns with the steady generalized Ekman solutions (Ekman, 1905; Welander, 1957; revisited by Choi et al., 2023), extending them to include time-dependency of the momentum (inertia term). The steady state solution is given by

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathit{u}}^{\ast }& =\underset{{\mathit{u}}_{\mathrm{e}}^{*}}{\underbrace{{\displaystyle \frac{\mathrm{2}\mathit{\tau }}{{\mathit{\rho }}_{\mathrm{0}}f{D}_{\mathrm{e}}^{*}j}}\tanh \left({\displaystyle \frac{jh}{D_{\mathrm{e}}^{*}}}\right)}}\\ & +\underset{{\mathit{u}}_{\mathrm{g}}^{*}}{\underbrace{i{\displaystyle \frac{g}{f}}{\displaystyle \frac{\partial n}{\partial \mathit{n}}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\cosh (jh/D_{\mathrm{e}}^{*})}}\right)}}\end{array}\end{array}$$

where $D_{\mathrm{e}}^{\ast }=\sqrt{\mathrm{2}{A}_z/f}$ is the traditional (steady) Ekman depth. In this steady-state solution, the first term represents the Ekman component (${\mathit{u}}_{\mathrm{e}}^{\ast }$), while the second term describes the geostrophic current and its bottom Ekman component (${\mathit{u}}_{\mathrm{g}}^{\ast }$). The solution of time-dependent Ekman theory, Eqs. (4) and (5), are almost identical to that of steady Ekman theory, Eq. (6), except that f in Eq. (6) is replaced by f−w_n. Dynamics of the time-dependent Ekman theory without the pressure gradient term, Eq. (4), was discussed in the study by Lilly and Elipot (2021). The time-dependent Ekman depth $D_{{\mathrm{e}}_n}$ becomes a function of the forcing frequency w_n. When w_n≪f (low frequency forcing), the dynamics are dominated by Coriolis force, governing classic steady Ekman and geostrophic currents. Conversely, when w_n≫f (high frequency forcing), Coriolis force becomes negligible relative to the inertia, and the system responds primarily with acceleration in the direction of the forcing. This response cannot be described by the classic steady Ekman theory, but is represented in the time-dependent solution, which better explains surface current responses under both low and high forcing frequency.

It is worth noting that resonance occurs when the frequency approaches the Coriolis frequency. In the limit of w_n→f (resonance condition), Eqs. (4) and (5) converge, respectively, to

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathit{u}}_{\mathrm{e}}={\displaystyle \frac{{\widehat{\mathit{\tau }}}_n}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{h}{A_z}}{e}^{-ift}\text{(8)}& {\displaystyle }{\mathit{u}}_{\mathrm{g}}=-{\widehat{\mathit{\eta }}}_n^{\prime }{\displaystyle \frac{g{h}^{\mathrm{2}}}{\mathrm{2}{A}_z}}{e}^{-ift}\end{array}$$

which are the exact solutions to Eqs. (2) and (3) in case of w_n=f, respectively. Consequently, when w_n→f, the velocity amplifies to specific constants rather than u→∞ because of finite depth and the frictional bottom boundary condition.

3 Data and methods

3.1 Data

In this study, the velocity fields were derived using altimetric sea surface height and the corresponding geostrophic current fields provided by the Copernicus Marine Environment Monitoring Service (CMEMS), as well as wind stress fields from the ERA5 reanalysis dataset provided by the European Centre for Medium-Range Weather Forecasts (ECMWF). This configuration is identical to that of Choi et al. (2023) and Globcurrent velocity fields (Rio et al., 2014). The altimetric sea surface height and geostrophic current fields were provided at a spatial resolution of 0.125° × 0.125° and a daily temporal resolution. The altimetry measurements were processed using low-pass filtering to remove high frequency components, such as tidal and inertial oscillations, and then gridded using an optimal interpolation technique (Saraceno et al., 2008; Pujol et al., 2016; Taburet et al., 2019). The wind stress fields were provided at a spatial resolution of 0.25° × 0.25° and an hourly temporal resolution. Park et al. (2018) demonstrated that wind fields from the ECMWF reanalysis dataset show the highest comparability with in-situ observations around the Korean Peninsula. Ocean bathymetry, with a spatial resolution of 15 arc-seconds, was obtained from the General Bathymetric Chart of the Oceans (GEBCO) and is used as input for the steady and time-dependent Ekman theory models. TOPEX/Poseidon global tidal model (TPXO9) is a high-resolution global tidal model from Oregon State University used to estimate tidal elevation. TPXO9 has a horizontal resolution of 0.033° × 0.033°. The TPXO9 is a sophisticated inverse model from Oregon State Tidal Inversion Software (OTIS) based on numerous altimetric and tide-gauge measurements (Egbert and Erofeeva, 2002).

In addition, we validate the velocity fields using in-situ velocity measurements from buoy (Eulerian) and surface drifter (Lagrangian) measurements. Buoy data provided by the Korea Hydrographic and Oceanographic Agency (KHOA) (available at http://www.khoa.go.kr/oceangrid, last access: 30 May 2025), which offers real-time and predictive information on oceanographic parameters (red triangle in Fig. 1). The KHOA buoy dataset covers 2012 to the present; however, due to missing and erroneous data before 2015, the observations from 2015 to 2023 were used in this study. Outliers were removed by excluding values beyond five standard deviations. Drifter data were obtained from the Jellyfish Information System (yellow lines and circles in Fig. 1). The drifter observations were conducted from 2015 to 2017 during summer seasons. We derived drifter velocities from the position changes by computing successive displacements at each timestep.

3.2 Methods

3.2.1 Surface current reconstruction

This study aims to estimate surface currents with high temporal resolution (hourly). To address this objective, the reconstruction process was designed to align the grids of the ERA5 wind stress and altimetric sea surface height fields. Therefore, ERA5 wind stress was spatially interpolated onto the grid of the altimetric sea surface height (η_altimetry), and η_altimetry was then temporally interpolated using the Fourier transform interpolation, which conserves the original spectrum of η_altimetry. It is worth noting that spiky errors are caused when linear interpolation is used in time (not shown here), so it is particularly meaningful to conserve the original spectrum of the data during temporal interpolation. For the pressure gradient forcing, as mentioned in Sect. 3.1, high frequency signals were removed in the η_altimetry during the gridding processes. To alleviate the absence of high frequency components in η_altimetry, sea surface height was reconstructed by simple summation of η_altimetry and tidal elevation η_tide based on TPXO9, so $\mathit{\eta }={\mathit{\eta }}_{\mathrm{altimetry}}+{\mathit{\eta }}_{\mathrm{tide}}$. This is an ad hoc approach caused by lack of spatiotemporal sampling ability of observations. It must be noted that the reconstructed η incorporates both low frequency and several high frequency components (tidal harmonics) but cannot consider entire spectrum of real sea surface height (will be discussed in Sect. 5.3).

The surface currents are estimated based on Eqs. (4) and (5) following the below steps: ${\widehat{\mathit{\tau }}}_n$ and ${\widehat{\mathit{\eta }}}_n^{\prime }$, depending on the frequency w_n, are obtained by Fourier transforming time series of τ^s and $\partial \mathit{\eta }/\partial \mathit{n}$. And then, $H_{{\mathrm{e}}_n}$ and $H_{{\mathrm{g}}_n}$ are multiplied by ${\widehat{\mathit{\tau }}}_n$ and ${\widehat{\mathit{\eta }}}_n^{\prime }$, respectively. Finally, inverse Fourier transforms were performed to $H_{{\mathrm{e}}_n}{\widehat{\mathit{\tau }}}_n$ and $H_{{\mathrm{g}}_n}{\widehat{\mathit{\eta }}}_n^{\prime }$, that derive the surface current components u_e and u_g, respectively, in the time domain.

3.2.2 Optimal eddy viscosity parameter estimation

It is worth noting that more generalized analytical solutions resolving variability in the eddy viscosity are known (Elipot and Gille, 2009; Wenegrat and McPhaden, 2016; and Lilly and Elipot, 2021). Nevertheless, we used the solution assuming constant eddy viscosity, Eqs. (4) and (5), because this study does not aim to examine the variability of eddy viscosity in the region, but to study influence of the inertia that was not considered in previous diagnostic velocity fields (Bonjean and Lagerloef, 2002; Rio et al., 2014; Choi et al., 2023). As mentioned above (Sect. 2), constant eddy viscosity was assumed for simplicity and consistency with the previous study but it cannot consider the seasonal variation of the eddy viscosity in the Yellow Sea (Hwang et al., 2014) influenced by the East Asian monsoon. To alleviate this issue, it is assumed that the eddy viscosity is a piecewise constant function, so the eddy viscosity is defined to be same constant value within a month. The value of the month was inversely estimated through sensitivity experiments by optimizing the eddy viscosity in u_e to minimize discrepancies with the buoy observations. The observation for u_e is defined as observed velocities minus altimetric geostrophic and tidal current components u_tide that was calculated by T_TIDE algorithm (Pawlowicz et al., 2002).

The eddy viscosity was optimized within the tested range of 0.001–0.2 m² s⁻¹ with 0.001 m² s⁻¹ interval. For each month, the analytical solution was computed using the predefined eddy viscosity. The Root-Mean-Square Error (RMSE) was calculated to measure the difference between the observed and estimated velocities and was then used as the objective function for optimizing the best-fit eddy viscosity estimation. Median values of the eddy viscosity estimated within each month (Fig. 3) are used to calculate velocity fields in this study.

4 Results

4.1 Estimation of eddy viscosity based on the optimization approach

Figure 2 presents the optimized vertical eddy viscosity (A_z) values, revealing distinct seasonal variability. The results indicate higher A_z values during winter (December–February), with a median value of approximately 0.17 m² s⁻¹ and a maximum reaching 0.2 m² s⁻¹. In contrast, summer (June–August) exhibits lower A_z values, with a median around 0.02 m² s⁻¹ and minimal fluctuations. This seasonal variation is closely linked to the East Asian monsoon system, where stronger winter north-westerly winds enhancing vertical mixing, rather than weaker southerly winds during summer (Hwang et al., 2014). The results align with previous parameterization studies (Csanady, 1981; Cushman-Roisin and Beckers, 2011), which suggest that Ekman depth – and consequently A_zi s strongly modulated by wind stress. Note that Csanady (1981) and Cushman-Roisin and Beckers (2011) parameterized the Ekman depth ($D_{\mathrm{e}}^{\ast }$) as a function of wind stress, where eddy viscosity can be expressed as $A_z=f{D}_{\mathrm{e}}^{\ast \mathrm{2}}/\mathrm{2}$. Additionally, Choi et al. (2023) showed that this parameterization effectively applies to the Yellow Sea, reinforcing the robustness of these estimates. The performance of this optimization approach was similar with Choi et al. (2023) using the Ekman depth parameterization, that will be shown and further assessed in Sect. 4.2 through comparisons with observational data.

Figure 2 Box plots for the optimal eddy viscosity parameter for each month, minimizing the Root-Mean-Square Error (RMSE) between surface velocities derived from Eq. (3) and buoy velocities. The whiskers (plus symbols) indicate variability outside the upper and lower quartiles, and the blue line in the blue box shows the median value, that were used as representative values to calculate surface current.

4.2 Verification of time-dependent Ekman theory model with observation and comparison with the steady Ekman theory model

The comparison between velocities estimated using the time-dependent Ekman theory model and observed velocities is presented in Fig. 3. Figure 3a and b show the comparison with buoy-observed velocities, showing correlation coefficients (R) of 0.70 and 0.66 for zonal and meridional components, respectively, and corresponding RMSE values of 0.15 m s⁻¹ for both components. Figure 4c and d show the comparison results with drifter-observed velocities. In this case, the correlation coefficients were higher than those for the buoy observations, with 0.78 for the zonal component and 0.82 for the meridional component. The RMSE for each component are 0.23 and 0.24 m s⁻¹, respectively. As shown in Fig. 1, the buoy observations (red triangle) are located in deeper ocean than the drifter trajectories (yellow circles and lines) concentrated in the shallower central and western Yellow Sea. The buoy location can be characterized by a stronger influence of the Kuroshio current branch (stronger geostrophic current component; Liu and Gan, 2012) and a weaker tidal current component (Hwang et al., 2014). This difference in measurement locations suggests that different hydrodynamic processes influence each observation dataset. Despite these regional differences, the model performs consistently across both observation locations, indicating its capability to represent surface current variability in the Yellow Sea. It is worth noting that the correlation with drifter observations is notably higher than that with buoy observations (Fig. 3), likely due to the stronger tidal components, which are periodic so relatively easier to estimate accurately using a reliable η_tide.

Figure 3 Comparison between time-dependent Ekman velocities with observations. Panels (a) and (b) show comparisons with buoy-observed velocities for the zonal and meridional components, respectively. Panels (c) and (d) display comparisons with drifter-observed velocities for the zonal and meridional components. The color scale indicates data density in logarithmic scale.

Figure 4 Comparison between observation (buoy and drifters) and velocity fields based on the time-dependent (left panels) and steady theory model (right panels). The upper panels compare observation with (a) time-dependent and (b) steady Ekman model velocities in hourly temporal resolution. The lower panels show comparisons for daily-averaged observation and (c) time-dependent and (d) steady models. The color scale indicates data density in logarithmic scale. Circles and triangles indicate zonal and meridional velocity components, respectively.

The steady Ekman theory model (u^∗) based on Choi et al. (2023) and the time-dependent Ekman theory model (u) proposed in this study were compared to assess their performance in reconstructing surface currents (Fig. 4). The results indicate that u much more accurately represents observed surface currents, particularly in capturing high frequency variability, whereas u^∗ shows notable limitations in this regard. This is reflected in the correlation coefficient of 0.78 and RMSE of 0.18 m s⁻¹, whereas the steady model shows a significantly lower correlation of −0.08 and a higher RMSE of 0.43 m s⁻¹ (Fig. 4a and b). Uncorrelated signals between the steady model and observation are because high frequency signals are not properly resolved by the governing equation ignoring inertia (Fig. 4c and d), discussed in Sect. 2.2. It is worth noting that the previous studies considered only low frequency signal (Bonjean and Lagerloef, 2002; Rio et al., 2014; Choi et al., 2023). These results indicate that the u properly represents high frequency variability in surface current, while the steady model has notable limitations.

Despite limitation of u^∗ in capturing high frequency variability, u^∗ remains effective in reconstructing low frequency surface current. Figure 4c and d compares the daily-averaged velocity from each model (u and u^∗) and that from observation. The steady model u^∗ shows a correlation of 0.53 and RMSE of 0.10 m s⁻¹, while time-dependent model u shows a correlation of 0.62 and RMSE of 0.09 m s⁻¹ that is slightly improved than steady model. This similarity is expected, as the time-dependent Ekman model (4) and (5) mathematically converges to the steady Ekman model (6) when w_n≪f. This explains why the steady Ekman theory can be used to simulate transport of the surface floating disturbances in the Yellow Sea (Choi et al., 2018, 2023), of which time scale is order of a month much longer than $\mathrm{1}/f$, so expected to be governed by the low frequency dynamics (typical steady Ekman dynamics and predominant geostrophic balance). Consequently, the ability of the time-dependent Ekman model to represent both low and high frequency components stands out as the most significant. While the steady model can resolve only low frequency variability. We expect that the approach in this study can be applicable to the regions like the Yellow Sea where stress effect is dominant.

5 Discussions

5.1 Limitation of the steady model

This study demonstrates that the steady Ekman theory model fails to represent high frequency variability in the observed velocities, whereas the time-dependent theory model provides a closer representation of the observations in both high and low frequencies. Figure 5 presents variance ellipses showing spatial patterns of surface current variability in observations (black solid ellipses), the steady theory models (u^∗; blue ellipses), and the time-dependent theory model (u; red dashed ellipses) within 0.5° × 0.5° bins. The colormap within the grid represents the number of observations. The variance ellipse represents current variability: its orientation indicates the dominant direction of variability, and the lengths of the major and minor axes correspond to the variances along those directions. These ellipses reveal the differences between observed and modeled surface current variability. The ellipses estimated from the u^∗ deviate from the observed patterns in both the magnitude and directional variability. This deviation becomes more pronounced toward the western coast of Korean Peninsula, where the ellipses appear enlarged and misaligned. In contrast, u produces ellipses that more closely match observed ellipses, both in size and orientation, suggesting its improved capability in representing current variability. It is worth noting that the drifter observations are predominantly collected during the summer, leading to a seasonal bias, but the buoy observations include data from all seasons. The observed variance ellipses exhibit patterns that closely resemble the tidal characteristics, particularly the M2 tidal component, as reported by Hwang et al. (2014). This consistency further reinforces the dominant influence of tides on current variability in the study area. This suggests that u effectively simulates tidal currents from the reconstructed η, shown in Sect. 3.2.1. Importantly, this underscores the role of the inertial term in representing tidal currents, as it enables the model to better represent high frequency variations.

Figure 5 Spatial distribution of observed (black solid line) and modeled variance ellipses by the steady Ekman model (u^∗; blue lines) and the time-dependent Ekman model (u; red dashed lines). The color scale indicates the number of observations in each bin. All variance ellipses are scaled consistently to represent the velocity magnitudes and directions within the respective grids. The steady model fails to resolve proper current variability, on the other hand, the variabilities are successfully considered by the time-dependent model.

Figure 6 shows Taylor diagrams for both the zonal and meridional velocity components, comparing the contributions of each current component, derived from u and u^∗, to observations. The diagram provides a statistical evaluation of model performance using three key metrics: correlation, Standard Deviation (STD), and Root-Mean-Square-Difference (RMSD). The black lined circle markers represent observed velocities from drifters and buoys. Red markers represent components derived from u, while blue markers indicate those from u^∗. When the wind-driven components from each model are compared, the time-dependent model (u_e; red squares in Fig. 7) shows a higher contribution than the steady models (${\mathit{u}}_{\mathrm{e}}^{\ast }$; blue squares in Fig. 6) so time-dependent model shows slightly improved agreement with observations. This implies that incorporating the inertial term enhances the model ability to represent wind-driven currents, particularly in response to high frequency wind forcing. The contribution of the time-dependency in the wind-driven component can be locally significant (e.g., typhoon event), that will be discussed in Sect. 5.2. The most significant limitation of the steady model is uncorrelated tidal components from ${\mathit{u}}_{\mathrm{tide}}^{\ast }$ (blue downward-pointing triangles), on the other hand, the time-dependent model successfully resolves the tidal components u_tide (red downward-pointing triangles). This alludes that the poor performance in the hourly steady model (Fig. 4b) is due to misrepresentation for tides. In Fig. 6, u_g (red upward-pointing triangles) is located closer to the observations with reasonable correlation coefficients (0.69 and 0.72 for the zonal and meridional components, respectively). In contrast, ${\mathit{u}}_{\mathrm{g}}^{\ast }$ (blue upward-pointing triangles) is uncorrelated with observation (−0.18 and −0.09 for the zonal and meridional components, respectively) due to the issues in tides. These results explain the limitation of the steady model ignoring inertia that plays a key role in tides.

Figure 6 Taylor diagram between observations (drifters and buoy, black lined circles), time-dependent Ekman components (red markers) and steady Ekman components (blue markers) for the zonal components (a) and meridional component (b). u (red circles) represents the total velocity, u_g (red upward-pointing triangles) indicates the pressure-gradient component, u_tide (red downward-pointing triangles) indicates tidal component, and u_e (red squares) indicates wind-driven component, all from the time-dependent Ekman model. u^∗ (blue circles) represents the total velocity, ${\mathit{u}}_{\mathrm{g}}^{\ast }$ (blue upward-pointing triangles) indicates the pressure-gradient component, ${\mathit{u}}_{\mathrm{tide}}^{\ast }$ (blue downward-pointing triangles) indicates tidal component, and ${\mathit{u}}_{\mathrm{e}}^{\ast }$ (blue squares) indicates wind-driven component, all from the steady Ekman model.

Figure 7 Taylor diagram comparing observed daily-averaged surface velocities (drifters and buoy, black-lined circles) with time-dependent Ekman model components (red markers) for (a) zonal and (b) meridional. u (red circles) represent the time-dependent Ekman total velocity. u_e (red squares) indicates wind-driven component, u_g (red upward-pointing triangles) indicates the pressure-gradient component, u_tide (red downward-pointing triangles) indicates tidal component. u_alt (yellow upward-pointing triangles) indicates the geostrophic current component derived from altimetry (η_altimetry).

It is worth noting that the contribution of each surface current component, explaining variance of the observation, can be quantified by the squared values of the correlation in Fig. 6. The pressure-gradient component (red upward-pointing triangles in Fig. 6), especially, tidal component (red downward-pointing triangles in Fig. 6), is the most dominant contributor to surface current variability, and other components have relatively minor role in explaining the variance of the observation. However, in the longer time scales than the periods of tidal components (e.g., daily averaged velocity fields; Fig. 7), contribution of the tidal components (red downward-pointing triangles in Fig. 7) dramatically decreases and those of the other components, such as geostrophic current (u_alt; yellow upward-pointing triangles in Fig. 7) and wind-driven current components (red squares in Fig. 7), become considerable. Note that the geostrophic current component u_alt is defined by the geostrophic balance equation and estimated by only altimetric sea surface height.

This is not surprising because the tidal components are purely periodic sinusoidal functions, whose means over their periods are intrinsically zero. Consequently, contribution of the velocity components varies with the time scale. For currents associated with the transport of the floating substances across the Yellow Sea, with a time scale on the order of months, non-tidal low frequency current components (geostrophic and Ekman currents) play the dominant role (Fig. 7), even though tides outperform the other components on the short time scale (Fig. 6). Based on the correlation shown in Fig. 7, Ekman current components contribute slightly more than geostrophic currents in the Yellow Sea, that is consistent with Choi et al. (2018).

5.2 Damped inertial oscillation caused by a typhoon

It is shown that tidal currents are the most dominant component in the variation of surface currents, but other components contribute to reconstructing more realistic surface currents in the study area (Fig. 6). In this section, cases where other components (wind- and friction-driven components including inertial oscillations) become locally significant are shown and discussed. Figure 8 presents the time series of surface current responses to Typhoon Chanthu in 2021 (KMA, 2022) at the buoy site, including wind stress (Fig. 8a and b), observed velocities from the buoy (black solid lines), and diagnostic velocities estimates derived from approaches with (time-dependent theory model u; red solid lines) and without (steady theory model u^∗; blue solid lines) the inertial term.

Figure 8 Time series comparison of buoy-observed velocities (black lines) and modeled velocities from the time-dependent (red lines) and steady (blue lines) Ekman theory models during Typhoon Chanthu in September 2021. The period when the typhoon passed over the buoy observation location is highlighted by blue shading. Upper panel (a) and (b) show time series of wind stress at the buoy site. Middle panel (c) and (d) show direct comparison between observed and modeled velocities. Lower panel (e) and (f) show comparison of de-tided velocity time series, where the red-dashed line represents the arbitrary exponential function using ${\mathit{\gamma }}_{\mathrm{0}}=A_z{\left(\mathit{\pi }/\left(\mathrm{2}h\right)\right)}^{\mathrm{2}}$ , with A_z corresponding to the September eddy viscosity shown in Fig. 2.

The M2 tidal constituent is dominant (Fig. 8c and d) and the time-dependent model properly captures the signals from M2 tidal component, on the other hand, steady model exhibits poor phase agreement with observation. In addition, considerable wind-driven component is also shown during the typhoon event (blue-shaded area in Fig. 8). Between September 16 and 20, which includes the typhoon-affected period, non-tidal components became even more dominant than tidal component; correlation coefficients between observed (black solid line in Fig. 8c and d) and reconstructed currents without tides (red solid line in Fig. 8e and f) reached 0.73, which is notably higher than correlation of tidal component, 0.45, based on T_TIDE. This indicates that, although tide generally dominates, the other (wind-driven) components play a significant role under strong wind events like typhoons.

Removing the tidal component from observation (u − u_tide) reveals clear identification of typhoon-induced inertial oscillations after the typhoon event (Fig. 8e and f). The steady model fails to represent the inertial oscillation, as it only reflects low frequency variations. This limitation is evident as the steady model neglecting inertia cannot consider the inertial oscillation, governed by balance between the Coriolis force and inertia. In contrast, the time-dependent model clearly resolves the inertial oscillation and also its damping (Fig. 8e and f). Following the abrupt perturbation from the typhoon, inertial oscillations are developed as the current responds to the typhoon and then gradually decay over time, reflecting dissipative processes in post-storm currents (Fig. 8e and f). Owing to the simplicity of the governing Eq. (1), the decay rate can be analytically determined by expanding boundary value problem (1) to initial-boundary value problem with initial condition $\mathit{u}\left(z,t=\mathrm{0}\right)={\mathit{u}}_{\mathrm{0}}\left(z\right)$ where u₀(z) is arbitrary initial velocity profile. Detailed mathematical derivations are in the Appendix. Considering an arbitrary initial condition to Eq. (1) yields additional velocity component u_i associated with the initial condition is given by

$$\begin{array}{}\text{(9)}& {\mathit{u}}_{\mathrm{i}}=\sum _{n=\mathrm{0}}^{\mathrm{\infty }}{e}^{-ift}{e}^{-{\mathit{\gamma }}_{n}t}{\mathit{U}}_n\cos \left(k_{n}z\right)\end{array}$$

where ${\mathit{\gamma }}_n=A_z{k}_n^{\mathrm{2}}$ indicates the decay rate of nth mode and ${\mathit{u}}_{\mathrm{i}}(zt=\mathrm{0})={\mathit{u}}_{\mathrm{0}}\left(z\right)-{\mathit{u}}_{\mathrm{e}}\left(z,\mathrm{0}\right)-{\mathit{u}}_{\mathrm{g}}\left(z,\mathrm{0}\right)=\sum _{n=\mathrm{0}}^{\mathrm{\infty }}{\mathit{U}}_n\cos \left(k_{n}z\right)$. The inertial oscillations and attenuations are described by e^−ift and $e^{-{\mathit{\gamma }}_{n}t}$ in Eq. (9), respectively. It is worth noting that the decay rate γ_n depends on the wavenumber k_n of vertical velocity profile. Therefore, it is difficult to determine precise decay rate of the surface current without information about the vertical velocity profile. Nevertheless, the slowest (minimum) decay rate occurs when n=0, so ${\mathit{\gamma }}_{\mathrm{0}}=A_z{k}_{\mathrm{0}}^{\mathrm{2}}=A_z{\left(\mathit{\pi }/\left(\mathrm{2}h\right)\right)}^{\mathrm{2}}$, which is not only applicable to surface current but also the most important value because its reciprocal indicates the longest time scale for the signals to arrive. In Fig. 8e and f, the red dashed line indicates exponential decay using the slowest decay rate and A_z value for September (shown in Fig. 2). It successfully explains decay of the inertial oscillation shown in the time-dependent theory model u.

The time-dependent model shows slightly stronger damping compared to the observation (Fig. 8e and f). Despite the overestimated damping, it must be highlighted that the time-dependent model successfully reconstructs the overall surface currents in the Yellow Sea (Fig. 3). The time-dependent model plausibly captures the damped inertial oscillation, that is totally ignored by the steady model (Fig. 8e and f). The simplicity of theory model provides valuable insight into the attenuation of inertial oscillations. It can be explained how vertical eddy viscosity and bottom friction attenuates signals that occurred in the past: the decay rate depends on the wavenumber of vertical velocity profile but the slowest damping rate can be applicable to surface current. Ultimately, this study underscores the advancement of diagnostic approaches in high frequency processes with the inclusion of the inertial term.

5.3 Resolved and unresolved spectrums by the time-dependent Ekman theory model

Figure 9 presents the spectrum of zonal (upper panels) and meridional (lower panels) velocity components. The left panels (Fig. 9a and c) show buoy-observed velocities while the right panels (Fig. 9b and d) show reconstructed velocities from the time-dependent Ekman model. The red vertical dashed line indicates the local Coriolis frequency ($f/\mathrm{2}\mathit{\pi }$), and the labeled peaks correspond to dominant tidal constituents: the diurnal (O1 and K1), and semidiurnal (M2, S2 and N2) tidal components. To highlight the dominant spectral patterns, band-averaged spectrum (bold solid lines) is also shown.

Figure 9 Spectrum amplitude of observed (left panels; a and c) and modeled (right panels; b and d) velocity for the zonal (upper panels; a and b) and meridional (lower panels; c and d) components. Blue lines represent the raw spectra. The bold solid lines represent the band-averaged spectra with 100-point window from observations (bold black lines in left panels) and the time-dependent model (bold red lines in right panels), respectively. Major tidal constituents (O1, K1, M2, S2, and N2) are indicated by yellow markers, and the vertical dashed red line denotes the local Coriolis frequency ($f/\mathrm{2}\mathit{\pi }$).

The time-dependent model (right panels) clearly represents tidal components, showing distinct spectral peaks at specific constituents in observation (left panels). In addition to low-frequency signals, the model spectrum also represents the resonant response near the inertial oscillation frequency. As discussed by Kim and Crawford (2014) and Elipot and Lumpkin (2008), resonance occurs when the frequency of external forcing is close to the $f/\mathrm{2}\mathit{\pi }$, amplifying the ocean's response to the forcing. This spectral peak near the $f/\mathrm{2}\mathit{\pi }$ is shown in the observation (black solid lines in left panels) and resolved by the model (red solid lines in right panels). Consequently, the time-dependent Ekman theory model u accurately resolves the resonant signals near the $f/\mathrm{2}\mathit{\pi }$ although the model underestimates high frequency signals.

This underestimation can be attributed to the uncertainties in the reconstruction of η using the ad hoc approach shown in Sect. 3.2.1. Since the pressure-gradient component is calculated by the reconstructed η, its accuracy directly impacts spectral characteristics. The spectrum of $\partial \mathit{\eta }/\partial x$ (Fig. 10a) and $\partial \mathit{\eta }/\partial y$ (Fig. 10b) reveals a limitation in the reconstruction approach, that considers only low frequency signals estimated from altimetry and several high frequency modes corresponding tidal harmonics. The simple summation of the low frequency component η_altimetry and the tidal component η_tide fails to fully capture high frequency variability and misses many intermediate and high frequency modes (Fig. 10). Consequently, the ad hoc η reconstruction is not sufficient to resolve entire spectrum. We expect that this limitation can be alleviated as availability for η observation increases, for example, by using blended sea surface height field by tide gauge data. Nonetheless, this study remains valuable as it demonstrates the ability to capture high frequency variability including tides and inertial oscillations, which are predominant in the study area.

Figure 10 Spectrum amplitude of the reconstructed sea surface height gradients: (a) $\partial \mathit{\eta }/\partial x$ and (b) $\partial \mathit{\eta }/\partial y$. Blue lines indicate raw spectra and black bold lines represent band-averaged spectra using a 100-point window. Major tidal constituents (O1, K1, M2, S2, and N2) are marked with yellow triangles, and the vertical red dashed line indicates the local Coriolis frequency ($f/\mathrm{2}\mathit{\pi }$). These spectra show limitation of the simple sea surface height reconstruction that resolves only low frequency signals from altimetry and several high frequency modes from tides. Many intermediate and higher frequency signals are not resolved.

6 Conclusion

The usefulness of the diagnostic surface velocity field was reported by many previous studies (Bonjean and Lagerloef, 2002; Rio et al., 2014; Choi et al., 2018, 2023), but steady models considering only geostrophic and Ekman current components have been used. To the best of our knowledge, this is one of the first studies developing the diagnostic velocity field considering the inertial terms in the momentum equations, along with Ubelmann et al. (2025). Although their response function approach, in time domain, is equivalent to our transfer function method, in frequency domain, based on the convolution theorem, there are key differences. To be specific, they used data-driven empirical modeling approaches for the response function, whereas we adopted a stricter formulation based on the analytical solution to the momentum equation. It is worth noting that Ubelmann et al. (2025) used a slab model, which is vertically averaged, for their rationale and discussion. On the other hand, the governing equations of this study are not vertically averaged. The fact that we generalized the time-dependent Ekman theory to incorporate pressure gradient term, thereby resolving tides, is our unique contribution.

This study evaluates the performance of the time-dependent Ekman theory model and compares the model with the steady Ekman theory model (explicitly and implicitly used by the previous studies) in the Yellow Sea. The time-dependent model successfully considers not only low frequency velocity components (geostrophic and Ekman currents) but also high frequency components (inertial oscillations and tides). The tidal currents outperform the other velocity components, such as geostrophic and Ekman current components, at short time scales described by the hourly diagnostic velocity fields (Fig. 6). On the other hand, the contribution of tides in longer time scales becomes negligible because averaged values over their periods are zero (Fig. 7). This implies that the dominant current component transporting a subject is determined by the time scale: surface floating substances crossing the Yellow Sea, whose time scale is order of months, are dominantly determined by low frequency velocity components including geostrophic and Ekman current components. Once position of the subject should be forecasted as quickly as possible on short time scale (e.g., search and rescue), tides become considerably important. This study also highlights the role of the simple analytical solutions improving our understanding for the high frequency surface current variability and provides a pathway for future advancements in surface current estimation.

Appendix A

This appendix provides the mathematical formulation of the expanded initial-boundary value problem given by

$$\begin{array}{}\text{(A1a)}& {\displaystyle \frac{\partial \mathit{u}}{\partial t}}+fi\mathit{u}=-g{\displaystyle \frac{\partial \mathit{\eta }}{\partial \mathit{n}}}+A_z{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{u}}{\partial z^{\mathrm{2}}}},\text{(A1b)}& {\displaystyle }{\left.A_z{\displaystyle \frac{\partial \mathit{u}}{\partial z}}\right|}_{z=\mathrm{0}}={\displaystyle \frac{{\mathit{\tau }}^{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{0}}}},\text{(A1c)}& {\displaystyle }{\left.\mathit{u}\right|}_{z=-h}=\mathrm{0}\text{(A1d)}& {\displaystyle }{\left.\mathit{u}\right|}_{t=\mathrm{0}}={\mathit{u}}_{\mathrm{0}}\left(z\right).\end{array}$$

It must be noted that the solution (4) and (5) discussed in Sect. 2 are from the boundary value problem (1), so there is no consideration of specific initial condition. The only difference from Eq. (1) is the inclusion of the initial condition (A1d), which defines an arbitrary initial velocity profile. The total velocity u is decomposed as $\mathit{u}={\mathit{u}}_{\mathrm{b}}+{\mathit{u}}_{\mathrm{i}}$, where u_b and u_i are the component associated with forcings and initial condition, respectively. The governing equation and boundary condition for u_b are identical to Eq. (1) so given by

$$\begin{array}{}\text{(A2a)}& {\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{b}}}{\partial t}}+fi{\mathit{u}}_{\mathrm{b}}=-g{\displaystyle \frac{\partial \mathit{\eta }}{\partial \mathit{n}}}+A_z{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{u}}_{\mathrm{b}}}{\partial z^{\mathrm{2}}}},\text{(A2b)}& {\displaystyle }{\left.A_z{\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{b}}}{\partial z}}\right|}_{z=\mathrm{0}}={\displaystyle \frac{{\mathit{\tau }}^{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{0}}}},\text{(A2c)}& {\displaystyle }{\left.{\mathit{u}}_{\mathrm{b}}\right|}_{z=-h}=\mathrm{0}\end{array}$$

Therefore, solution to Eq. (A2) is the time-dependent theory model ${\mathit{u}}_{\mathrm{b}}={\mathit{u}}_{\mathrm{e}}+{\mathit{u}}_{\mathrm{g}}$ where u_e and u_g are given by Eqs. (4) and (5). Subtracting Eqs. (A1) and (A2) yields the governing equation and boundary conditions for u_i given by

$$\begin{array}{}\text{(A3a)}& {\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{i}}}{\partial t}}+fi{\mathit{u}}_{\mathrm{i}}=A_z{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{u}}_{\mathrm{i}}}{\partial z^{\mathrm{2}}}},\text{(A3b)}& {\displaystyle }{\left.A_z{\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{i}}}{\partial z}}\right|}_{z=\mathrm{0}}=\mathrm{0},\text{(A3c)}& {\displaystyle }{\left.{\mathit{u}}_{\mathrm{i}}\right|}_{z=-h}=\mathrm{0}.\end{array}$$

The initial condition for u_i can be obtained by rewriting Eq. (A1d) and given by

$$\begin{array}{}\text{(A4)}& {\mathit{u}}_{\mathrm{i}}(zt=\mathrm{0})={\mathit{u}}_{\mathrm{0}}\left(z\right)-{\mathit{u}}_{\mathrm{b}}\left(z,\mathrm{0}\right).\end{array}$$

Consequently, Eqs. (A3) and (A4) become a problem without forcing terms, ${\mathit{\tau }}^{\mathrm{s}}/{\mathit{\rho }}_{\mathrm{0}}$ and $-g\partial \mathit{\eta }/\partial \mathit{n}$. Solution to Eq. (A3) with the initial condition (A4) is given by Eq. (9) and can be obtained by typical approach to solve the heat equation.