Geostrophic Spirals Generated by the Horizontal Diffusion of Vortex Stretching in the Yellow Sea

Current spirals exist in the ocean. They are important for understanding ocean dynamics and the transport of chemicals, contaminants (such as radioactive waste and hypertoxic materials) and nutrients. There are two types of spirals: one is the classic Ekman spiral (Ekman, 1905; Price et al., 1987) in the Ekman layer, and the other is the beta spiral (Stommel and Schott, 1977) in the geostrophic interior of the ocean. These two kinds of spirals are controlled by different dynamics. The Ekman spiral is the balance between the local surface wind stress and vertical eddy viscosity, while the beta spiral is primarily generated by large-scale vorticity balance.

The beta spiral is a terminology that describes the rotation of the geostrophic velocity with increasing depth in the ocean interior. Historically, however, before the beta spiral was discovered, the law of parallel solenoids (Neumann and Pierson, 1966) was the dominant theory. According to this theory, ocean currents should be parallel not only to isobars but also isopycnals. The departure from this theoretical framework of physical oceanography began with the pioneering work of (Stommel and Schott, 1977). They laid down the theoretical foundation of the beta spiral and demonstrated the existence of such spirals using hydrographic data obtained from the North Atlantic Ocean. Their results showed that vertical advection and stratification play important roles in determining the observed structure of beta spirals. The basic dynamic constraint for beta spirals is geostrophy. The basic assumption made in the study of the classic beta spiral is no cross-isopycnal mixing (e.g., potential vorticity is conserved, and the relative vorticity is much smaller). Therefore, the typical beta spiral only exists in the ocean interior, where current and diffusion processes are weak. In the interior ocean, upwelling/downwelling is associated with Ekman pumping due to basin-scale wind stress curl and the beta effect. Negative wind stress curl in the subtropical basin leads to clockwise beta spirals. Spiral calculations were extended from the North Atlantic Ocean to the global oceans (Schott and Stommel, 1978; Behringer and Stommel, 1980; Cui et al., 1991), similar to (Stommel and Schott, 1977). In a recent study by (Yang et al., 2018), it was found that the topographic beta effect associated with the Kuroshio intrusion-induced upwelling northeast of Taiwan Island is essential in causing velocity spirals, which plays a vital role in regulating onshore intrusion flow and the cross-shelf transport of energy and nutrients.

The beta spiral dynamics were extended to the so-called beta spiral method, which has been used to calculate absolute velocity via hydrographic data (Davis, 1978; Beheinfer, 1979; Schott and Zantopp, 1979; Bigg, 1985). This method also provides additional constraints for inverse models (Wunsch, 1978). When cooling is very strong and dominant, convective overturning and re-stratification can generate an anticlockwise cooling spiral. Cooling spirals have been documented in the subpolar ocean (Schott and Stommel, 1978) and western boundary current extensions (Spall, 1992). (Spall, 1992) showed that parcels in regions with cooling spirals are forced to cross from south to north, which opposes the classic beta spiral experienced in the interior of the subtropical ocean. Cooling spirals are dynamically important for understanding the subtropical gyre recirculation. Other processes can also induce geostrophic velocity spirals, such as horizontal and vertical eddy diffusion. However, such dynamic issues still remain unexplored. In fact, although velocity spirals in the open ocean have been well documented in previous studies, whether they exist in coastal seas is not clear. Our study here is focused on the characteristics and dynamics for the newly found spirals in the Yellow Sea (Fig. 1). In this paper, we answer the following two questions: (1) Can geostrophic spirals be observed in the coastal sea, where synoptic-scale ocean dynamics are complex?

(2) Can other processes (e.g., diffusion processes) cause geostrophic spirals besides those caused by Ekman pumping and the cooling effect?

In-situ observations were employed to demonstrate the dynamics of spirals in the Yellow Sea. Comprehensive hydrographic observations and objectively analyzed data were used to analyze the dynamics. The rest of this paper is organized as follows: In section 2, the observation program and geostrophic spiral theory are briefly introduced. In section 3, we use mooring observations in the western and northern Yellow Sea (Yu et al., 2010; Lin et al., 2011) to demonstrate the structure of the velocity spiral in coastal seas and explore its physical mechanisms. Section 4 diagnoses the dynamics of the horizontal velocity spirals in terms of vertical advection, wintertime surface cooling, and diffusion processes. Finally, we draw conclusions in section 5.

2.1.1. In-situ hydrodynamic observations

Three Chinese research vessels participated in simultaneous joint investigations. The First Institute of Oceanography, affiliated to the State Oceanic Administration and the Ocean University of China, coordinated a survey from 21 December 2006 to 14 February 2007. The moorings located at the sea bottom were equipped with a tide gauge, 300-KHz Teledyne RD Instruments (RDI) upward-looking acoustic Doppler current profiler (ADCP) from M1 to M5 (solid/empty black squares in Fig. 1), and a 250-KHz Sontek ADCP at M6 (black square in Fig. 1). Unfortunately, M1 and M3 (empty black squares in Fig. 1) were not retrieved normally due to frequent fishing activities in this region. Detailed information of the moorings is summarized in Table 1. The water depths from M1 to M5 range from 50 m to 73 m, with horizontal topographic gradients. The CTD casts in the western and northern Yellow Sea during the winter of 2006/07 were also conducted on the same cruise (black dots and lines in Fig. 1; Yu et al., 2010; Lin et al., 2011). For current measurements, the sampling interval was set to 5 min, with a vertical bin size of 2 m, while the vertically averaged bins for temperature and salinity were set to 1 m. All instruments were calibrated before and after the deployments with standard procedures.

2.1.2. High-resolution wind stress from data assimilation

The wind stress data based on the Weather Research and Forecasting Model output are used to calculate the surface Ekman pumping rate. The model assimilates wind obser-vations at land and coastal meteorological stations (Wu et al., 2011), with a 3-h temporal resolution and 0.1^° spatial resolution. These outputs have been operated to predict the weather and long-term environmental changes in the Yellow Sea and the East China Sea.

2.1.3. Surface heat flux

Surface heat fluxes from different heat flux datasets are subject to great bias and uncertainties (Song and Yu., 2013). To reduce the biases of heat fluxes, five full-flux datasets were used to obtain average surface fluxes, including three newly updated atmospheric reanalyses and two analyzed datasets (Table 2). The monthly mean net surface heat fluxes in December 2006 and January 2007 are used here. The three latest reanalyses include the NCEP Climate Forecast System Reanalysis (CFSR; Saha et al., 2010), the Modern Era Retrospective Analysis for Research and Applications (MERRA) from NASA (Bosilovich, 2008; Rienecker et al., 2011), and ECMWF's ERA-Interim project (Dee and Uppala, 2009). Two globally analyzed heat flux datasets are used in this study. The WHOI's Objective Analyzed Air-sea Fluxes (OAFlux) project (Yu and Weller, 2007) is combined with the International Satellite Cloud Climatology Project (ISCCP; Zhang et al., 2004) to obtain the net surface heat flux (Q _NET). Another analyzed net surface heat flux is derived from the version 2 forcing for the Co-ordinated Ocean-Ice Reference Experiments (CORE.2; Large and Yeager, 2009). Detailed information on these heat flux estimates has been described in our previous study (Song and Yu., 2013).

2.1.4. Mesoscale eddy diffusivity estimate

The mesoscale eddy diffusivities——namely, the lateral mixing by mesoscale eddies——were borrowed from (Klocker and Abernathey, 2014) (hereafter, KA14). These diffusivities are closely associated with the mixing length, eddy size and propagation relative to the mean flow. The eddies are obtained using AVISO altimetry data. The suppressed (i.e., without mean flow) and unsuppressed (i.e., with mean flow) diffusivities are 2752 and 2881 m² s^-1, respectively, in the region bounded by the red frame (Fig. 1). The suppression effect is not significant in the Yellow Sea, since the mean flow is weak (only ∼3 cm s^-1). The red framed region in Fig. 1, with the dashed lines, was chosen to locate the moorings in the middle and reduce the SSH uncertainties near the coast.

Figure 1.  The topography of the Yellow Sea and Bohai Sea and the mooring sites (labeled with black squares). Stations M1 (lost), M2 (normal), M3 (lost), M4 (normal), and M5 (normal) reside in the western Yellow Sea, and station M6 (normal) is located in the northern Yellow Sea. The solid lines and black dots represent parts of the synchronized hydrographic observations on the cruises via CTD casts. The contours show the bathymetry (m). The mean mesoscale eddy diffusivities (Klocker and Abernathey, 2014) were averaged in the region bounded by the dashed red line.

2.2.1 Basic dynamics

The dynamics of geostrophic spirals have been discussed in detail by (Stommel and Schott, 1977) and subsequent studies (Spall, 1992; Huang, 2010). Consider the steady density equation: \begin{equation} u\dfrac{\partial\rho}{\partial x}+v\dfrac{\partial\rho}{\partial y}+w\dfrac{\partial\rho}{\partial z}= k_{\rm h}\nabla^2\rho+k_{\rm v}\dfrac{\partial^2\rho}{\partial z^2}+\dot{q} , \ \ (1)\end{equation} where u, v and w represent the zonal, meridional and vertical velocities, respectively; ρ represents the density of sea water; k _h and k _v represent horizontal and vertical eddy diffusivities, respectively; and $\dot q$ represents the net surface heat flux. The thermal wind relation is: \begin{eqnarray} \dfrac{\partial u}{\partial z}&=&\gamma\dfrac{\partial\rho}{\partial y} ;\ \ (2)\\ \dfrac{\partial v}{\partial z}&=&-\gamma\dfrac{\partial\rho}{\partial x} . \ \ (3)\end{eqnarray} Here, γ=g/(fρ) is an assumed constant associated with the geostrophic parameters, g represents the gravitational acceleration, and f represents the Coriolis force. The velocity angle θ= arctan(u/v) ranges from 0^° to 360^°, with an eastward velocity zero. Its derivation in the vertical direction is \begin{equation} \dfrac{{\rm d}\theta}{{\rm d}z}=\dfrac{uv'-u'v}{u^2+v^2} . \ \ (4)\end{equation}

Substituting the thermal wind relation [Eqs. (2) and (3)] into the first and second terms on the lhs of Eq. (1), we obtain: \begin{equation} u\dfrac{\partial\rho}{\partial x}+v\dfrac{\partial\rho}{\partial y}=-\dfrac{1}{\gamma}\left(u\dfrac{\partial v}{\partial z} -v\dfrac{\partial u}{\partial z}\right)=-\dfrac{1}{\gamma}(u^2+v^2)\dfrac{{\rm d}\theta}{{\rm d}z} . \ \ (5)\end{equation}

Substituting the thermal wind relation into the first term on the rhs of Eq. (1), we obtain the horizontal diffusion of vortex stretching: \begin{equation} k_h\nabla^2\rho=-\dfrac{1}{\gamma}k_h\left[\dfrac{\partial}{\partial z}\left(\dfrac{\partial u}{\partial y}-\dfrac{\partial v}{\partial x}\right)\right] =-\dfrac{1}{\gamma}k_h\xi_z . \ \ (6)\end{equation}

Here, $\xi_z$ represent the vortex stretching. Replacing the terms in Eq. (1) with Eqs. (5) and (6) yields the control equation for the velocity spiral: \begin{equation} \underbrace{\dfrac{{\rm d}\theta}{{\rm d}z}}\limits_{\rm A}= \dfrac{\gamma}{u^2+v^2}(\underbrace{w\rho_z}\limits_{\rm B}-\underbrace{\dot{q}}\limits_{\rm C}+ \underbrace{\dfrac{1}{\gamma}k_h\xi_z}\limits_{\rm D}-\underbrace{k_v\rho_{zz}}\limits_{\rm E}) . \ \ (7)\end{equation} where ρ_z and ρ_zz represent the first and second derivative of the sea water density ρ. Physically, the velocity rotation (Term A) with depth is the result of the competition between terms on the rhs of Eq. (7), which includes the vertical advection of density (Term B), surface buoyancy changes due to heating/cooling (Term C), and the horizontal (Term D) and vertical (Term E) eddy diffusion processes. In this paper, the rotational budget of Eq. (7) was used to diagnose the mechanism that caused the velocity spirals.

2.2.2. Methods for dynamic calculations

The vertical advection w of Term A on the rhs of Eq. (7) is the vertical velocity, which can be estimated from the net Ekman pumping velocity in the surface/bottom boundary layers (w _EKS and w _EKB, respectively) and topography-induced upwelling (w _BFT).

The surface Ekman pumping rate (w _EKS) is associated with the surface wind stress and beta effect: \begin{equation} w_{\rm EKS}=\dfrac{1}{\rho f}{\rm curl}(\tau)+\beta\dfrac{\tau_x}{f^2\rho} , \ \ (8)\end{equation} where τ represents the wind stress according to (Large and Yeager, 2009), f is the Coriolis parameter, and β= df/ dy. High-resolution wind stress (Wu et al., 2011) during the same period of observation was used for the calculation.

Across the sea floor, bottom friction also gives rise to a vertical pumping similar to Ekman pumping in the upper ocean, i.e., the divergence and convergence of boundary layer transport induces a compensated vertical velocity (Pedlosky, 1979; Cushman-Roisin and MalaČiČ, 1997). Although vertical Ekman pumping is a complicated function of geostrophic vorticity, with stress-dependent eddy diffusivity (Cushman-Roisin and MalaČiČ, 1997), the following formula can be used to diagnose the bottom Ekman pumping velocity: \begin{equation} w_{\rm EKB}=\sqrt{\dfrac{v}{2f}}\xi_{\rm B} , \ \ (9)\end{equation} where v represents the mean eddy diffusivity in the bottom Ekman layer and $\xi_\rm B$ represents the bottom vorticity.

The vertical velocity (w _BFT) is induced by barotropic currents occurring over the bottom topography: \begin{equation} w_{\rm BFT}=-{U}\cdot \nabla h , \ \ (10)\end{equation} where U represents the barotropic current and h represents the water depth.

The vertical velocity integrated over the whole depth was used as the mean vertical velocity for the estimation of the spiral structure: \begin{equation} w=\dfrac{\int_{-H}^{-(H-d_{\rm EKB})}w_{\rm EKB}{\rm d}z+\int_{-H}^0w_{\rm BFT}{\rm d}z +\int_{-d_{\rm EKS}}^0w_{\rm EKS}{\rm d}z}{H+d_{\rm EKB}+d_{\rm EKS}} . \ \ (11)\end{equation} Here, d _EKB and d _EKS represent the depths of the bottom and surface Ekman layers, respectively. According to (Pacanowski and Philander, 1981), the Richardson number (Ri)-dependent vertical eddy diffusivity can be estimated from buoyancy frequency and velocity shear (Song and Yu., 2013): \begin{equation} k_v\propto \dfrac{1}{Ri^n}\propto \dfrac{(\bar {u}_z)^2}{T_z} . \ \ (12)\end{equation}

Term D on the rhs of Eq. (7) represents the dynamic contribution associated with horizontal diffusion of vortex stretching, where $\xi_z$ represents the vertical gradient of the relative vorticity. Assuming that the geostrophic eddy is axisymmetric, the contribution due to the meridional gradient of the zonal velocity -∂ u/∂ y is the same as that from the other term ∂ v/∂ x. Since only the term ∂ v/∂ x is available from in-situ observations, we make use of the working assumption for a vortex: \begin{equation} \xi=2\dfrac{\partial v}{\partial x} . \ \ (13)\end{equation}

Accordingly, the northward currents at M2, M4 and M5 can be used to calculate the vorticity at M4 from the following relation: \begin{equation} \xi=2\dfrac{v_5-\frac{1}{2}(v_2+v_4)}{2L_{M_4\to M_5}} . \ \ (14)\end{equation}

The strategy for diagnosing the observational spiral dynamics is to balance the terms in Eq. (7). Given the observations (Fig. 1; Tables 1 and 2), Terms A, B, C and E can be calculated directly; however, Term D is difficult to estimate directly from the observations. For stations M2, M5 and M6, the relative vorticity $\xi_z$ is absent, and therefore Term D is considered to be the residual. Using the relative vorticity $\xi_z$ [Eq. (14)] for M4 and the mesoscale eddy diffusivities of KA14, we can obtain a self-balanced system. In addition, using the relative vorticity $\xi_z$ and residual Term D of M4, the eddy diffusivity can be inversely derived and compared with the results of KA14 to see if it is physically reasonable.

Geostrophy is the precondition for spirals in Eq. (7). The nondimensional Rossby number Ro=U/(2Ω L _basin) for the study case is approximately 4× 10^-4 by estimating the mean velocity, Earth's rotation (Ω) and basin scale of the Yellow Sea (L _basin). The Rossby number is quite small, and the inertial effect can be ignored. The geostrophy is diagnosed in terms of the thermal wind relation using hydrological and current observations. The existence of a density front is an important factor for velocity shear in the vertical direction and, therefore, contributes to the spirals. Figure 2 shows the dynamic balance between the left and right sides of Eqs. (2) and (3). The vertical mean (e.g., $\overline{\partial u/\partial z}$) was obtained by averaging the vertical bins. Figure 2a shows that the y-axis $\overline{\partial u/\partial z}$ and x-axis $\overline{\gamma(\partial\rho/\partial y)}$ are nearly balanced at all mooring sites. Figure 2b illustrates the same feature. However, Fig. 2 also indicates that the velocity shear and density gradient are not exactly the same, with slight differences in magnitude. This is likely due to the effect of the inertial, nonlinear and frictional terms of the momentum equation. Overall, this evidence supports that geostrophic dynamics and Eq. (7) can be used to analyze geostrophic spirals. The magnitudes of $\overline{\partial v/\partial z}$ and $\overline{-\gamma(\partial\rho/\partial x)}$ at M2 (Fig. 2b) are approximately three to four times greater than those at other sites, which is due to the strong density front between the cold fresh water along the coast and the warm salty water as a result of the Yellow Sea Warm Current. The significant temperature and density fronts using our hydrologic and satellite observations have already been shown in (Lin et al., 2011); therefore, they are not shown here.

Figure 2.  The thermal wind relation diagnostic according to the observations (Fig. 1). Panel (a) shows the diagnostic of Eq. (2) for M2, M4, M5 and M6, where the x-axis represents the vertical mean meridional density gradient ($\overline{\gamma(\partial\rho/\partial y)}$; units: s^-1), and the y-axis represents the vertical mean zonal velocity shear ($\overline{\partial u/\partial z}$; units: s^-1). Panel (b) is the same as (a) but for Eq. (3).

Before presenting the velocity rotations, the Ri-dependent eddy diffusivity [Eq. (12)] is used to estimate the surface and bottom Ekman depths (i.e., the e-folding depth). Figure 3 shows the logarithmic vertical eddy diffusivity (m² s^-1) in terms of hourly velocity shear and vertical buoyancy frequency. The hourly velocity data are used to calculate the diffusivity since tidal mixing is quite significant in the Yellow Sea; through this method, tidal mixing is implicitly included. Strong mixing was observed in the surface and bottom layers. Mixing at the surface can be caused by significant wind and wave stirrings and cooling in the wintertime, and mixing at the bottom is associated with possible topographic and tidal effects. However, for the bottom layer, the maximum eddy diffusivity is not right along the bottom; instead, it is approximately 15 m above the bottom. The in-situ microstructure measurements indicate similar results (personal communication with Jiwei TIAN). The reason for such a phenomenon remains to be explored in future studies.

Figure 3.  The Ri-dependent vertical eddy diffusivity [$\lg k_v$; units: m² s^-1] in terms of buoyancy frequency and velocity shear at M2, M4, M5 and M6. The colored panels to the left represent the hourly eddy viscosities, while those to the right show the time means (×10^-3; units: m² s^-1).

The average eddy diffusivity within 20 m of the surface and bottom layers reaches 10^-3 m² s^-1 in magnitude and is listed in Table 3. The largest vertical eddy diffusivity (approximately 3.2× 10^-3 m² s^-1) is observed at the M4 station, while the largest average eddy diffusivity (1.8× 10^-3 m² s^-1) is at the M5 station. Eddy diffusivity in the Yellow Sea can be quite large, and it is likely caused by strong wintertime wind stirring and cooling, surface waves, topographic effects, and tidal dissipation. With the vertical mean viscosity and Coriolis parameter, the Ekman layer depth is $\sqrt{2v/f}$. Due to strong momentum dissipation in the Ekman layers, these layers are quite thin. In fact, the maximum surface e-folding depth is approximately 9 m at M4, and the minimum is only 2 m at M6 (Table 3). The maximum bottom e-folding depth is approximately 7 m at M5, and the minimum is 1.3 m at M4. To avoid contamination due to a strong Ekman effect in the surface and bottom layers, our spiral analysis is confined to the interior water column, which is 20 m away from the surface and bottom layers. Removed from these thin layers, the non-constant nature of vertical viscosity can be dynamically important; thus, the non-constant vertical eddy diffusivity diagnosed from the data is used to calculate Term E in Eq. (7) in the following sections.

The barotropic characteristics of currents and their dynamics at the mooring sites were discussed by (Lin et al., 2011); however, their baroclinic structures have not been fully understood yet. In general, geostrophic spirals are discussed within the framework of a steady state with a small Rossby number, and the effects of tides are filtered from the observation data by applying Lanczos low-pass-filtering (Walters and Heston, 1982), similar to (Yu et al., 2010). The cut-off frequency of the low-pass filtering is set to 1/38 of an hour. The long-term mean currents were obtained by averaging the velocity after the low-pass filtering.

Figure 4 shows the three-dimensional mean currents at mooring sites M2, M4, M5 and M6. Remarkable clockwise rotations with increasing depth can be identified from these sites, particularly at stations M2 and M5. The barotropic currents, baroclinic currents and simulated one-dimensional Ekman currents (Price and Sundermeyer, 1999) forced by the in-situ wind observations have been calculated in our previous study (Lin et al., 2011). In summary, the magnitude of the mean baroclinic current is half that of the mean barotropic current, while that of the mean Ekman current in the interior (away from the boundary layers) is no more than 10% of the baroclinic current magnitude. Although the spirals are not as smooth as the curves shown by (Schott and Stommel, 1978), which are computed by smoothing over a moderate number of stations, there are marked spirals at all four moorings. The mean spiral of the four moorings shows a net displacement to the right of the wind and a clockwise rotation with depth at a mean rate of 1.7^° m^-1. Figure 5 shows the mean currents in the interior depth and the sea surface wind vectors at the four moorings. The magnitude and direction of the current do not seem to obey the law predicted by the classic Ekman theory. At the M2 and M5 stations, there is a false appearance of an Ekman spiral. In light of the Ekman layer theory, the current should be rather weak when the angle between the wind direction and velocity vector is close to 135^° (e.g., M5). However, the currents at the M2 and M5 stations are 2 cm s^-1 and 8 cm s^-1, respectively, when the angle is approximately 135^°. Compared with a one-dimensional Ekman layer model, we concluded that the currents at this depth are not Ekman currents (Lin et al., 2011). Therefore, the observed velocity rotations cannot be explained in terms of local Ekman dynamics. Interestingly, at stations M4 and M6, the currents in the interior (i.e., 20 m away from the surface and bottom layers) flow in opposite directions compared to the surface wind (Fig. 5). Apparently, the observed velocity rotation at these locations cannot be explained in terms of local wind forcing, and we attempt to explain these long-term mean spirals using the geostrophic spiral associated with other terms in the large-scale ocean dynamics [see Eq. (7)].

Figure 4.  The structure of the absolute mean three-dimensional velocity spirals at M2, M4 and M5 in the Yellow Sea (Fig. 1). (b) As in (a) but for M6.

Figure 5.  The velocity rotations at mooring sites M2, M4, M5 and M6 (a view of the velocity profile from the surface) and velocity hodographs at the surface and bottom Ekman layers (20 m; not included). The blue arrows represent the mean wind vectors during the time of the mooring observations. The dashed orange circles with numbers represent the velocity magnitudes, and the numbers (red) denote the depths of the vectors.

Figure 6 shows the surface Ekman pumping rate [Eq. (8)] and surface cooling (Q). Based on all available observations, the contributions due to velocity rotations (Term A), mean vertical advection associated with Ekman pumping and topography-induced upwelling/downwelling (Term B), cooling effects (Term C) and vertical diffusion (Term E) at the four moorings are shown in Fig. 7. The horizontal diffusion of vortex stretching (Term D) is taken as a residual term (colored stars in Fig. 7). Figure 7 shows the dynamic relationship of Term A with Terms B, C, D and E.

Figure 6.  (a) The surface Ekman pumping rate (w _EKS; units: × 10^-7 m s^-1) in Eq. (8) calculated via high-resolution wind assimilations based on the WRF model (see the data description in section 2; Wu et al., 2011) and (b) the net surface air-sea heat fluxes (Q _NET; units: W m^-2; negative upward flux) from the OAFlux+ISCCP dataset.

In the subtropical basin interior, negative Ekman pumping dominates the dynamic balance of the clockwise beta spiral (Stommel and Schott, 1977). The vertical advection of density contributes to the maintenance of the spiral. Figure 6a shows that, in the Yellow Sea (M2, M4 and M5), surface Ekman pumping is generally negative, except for a small region that has positive values in the northern Yellow Sea, which is where M6 is located. The Ekman pumping rate is calculated from the outputs of the high-resolution models with data assimilations (Wu et al., 2011), and the stress pattern from QuikSCAT (Liu, 2002) shows the same results (not shown in this paper). The transition from negative (M2, M4 and M5) Ekman pumping to positive (M6) pumping provides a good opportunity to test the role of surface Ekman pumping. Figures 6 and 7 show that, although surface Ekman pumping has different signs, the velocity displays the same clockwise rotations with increasing depth at moorings M2, M4, M5 and M6. Thus, the surface Ekman pumping effect should not be the dominant dynamic of our observations.

The bottom Ekman pumping velocity, in terms of eddy diffusivity and geostrophic vortex (w _EKB), and the upwelling/ downwelling associated with flows over the topography (w _BFT), are commonly assumed to be small and negligible about the traditional beta spiral in the open ocean. However, in coastal seas, due to strong bottom friction and a steep topographic gradient, these terms can be larger or at least comparable to the surface Ekman pumping rate. Equation (9) is used to calculate the bottom Ekman pumping velocity with the Ri-dependent eddy diffusivity averaged in the bottom layer (Table 3) and the observed geostrophic vorticity ($\xi_\rm B$). The rectilinear moorings at a latitude of 34.4^°N (Table 1) provide a favorable opportunity to estimate the geostrophic vorticity [Eqs. (13) and (14)].

The geostrophic vorticity was obtained from the observations at only one station; therefore, these vorticity data were used to calculate bottom Ekman pumping at all moorings. Estimates of bottom Ekman pumping at the four sites are listed in Table 4. The hourly data of the barotropic current are used to estimate the vertical velocity (w _BFT), and the eventual value of the vertical velocity is obtained by averaging the hourly calculations (∑_n=1^Nw _BFT), where N represents the number of observations. Table 4 lists the number of parameters associated with topography-induced upwelling, including the barotropic velocity, current magnitude, water depth, and topographic gradient. It shows that the zonal topographic gradient (h_x) dominates the topographic effect compared to the meridional gradient (h_y) at M2, M4, and M5, while the meridional topographic gradient is comparable to the zonal gradient at M6, which is where the meridional valley exists (Fig. 1). Due to both the zonal and meridional topographic gradients, upwelling at M6 is extremely large, with a value of 6.5× 10^-2 cm s^-1. Table 4 shows that the velocity caused by the flow over the bottom is at least two orders larger than the value for surface and bottom Ekman pumping. Therefore, the mean vertical velocity is dominated by topography-induced vertical velocity [Eq. (11)].

The net contribution of vertical advection is eventually estimated by the mean current, vertical velocity [Eq. (11)] and vertical density stratification. Figure 7 demonstrates that the magnitude of the advection term (Term B) associated with such a weak stratification is quite small. The vertical advection terms contributing to rotation range from 10^-3 to 10^-2 (units: ^° m^-1), which is at least one order smaller than the observed velocity rotation (Fig. 7). Therefore, there is a very weak dynamic correlation between vertical advection and subsurface velocity rotation in our study.

The net surface heat flux (Q _NET) in the Yellow Sea is negative throughout the whole basin, as indicated in the current flux datasets (Table 5). The magnitudes of the heat fluxes at the four moorings are different, ranging from -118 W m^-2 (M2 station via MERRA) to -236 W m^-2 (M6 station via CFSR; a negative sign indicates an upward flux), which reflects the poorly constrained quality of heat fluxes entering/leaving the ocean among the different heat flux products (Song and Yu., 2013). Although the estimated surface heat flux is contaminated with great biases in the flux magnitudes, there is no question that the Yellow Sea is subject to strong heat loss in the winter. When warm water in the Yellow Sea encounters cold air from inland regions and high latitudes, a huge amount of heat is released in the forms of latent and sensible heat fluxes via evaporation and convective processes. The Term C in Eq. (7) can be expressed as \begin{equation} \dot{q}=\dfrac{\alpha Q_{\rm NET}}{hc_p} , \ \ (15)\end{equation} where α represents the thermal expansion coefficient, Q _NET represents the net surface heat flux, c_p represents the specific heat capacity under constant pressure, and h represents the mixed layer depth (MLD). However, under strong wintertime wind stirring and cold-air cooling, the MLD can reach the bottom based on a threshold value of 0.03 kg m^-3 for the density difference (de Boyer Montégut, 2004). Therefore, h is referred to as the total water depth.

The mean net heat flux (Q _NET) from the five flux products (CFSR, CORE.2, ERA-Interim, MERRA and OAFlux+ ISCCP) indicates a heat loss of -165 38 W m^-2 at M2, -168 39 W m^-2 at M4, -173 37 W m^-2 at M5, and -188 42 W m^-2 at M6. The errors after the numbers represent the standard deviations (STDs) among the five heat fluxes, which is approximately 22% of the mean. Figure 7 shows that the vertical distribution of the velocity rotations due to the cooling effect at the four mooring sites has the same direction but different magnitudes. Winter cooling in the Yellow Sea can induce an anticlockwise rotation, with magnitudes of 1.6 0.36^° m^-1, 1.9 0.45^° m^-1, 0.15 0.03^° m^-1 and 1.9 0.44^° m^-1 at the four moorings. These errors are due to uncertainties in the different heat flux datasets listed in Table 5. Cooling at M5 contributes to a relatively weaker anticlockwise rotation compared with the other stations due to deeper waters and larger absolute velocities [Eq. (7)]. The observations indicate that surface cooling can induce a relatively strong rotation, but it cannot generate a clockwise current spiral. Therefore, other physical processes are required to counterbalance negative cooling spirals and build up positive velocity spirals.

Vertical mixing (-k_vρ_zz) can also influence the velocity structure in the vertical direction [Term E of Eq. (7)], particularly in the deep ocean (Schott and Zantopp, 1980). Although vertical stratification is weak in the wintertime, the contribution of vertical diffusion to spirals cannot be ignored due to overwhelmingly large vertical mixing (Fig. 3). Vertical eddy diffusivity based on the Ri reaches 10^-3 m² s^-1 in magnitude, which is much larger than the global mean value of 10^-4 m² s^-1 (Munk and Wunsch, 1998). The depth-dependent eddy diffusivity is used to calculate Term E in Eq. (7). Figure 7 shows that the vertical diffusion (empty squares) at all stations is nearly zero. From the above analysis, there are no significant processes (e.g., vertical advection, cooling and vertical diffusion) that can induce the clockwise spirals observed in this region.

The vertical mean spiral rate diagnosed from the observations at M4 is approximately 2.7^° m^-1, while the cooling spiral rate is -1.9^° m^-1. The contributions of the mean vertical advection and vertical diffusion are close to zero (Fig. 7). So far, the spiral budget has not been balanced; thus, the only term left behind, which is the horizontal diffusion of vortex stretching [Term D on the rhs of Eq. (7)], should be explored, and it may play a dominant role in closing the rotation budget. The residual term for the horizontal diffusion of vortex stretching is 4.6^° m^-1, on average, which indicates that the contributions of cooling and vertical diffusion should be counterbalanced by vortex stretching diffusion. Using the relative vorticity at M4 and the residual from the spiral budget, we can estimate the horizontal eddy diffusivity k_h required for a balanced spiral budget. The mean meridional currents at stations M2, M4 and M5 are 0.63 cm s^-1, 2.52 cm s^-1 and 3.56 cm s^-1, respectively; thus, there is a significant cyclonic geostrophic vortex in the mean state. A positive relative vorticity can help generate and maintain the observed spirals. The magnitude of the relative vorticity is 7× 10^-7 s^-1, based on Eq. (13). Although this value is smaller than that of the planetary vorticity, the role of vortex stretching diffusion is essential in balancing the cooling effect and vertical diffusion as a result of a smaller velocity magnitude and relatively larger horizontal eddy diffusivity (KA14).

Figure 7.  The velocity spiral of Term A ( dθ/ dz) versus Ekman pumping [Term B, (γ wρ_z)/(u²+v²); empty circles], the cooling effect [Term C, $(\gamma \dot q)/(u^2+v^2)$; crosses], the horizontal diffusion of vortex stretching [Term D, $(k_h\xi_z)/(u^2+v^2)$; stars] and the vertical diffusion [Term E, -(γ k_vρ_zz)/(u²+v²); empty squares]. Note the legends in the top-left corner. The solid line denotes the polyfit of Term A versus Term D for all mooring observations. Note the different colors for the different moorings in the bottom right.

The horizontal eddy diffusivity required for balancing a spiral can be estimated as: \begin{equation} k_h=\dfrac{D_{\rm RES}}{\xi_z}(u^2+v^2) , \ \ (16)\end{equation} where D _RES represents the residual term D in Eq. (7). To balance the residual term (4.6^° m^-1), the mean horizontal eddy diffusivity required is estimated to be 2870 m² s^-1. However, the corresponding error bar is 1320 m² s^-1, due to the uncertainties in the net air-sea heat flux (Table 5). Although the error is approximately 50% of the estimated diffusivity, the diffusivity is in good agreement with that of KA14, which has values of 2752 (suppressed) and 2881 (unsuppressed) m² s^-1. This indicates that the budget for the spiral terms [Eq. (7)] can be dynamically self-balanced with a physically reasonable estimate of k_h in terms of mesoscale eddies. Using the suppressed and unsuppressed numbers and the relative vorticity $\xi_z$ at M4, Term E is shown in Fig. 7 as a cross with a circle and triangle. The calculated value of Term E for M4 is quite close to the residual estimate (red stars in Fig. 7). The thick, solid line in Fig. 7 represents the polyfit line of the residual terms at the four moorings. Surprisingly, the estimate via the suppressed estimate (a cross with a circle) is right on the line. Therefore, the main balanced terms that contribute to the observed spirals are the net effects of cooling and diffusion. It is likely that the vortex stretching diffusion processes counterbalance the cooling spiral and generate the observed spirals in the Yellow Sea.

Wintertime in-situ observations indicated the existence of current spirals in the water column beyond the surface and bottom boundary layers, with a mean rate of 1.7^° m^-1. The hydrological observations show that the spirals are geostrophically balanced according to the diagnostic of the thermal wind relation. However, the observed spirals in the Yellow Sea are different from the beta spiral in the open ocean (Stommel and Schott, 1977) and the pure cooling spirals along the western boundary extensions (Spall, 1992). The spirals in the Yellow Sea are balanced by both the cooling effect and diffusion processes. Dynamically, spirals are regulated by a counterbalance between surface cooling and the horizontal diffusion of vortex stretching. Intensified surface cooling and wind stirring in winter give rise to weak stratification. As a result, the vertical advection term is quite weak. Horizontal eddy diffusivity can be diagnosed by regarding the horizontal diffusion of vortex stretching as a residual and estimating the relative vorticity based on rectilinear moorings. Eventually, the eddy diffusivity, 2870 1320 m² s^-1, is required to balance the dynamic system. The residual viscosity number agrees well with that of KA14.

The dynamics of the spirals, such as the classical beta spiral and cooling spiral, are controlled and generated by changes in the vertical density. The classic beta spiral and pure cooling spiral in the western boundary extensions are dominated by the vertical advection of density and cooling-induced density modifications, respectively. However, in this study, the density changes are closely related to the cooling and diffusion processes, when the vertical advection is weak. Dynamically, the horizontal diffusion of vortex stretching is the potential factor to determine the spiral by controlling the vertical density.

In this paper, we only consider the time-mean velocity spirals. In fact, the velocity vectors rotate at all times during the observation period. Although clockwise rotations are dominant throughout the whole period, there are occasional anticlockwise rotating spirals. The dynamic details on how surface cooling associated with cold-air outbreaks and eddies modify the rotating direction of velocity spirals are left for future studies. Additional long-term observations in this region are currently underway, and we hope that these new observations will reveal more dynamic details of spirals in the coastal ocean.