Journal of Oceanology and Limnology   2021, Vol. 39 issue(1): 64-78     PDF
Institute of Oceanology, Chinese Academy of Sciences

Article Information

BI Congcong, YAO Zhigang, BAO Xianwen, ZHANG Cong, DING Yang, LIU Xihui, GUO Junru
The sensitivity of numerical simulation to vertical mixing parameterization schemes: a case study for the Yellow Sea Cold Water Mass
Journal of Oceanology and Limnology, 39(1): 64-78

Article History

Received Oct. 10, 2019
accepted in principle Nov. 12, 2019
accepted for publication Dec. 3, 2019
The sensitivity of numerical simulation to vertical mixing parameterization schemes: a case study for the Yellow Sea Cold Water Mass
Congcong BI1,2, Zhigang YAO1,2, Xianwen BAO1,2, Cong ZHANG3,4, Yang DING2, Xihui LIU5, Junru GUO6     
1 College of Oceanic and Atmospheric Sciences, Ocean University of China, Qingdao 266100, China;
2 Physical Oceanography Laboratory/CIMST, Ocean University of China and Qingdao National Laboratory for Marine Science and Technology, Qingdao 266100, China;
3 The First Institute of Oceanography, Ministry of Natural Resources, Qingdao 266061, China;
4 Laboratory for Regional Oceanography and Numerical Modeling, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266000, China;
5 China National Environmental Monitoring Centre, Beijing 100012, China;
6 Operational Oceanography Institution, Dalian Ocean University, Dalian 116023, China
Abstract: The vertical mixing parameterization scheme, by providing the effects of some explicitly missed physical processes and more importantly closing the energy budgets, is a critical model component and therefore imposes significant impacts on model performance. The Yellow Sea Cold Water Mass (YSCWM), as the most striking and unique phenomenon in the Yellow Sea during summer, is dramatically affected by vertical mixing process during its each stage and therefore seriously sensitive to the proper choice of parameterization scheme. In this paper, a hindcast of YSCWM in winter of 2006 was implemented by using the Regional Ocean Modeling System (ROMS). Three popular parameterization schemes, including the level 2.5 Mellor-Yamada closure (M-Y 2.5), Generic Length Scale closure (GLS) and K-Profile Parameterization (KPP), were tested and compared with each other by conducting a series of sensitivity model experiments. The influence of different parameterization schemes on modeling the YSCWM was then carefully examined and assessed based on these model experiments. Although reasonable thermal structure and its seasonal variation were well reproduced by all schemes, considerable differences could still be found among all experiments. A warmer and spatially smaller simulation of YSCWM, with very strong thermocline, appeared in M-Y 2.5 experiment, while a spatially larger YSCWM with shallow mixed layer was found in GLS and KPP schemes. Among all the experiments, the discrepancy, indicated by core temperature, appeared since spring, and grew gradually by the end of November. Additional experiments also confirmed that the increase of background diffusivity could effectively weaken the YSCWM, in either strength or coverage. Surface wave, another contributor in upper layer, was found responsible for the shrinkage of YSCWM coverage. The treatment of wave effect as an additional turbulence production term in prognostic equation was shown to be more superior to the strategy of directly increasing diffusivity for a coastal region.
Keywords: Yellow Sea Cold Water Mass (YSCWM)    vertical mixing parameterizations    thermocline    background diffusivity    surface wave induced mixing    

The Yellow Sea Cold Water Mass (YSCWM), as the most striking phenomenon in the summer Yellow Sea (YS), occupies the central trough of the YS and distinguishes itself from surrounding waters by much lower temperature. The YSCWM, by its vast body and distinct properties, dominates the summer hydrographical fields as well as corresponding circulation patterns (Ho et al., 1959; Guan, 1963; Zhang and Yang, 1996). The YSCWM was convinced to be essentially the remnant of cold water from previous winter under joint effects of surface warming, vertical mixing, and other factors such as bathymetry, and it experiences apparent seasonal cycle with its peak in summer (Yu et al., 2006; Yao et al., 2012). Vertical mixing processes, indicated by thermocline as well as front structures, may play critical role in development of the YSCWM (Ho et al., 1959; Guan, 1963; Zhao, 1985, 1987; Yu et al., 2006; Yao et al., 2012). The significance of YSCWM is not limited within only hydrographical aspects, but also proved in biochemical and other fields in the YS. For instance, constant accumulations of nutrients, such as inorganic nitrogen, phosphate, and silicate, are found in bottom cold waters (Wang, 2000), while peaks of dissolved oxygen, pH, and chlorophyll-a in vertical usually exist at the base of thermocline (Xin et al., 2015). Fish activity also shows a close relationship with the YSCWM (Cho, 1982; Zhang et al., 1983; Liu and Ning, 2011).

However, our knowledge of YSCWM remains limited by rare observation and considerable mismatches in numerical solution, either overestimation (Ren and Zhan, 2005; Lü et al., 2010; Han and Yuan, 2014; Li, 2016), characterized by much colder and larger YSCWM, or underestimation (Song et al., 2009) of the YSCWM. Discrepancies such as in spatial coverage and core location of the YSCWM in previous studies (Moon et al., 2009; Yang et al., 2014; Yu and Guo, 2018; Zhu et al., 2018) further confirm challenges in YSCWM modelling with good accuracy, and one of the major difficulties comes from parameterization of vertical mixing processes, which is responsible for vertical transfers of heat, salt as well as momentum flux and thus plays essential role in development of the YSCWM (Zhao, 1985, 1989; Su and Su, 1996). The vertical mixing schemes is essential for model energy budgets, and a variety of schemes, such as K-Profile Parameterization (KPP) (Large et al., 1994) and two-equation model (Mellor and Yamada, 1982; Umlauf and Burchard, 2003), have been developed based on different theoretical and physical bases, thus particularly applicable to specific regions or processes (Wijesekera et al., 2003; Durski et al., 2004; Li et al., 2005). Furthermore, additional contributors of surface waves (Qiao et al., 2004, 2010; Ghantous and Babanin, 2014), Langmuir circulation (Smyth et al., 2002; Noh et al., 2016) and near-inertial waves (Jochum et al., 2013) are recently proposed for the improvement of original schemes. Although reasonable result is obtained in the open ocean, the performance of vertical mixing scheme may be considerably compromised or even dubious in coastal region, especially in cases like summer YS, where water column ranges from well mixed to highly stratified over a few kilometers and surface/bottom boundary layers could interact with each other. In addition, coastal model is rarely with the accuracy necessary for vertical mixing estimates when directly compared with small-scale turbulence measurements in the field. Situation may be even worse in summer YS considering the facts that: 1) stratification ranges over several orders of magnitude; 2) surface/bottom turbulent boundary layers may interact with each other; 3) horizontal density gradients interact with vertical mixing in frontal zones; 4) stratification may interact with bathymetric variation. The dynamics associated with YSCWM is complex, but vertical mixing is definitely one of the most important factors, so it is necessary to carefully examine and assess the influence of vertical mixing schemes on the YSCWM simulation, which, however, has not been made yet so far.

The detailed understanding of how vertical mixing schemes give different solutions could facilitate our knowledge, and help to determine what processes or mechanisms could be critical to the development of YSCWM. The objective of this paper, by clarifying the differences of simulated YSCWM due to mixing parameterizations and their possible causes, aims to reveal how future theoretical, laboratory, highresolution numerical, or field work can most profitably advance our understanding of the YSCWM as well as vertical mixing processes, and lead to improvement in simulating the YSCWM.

2 METHODOLOGY 2.1 Brief description of the model and vertical mixing parameterizations

The three-dimensional hydrodynamic model in this study was adapted from Regional Ocean Modeling System (ROMS) (Shchepetkin and McWilliams, 2003, 2005), which is a non-linear, baroclinic circulation model discrete from primitive NavierStokes equations, and provides more accurate numerical algorithms than its predecessor. ROMS contains a variety of methods for setting vertical mixing coefficients, covering three popular turbulent mixing schemes of KPP, the level 2.5 Mellor-Yamada closure (M-Y 2.5) (Mellor and Yamada, 1982) and the Generic Length Scale (GLS) scheme (Umlauf and Burchard, 2003), which serves as the basis of our sensitivity experiments and assessment.

The KPP scheme of Large et al. (1994) matches separate parameterizations for vertical mixing of the surface boundary layer and the ocean interior. A formulation based on boundary layer similarity theory is applied above a calculated boundary layer depth, at the base of which the parameterization is matched with mixing formulations to account for local shear, double diffusive mixing, and internal wave effects. The KPP scheme was adapted for the shallow sea application by appending a bottom boundary layer formulation to its original version, which was developed for the open ocean. The GLS scheme is a generic two-equation turbulence model, where the vertical mixing coefficient comes from the turbulent kinetic energy (TKE), a length scale (l) and the stability function. TKE and l can be solved from two prognostic equations which explicitly describe the energetics of mixing, with one for the evolution of turbulent kinetic energy (TKE), and the other for a length-scale determining parameter (ψ), such that the GLS scheme could behave like a variety of traditional two-equation turbulence closures, such as the k-kl, k-ω, gen scheme and so on, by properly setting parameters. The M-Y 2.5 model of ROMS, modified following Galperin et al. (1988), is essentially identical to the k-kl of GLS scheme, except differences in several details (Warner et al., 2005).

Ocean surface waves may significantly contribute to the vertical mixing in upper ocean, and hence the thermal structure of upper layers could be improved by introducing wave-induced mixing effects, especially for sea surface temperature and mixed layer depth (Craig and Banner, 1994; Mellor, 2003; Qiao et al., 2004; Babanin et al., 2009). Therefore, two additional experiments were conducted to examine the effects of surface wave, differing in how non-breaking wave-induced turbulence was handled in model. The first experiment is to modify Reynolds stress (flux) by an additional term of vertical kinematic viscosity (diffusivity) Bv, following Qiao et al. (2004). For a monochromatic wave in finite water depth, the Bv coefficient can be expressed in a simpler form (Qiao et al., 2010),


where α is a constant; A, k, and ω denote the surface wave amplitude, number, and angular frequency; h and z represent the water depth and vertical coordinate depth, respectively. Bv from Eq.1 is directly added to original mixing coefficients. The other experiment considers wave dissipation as an extra source of turbulence production, and makes an additional term PW due to wave orbital motion in turbulence closure equations, following Ghantous and Babanin (2014),


where b1 is a dimensionless parameter, H is surface wave height, and k, ωp, z are similar to those in Eq.1.

2.2 Model configuration and numerical experiments

Our model domain covers the East China Sea region (117.40°E–134.25°E, 23.00°N–41.20°N) with a spatial resolution of 3 minutes in horizontal and 30 layers in vertical. The bathymetry comes from three parts: measurements of bathymetry close to the coast of China mainland, a bathymetry database from Sung Kyun Kwan University (Korea) for other region within the first island chain, and the ETOPO1 dataset for the remaining. Bathymetry is properly smoothed for model stability, with its range set to 7 m and 6 000 m (Fig. 1a). Tidal contribution is introduced along the open boundary with 8 tidal constituents (M2, S2, K1, O1, N2, K2, P1, Q1) from TPXO7-atlas solution of OSU Tidal Inversion Software (OTIS, The atmospheric forcing is exerted in form of bulk formulas, driven by 3-6 hourly ERA-Interim Reanalysis dataset of winds (10 m), air temperature (2 m), sea level pressure, and shortwave and longwave radiation. The model is oneway nested within 1/12 degree global HYCOM/ NCODA solution in terms of initial and boundary inputs. Surface wave parameters are provided by 3 hourly WAWEWATCH III solution from National Ocean and Atmosphere Administration (NOAA). Discharges of four major rivers, including the Changjiang (Yangtze) River, Huanghe (Yellow) River, Liaohe and Yalu Rivers, are included to provide buoyancy forcing, with monthly discharge extracted from Chinese River Sediment Bulletin (

Fig.1 Model domain and topography (in m) (a), and locations of hydrographic survey sites (black dots) and observed bottom temperature (color image, with the 8 ℃ and 10 ℃ isotherms highlighted in bold black lines) in the NYS in summer of 2006 (b) In (a), the NYS region is outlined by a red rectangular; in (b), two selected sections (the D-C and 37°N sections) are defined as red solid lines connected the sites; grey contours denote isobaths in meters; and the black pentagram marks geographical location of climatological cold core for YSCWM in NYS.

Model was first integrated over 3 years for spin-up, from January 1, 2003 to December 31, 2005 using M-Y 2.5 scheme, and then, sensitivity experiments, only differing in vertical mixing parameterization, were implemented for another year of 2006 for further analysis. Details of all experiments are listed in Table 1. Among all, Exp1–Exp5 adopt background diffusivity of 0.01 cm2/s with different schemes, and a background diffusivity of 0.1 cm2/s with M-Y 2.5 scheme for Exp6, with two wave-induced mixing experiments in Exp7 and Exp8. The setting for the other parameters follows ROMS's suggestion (as shown in Appendix A).

Table 1 Details of experiment settings

As the YSCWM is the focus of our interest, the associated indexes of YSCWM, such as spatial distribution, strength, temporal variation of cold waters as well as accompanied thermocline and frontal structure, will be used for further analysis of parameterization assessment. A comprehensive survey in summer of 2006 (07/23–08/06), covering the northern Yellow Sea (NYS) region west of 124°E and north of 37°N, provides ground truth from highresolution Conductivity, Temperature and Depth (CTD) profiles of 120 stations, and hence our analysis and comparison will be made mainly focused on the NYS region.

3 RESULT 3.1 Thermal structure of the YSCWM in summer

In summer of 2006, the main body of northern YSCWM, bounded by isotherm of 8 ℃, mainly occupies the central trough of NYS west of 123.5°E and deeper than 40 m, accounting for ~25% of the total NYS area, in an elongated oval shape running northwestward (Fig. 1b). The isotherms roughly follow the isobaths, indicating the importance of topography in development of the YSCWM, when some discrepancies are also noticeable in eastern NYS. The core of YSCWM is even colder than 6 ℃, leading to significant difference highlighted by sharp temperature gradient across frontal zones, among which the most striking features can be identified to south of Liaodong Peninsular (hereafter the northern front) and north of Shandong Peninsular (the southern front). Along longitude of 122.12°E, across the cold core of climatological YSCWM in NYS, a more sharp gradient structure exists in northern front other than its southern counterpart.

Results from sensitivity experiments are presented in Fig. 2 for bottom temperature, zooming in Bohai Sea and the YS region, and derived from model solution over period of cruise in 2006. The YSCWM, in all experiments, is reasonably reproduced with acceptable accuracy, confirmed by the resemblance of modelled result to its observed counterpart. However, discrepancies, among all experiments, are also moderately apparent for the YSCWM, indexed by its spatial coverage, shape as well as location and temperature of cold core. An overestimation of YSCWM is obtained in GLS (Exp2–4) and KPP (Exp5) schemes experiments, highlighted by increased spatial coverage accounting for ~1/3 of total NYS and core temperature of less than 6 ℃, when k-kl and k-ω experiments (Exp2–3) present a cold core structure of far too eastward shift relative to its normal location. In contrast, an underestimation of YSCWM appears in the M-Y 2.5 closure experiment (Exp1) characterized by decreasing spatial coverage and increased core temperature of ~7 ℃, as a result of which a double-core structure instead occupies the central trough of northern and southern YS, respectively. By enhancing background diffusivity, the decrease of YSCWM is even more apparent, as shown in Exp6, and its coverage, bounded by isotherm of 8 ℃, accounts for ~11% of total NYS, in contrast to the ~25% of observational counterpart. The effects of wave-induced mixing (Exp7–8) are similar in weakening the YSCWM, among which direct modification of vertical mixing coefficients by Bv (Exp7) is responsible for the smallest coverage of YSCWM among all experiments. It should be admitted that some model/observation discrepancies still exist in our modelling results, for example, the running of isotherms slightly differs from its observational part along eastern boundary of YSCWM, by paralleling to the isobaths. In addition, even larger mismatch can also be expected for bottom temperature in shallow regions like Bohai Sea and some coastal area, possibly due to their smaller heat contents.

Fig.2 Bottom temperature (℃) from all experiments in Bohai Sea and the YS in summer of 2006 (time averaged results over the cruise period)

Further inspection is also made along two popular historical transects, i.e. Dalian-Chengshantou (D-C) and 37°N sections as illustrated in Fig. 1b. As shown in Figs. 3a & 4a, a giant thermocline, separating cold YSCWM in bottom and warm water at surface in vertical, and sharp fronts, isolating cold YSCWM from surrounding coastal waters in horizontal, together form a cap-shaped structure that dominates the entire trough region and is consistent with previous findings (Miao et al., 1991). The surface mixed layer reaches a depth of ~10 m before experiencing dramatic changes within the thermocline (10–25 m), below which a very minor thermal variation exists inside the YSCWM. It is noteworthy that considerable fluctuations can be found for the thermocline in either its depth or intensity, indicating many other factors, such as topographic features, advection, and internal waves, may be responsible for temporal and spatial variations of YSCWM. The intensity of thermocline associated with YSCWM, quantified by vertical thermal gradient, can reach an average of ~1 ℃/m along two sections, which reflects the effects of turbulent mixing and other processes.

Fig.3 Comparison of the observed (a) and simulated (b–i) temperature (℃) from all experiments along the D-C section in summer of 2006 Locations of observational data are marked by black dots in (a).
Fig.4 Same as Fig. 3, but for 37°N section

Along D-C section (Fig. 3), the dome-shaped isotherms, forming jointly by thermocline in vertical and fronts in horizontal, are clearly seen from all experiments, which confirm the ability of available mixing schemes. However, the fluctuations of small scales are absent from all numerical solutions, indicating that some processes or mechanisms may be missed in our parameterization schemes or even the numeric model itself. For example, under hydrostatic approximation adopted by ROMS as well as most circulation models, interval wave cannot be well resolved, which may be responsible for absent smallscale fluctuations as shown in Fig. 3 and beyond the scope of this paper. Another index to be examined is the thermocline, the transitional zone where YSCWM competes with surface warmer water by vertical mixing processes. The thermocline in M-Y 2.5 (Exp1) lies between 10–15 m, the thickness of which is significantly smaller than the observed with the largest thermal gradient among all experiments, and upper bound of YSCWM could reach up to ~15 m. The thermal gradient is reduced in GLS (Exp2–4) and KPP (Exp5) scenarios, supported by expanded transitional zone at 5–15 m as well as a shallower upper mixed layer. Despite difference in upper reach of YSCWM, which falls to the depth of 15–20 m in k-kl and k-ω (Exp2–3), the GLS and KPP schemes (Exp2–5) lead to a larger extent in horizontal, relative to the ground truth. By increasing background diffusivity (Exp6), the transition from warm surface to cold YSCWM becomes smoother, supported by a considerably weakened thermal gradient and decreasing trend of YSCWM top bounds. Additional wave-induced mixing apparently extends surface mixed layer to a depth below 10 m, as revealed in both experiments (Exp7–8), and dramatically weakens YSCWM, in either spatial extent or its strength. Vertical thermal gradient experiences minor change in Exp8, in contrast to considerable decrease in Exp7. The results, along parallel of 37°N, also show similar trend as that along D-C section, although slight difference may exist (Fig. 4), and therefore, the repetitious detail for 37°N section is not given here. Overall, the experiment of enhanced background diffusivity (Exp6) gives the best agreement of spatial extent, bounded by isotherm of 10 ℃, and location of cold core to the observed along D-C section, but fails in a higher core temperature. As for section of 37°N, Exp5 stands out in the shape of YSCWM and thermal gradient when Exp8 emerges by its tidal front structure.

3.2 Seasonal variation of the YSCWM

Seasonal cycle of the YSCWM is further examined to illustrate how the temporal variation of YSCWM is modulated by parameterization schemes. The cold core of climatological YSCWM in NYS is taken as our analysis site, roughly at 122.12°E, 38.22°N (Li et al., 2015), with results presented in Fig. 5. Thermocline structure appears in mid-April, and grows gradually in May before rapid increase of thermal gradient due to solar heating in summer, when peak temperature of more than 26 ℃ at surface is reached in mid-August. The trend of surface temperature is reversed in autumn, dropping to ~16 ℃ in early November. In contrast to the surface, characterized by dramatic fluctuation, temporal variation is relatively small for lower depth. A slow increasing trend dominates the lower levels until November when it is accelerated. The divergence of numerical integration appears since late April, caused by different mixing schemes. In M-Y 2.5 case (Exp1), the lower cold water warms rapidly until mid-June, and its core temperature rises up to 8 ℃ during early September, after which the increase slows down, and it decays to 10 ℃ finally in early November when a well-mixed water column gradually forms later. The layer thickness of YSCWM basically stabilizes during summer season, with its upper reach up to 25 m mostly. The k-ω (Exp3) and k-kl (Exp2) cases demonstrate YSCWM temporal evolution of different temperature. The lower layer slowly warms until July, and with the deepening of surface mixed layer as well as the great reduction of thermocline since August, its temperature eventually rises up to 8 ℃ by the end of September and 10 ℃ by early November, respectively. The gen (Exp4) and KPP (Exp5) cases show a similar temporal curve of temperature of lower layers, which is much colder than that in other experiments and remains less than 8 ℃ until November. By enhancing background diffusivity (Exp6), the YSCWM is significantly eroded, confirmed by reduced thermal gradient and one month shorter time window of YSCWM, especially after mid-June when the decay accelerates. The wave-induced mixing facilitates the downward transfer of momentum and heat flux, revealed by deepening of surface mixed layer. Reduced thermal gradient as well as more rapid warming of bottom cold water could also be found in Exp7. Conversely, by modification of turbulence production term in Exp8, only slight warming trend is found in YSCWM and does not impose much influence on its time window.

Fig.5 Temporal evolution of temperature (℃) in 2006 at site of 122.12°E, 38.22°N, location of climatological cold core for YSCWM in NYS, from all experiments

Overall, the divergence of model results, among all experiments, appears since spring and grows continuously until late November, a transitional period before two-layer structure is finally destroyed in December, due to massive heat losses from the surface. Hydrographic condition of this transitional stage is thus examined, with monthly bottom temperature of November in 2006 shown in Fig. 6. The difference is found to be particularly striking in regions deeper than 30 m. The GLS (Exp2–4) and KPP (Exp5) cases are highlighted by spatial coverage of cold water, which, bounded by isotherm of 10 ℃, extends from 34°N to 39°N meridionally. Especially in gen (Exp4) and KPP (Exp5) cases, abnormally cold waters of less than 8 ℃ occupy the bottom trough of YS, which is inconsistent with our previous knowledge (Editorial Board for Marine Atlas, 1992; Yu et al., 2006) and suggests vertical mixing by these two schemes is too weak. Oppositely, extremely small cold water area, in a box of roughly 1°×1°, is presented in wave-induced mixing experiment by Bv coefficient (Exp7), indicating an over-mixing case in vertical, when the solution from Exp8 is much better by expression of wave effects in turbulence production form.

Fig.6 Bottom temperature (℃) from all experiments in Bohai Sea and the YS in November of 2006 (monthly averaged results)
3.3 Statistical characteristics and quantitative assessment

Statistical indexes, defined by characteristics associated with YSCWM, are extracted and analyzed for the YSCWM of summer in 2006, to provide quantitative assessment, consisting of core temperature, spatial coverage percentage and upper reach of the YSCWM, duration, intensity, depth and thickness of the thermocline as well as location and intensity of the fronts, as listed in Table 2. Parameters for the thermocline are extracted at core location of northern YSCWM, and those for fronts are taken along longitude of 122.12°E. The M-Y 2.5 scheme is defined by a moderately weakened cold core and too sharp thermal gradient, when GLS and KPP schemes present a YSCWM with too large spatial coverage in either horizontal or vertical, and a shallower mixed layer. The thermal gradient (intensity) and thickness of thermocline varies significantly, by increasing background diffusivity, leading to a reduced YSCWM by area percentage and top boundary depth. The effects of wave-enhanced mixing are illustrated by deepening of surface mixed layer and reduced coverage of cold water, when the strategy of Bv leads to severely weakened stratification.

Table 2 Comparison of the observed and simulated statistical characteristics of northern YSCWM in summer of 2006

Responses of numerical model to vertical mixing parameterization schemes could be very complex, and some indexes may be improved when the others deteriorate, as a result, it is difficult to assess each scheme by its total performance. Here a Skill Score (SS) parameter (Murphy, 1988) is adopted as an indicator of model performance,


where X, F are observational data and corresponding simulation results respectively, n is sample number, and X is the mean value of X. According to its definition, a larger SS corresponds to a smaller model/ data error, with a maximum of 1 for perfect agreement. SS parameter is calculated for all experiments based on cruise data in 2006 (listed in Table 3). The M-Y 2.5 scheme provides a good model/observation agreement with a SS of 0.81, which could be optimized to 0.83 by increasing background diffusivity. The GLS of k-kl and k-ω, as well as the KPP scheme follows with an acceptable SS of 0.80, 0.79, and 0.76, when a depressing value of 0.71 comes from the GLS of gen. In addition, model performance confirms the superior of treatment of wave-induced mixing as an extra source of turbulence production than direct increase of mixing coefficient, supported by SS dropping from 0.75 to 0.58. However, the wave-induced mixing effects, expressed by above two strategies, do not significantly improve model performance as expected and possible explanations will be discussed later.

Table 3 Skill Scores from all experiments in summer of 2006

Despite of various parameterization schemes, the basic idea is to modulate turbulent mixing coefficient, by which effects of vertical mixing are introduced in numerical model. Therefore, it is necessary to examine how turbulent mixing coefficient responds to different mixing schemes. The vertical profiles of turbulent mixing coefficient, averaged over period of 2006 cruise, are shown for all experiments at site of 122.12°E, 38.22°N in Figs. 7 & 8.

Fig.7 Profiles of diffusivity (a, c, e) and vertical density gradient (b, d, f) from Exp1 to Exp5 at site of 122.12°E, 38.22°N in summer of 2006, with entire profile in the 1st row, upper mixed layer in the 2nd row and thermocline in the 3rd row, respectively
Fig.8 Same as Fig. 7, but for comparison of Exp1 and Exp6–8

The coefficient experiences significant variation with depth, which also largely depends on its vertical mixing scheme (Fig. 7). In upper mixed layer (≤4 m), the M-Y 2.5 and k-kl schemes are highlighted by strong vertical mixing, when other schemes like k-ω and KPP lead to a weak condition, as a result of which a shallow upper mixed layer appears. Moving downward to the thermocline of approximate 4–30 m depth, where stratification is strong and vertical exchanges are efficiently suppressed, extremely weak mixing condition could dominate the entire zone for all experiments, revealed by dramatically reduced diffusivity of less than 0.2 cm2/s. Among all, mixing coefficient of M-Y 2.5 and gen schemes even is reduced to background diffusivity of 0.01 cm2/s, corresponding to sharp thermal gradient and thinner thermocline, while that of k-ω and k-kl is relatively big, leading to weak thermocline with increased thickness. It is interesting that the coefficient minimum by KPP scheme appears at 10–20 m, differing from pycnocline in depth of 3–15 m. Below thermocline, tidal mixing process dominates the remaining water column, where intensive stirring by bottom friction makes water well mixed and density difference almost vanishes. Among all, the diffusivity from gen scheme stands out, followed by that of k-ω, k-kl and M-Y 2.5, when the minimum occurs in KPP case.

Experiment results for sensitivity to background diffusivity and wave effects are shown in Fig. 8, which are based on M-Y 2.5 scheme. By increasing the background, lower limit is largely lifted for diffusivity within pycnocline, from 0.01 cm2/s to 0.1 cm2/s, where vertical mixing is supposed to be greatly suppressed by buoyancy, but it contributes little outside of the pycnocline. Despite limited level of mixing variation, the effects of increasing background diffusivity could gradually accumulate with time, which should be responsible for density gradient spiking nearby the upper bound (~10 m) but gradual decrease toward the lower bound (10–25 m) of thermocline. The contribution of waves, referring to Bv and PW formulas (Eq.1 and Eq.2), is basically confined within the upper layers, and decays quickly with depth, revealed by model-calculated Bv profile. In strategy of Bv, the suppression by stratification to mixing is not considered in either complete form of wave number spectrum (Qiao et al., 2004) or its simplification for monochromatic wave only (Qiao et al., 2010), thus over mixing can be expected for coastal region where thermocline is strong and shallow in summer, as shown in Exp7. The other case, by strategy of PW, could be well constrained by stability function and other factors, and hence provides a reasonable level for mixing. From this point of view, the second strategy of PW is superior for coastal region than the first, which still needs more efforts before it is applicable to shallow regions.

In summary, the YSCWM could be considerably sensitive to mixing parameterization schemes, so carefully testing as well as optimization of turbulence model may be needed for the simulation of YSCWM. In addition, other factors, such as advection, surface forcing frequency, and absorption profile of shortwave radiation, may also contribute to the simulation, which still needs more efforts in future.


The study focuses on finding points of difference among various vertical mixing schemes, by hindcasting the YSCWM in summer of 2006 using ROMS model. Three schemes, including the M-Y 2.5 closure, KPP as well as GLS schemes, are carefully examined and assessed, motivated by the fact that all schemes, as originally formulated, inaccurately represent mixing in nearshore region where surface and bottom boundary layers, topography may interact strongly. The results are as follows:

1) The YSCWM structure could be reasonably reproduced by the M-Y 2.5, GLS and KPP mixing schemes. The M-Y 2.5 scheme is characterized by a moderately weak YSCWM, with sharp thermal gradient above, when GLS and KPP are responsible for a too strong YSCWM, with a shallower thermocline;

2) The major difference of temperature evolution, among all experiments, occurs in the lower mixed layer, and it appears since spring, growing gradually until late November;

3) The thermal gradient of thermocline could be largely eroded, by a larger background diffusivity, when the YSCWM can be significantly weakened, in either coverage or intensity;

4) The addition of wave induced mixing into numerical model, by deepening the upper mixed layer and weakening the lower YSCWM, cannot improve the model performance as expected, at least for our YSCWM case. However, the strategy of PW, by enhancing turbulence production, is superior in coastal region than that of Bv, by direct increase of diffusivity, which may overestimate vertical mixing near seasonal pycnocline by ignoring the suppression effect of stratification to turbulence, so that still needs more caution and efforts before it is applicable to shallow regions.

It would be useful if the differences found by model performance could suggest a superior to another, and facilitate further studies. However, a conclusive statement cannot be made for coastal applications, or even the YSCWM simulation, considering the limits of each schemes and complexity of vertical mixing in coastal region. More efforts are still needed, either to improve the YSCWM simulation, or optimizing our mixing parameterization. Hopefully the points of difference noted in this study could help direct such future efforts.


According to "Regulations on the Scope of State Secrets in Oceanic Work", observational data used in this study belong to the category of high confidentiality and cannot be shared publicly over a long period of time. All analysis results mentioned above have been shown in the form of figures and tables in this paper. Data are however available from the corresponding author ZGY upon reasonable request.


We thank the Data and Simulation Center of the Physical Oceanography Laboratory, Ocean University of China, for providing super computers for all the model simulations. We thank the National Science & Technology Infrastructure-the Dalian branch of National Marine Scientific Data Center ( for providing data support. We also thank the National Marine Data and Information Service of China, Sung Kyun Kwan University (Korea), OSU Tidal Data Inversion, European Centre for Medium-Range Weather Forecasts (ECMWF), HYCOM Consortium, National Ocean and Atmosphere Administration (NOAA), and Bureau of Hydrology, Ministry of Water Resources of China for providing valuable data.

APPENDIX Appendix A Parameter values for vertical mixing used in this paper

The mixing associated parameters, except for background diffusivity, are set as suggested in ROMS. Specifically, background vertical mixing coefficient for momentum, TKE and ψ are set to 1.0×10-5, 5.0×10-6, 5.0×10-6 m2/s, respectively; choices of various parameters in GLS scheme are listed in Table A1, where p, m, n and cμ0 are employed to establish the generic parameter ψ (ψ=(cμ0)p(TKE)mln), TKEmin and ψmin are minimum value of TKE and ψ, c1, c2, c3 are coefficients in the equation for ψ, and σ is constant Schmidt number.

Table A1 Suggested parameter values for various turbulence closure models from GLS scheme
Babanin A V, Ganopolski A, Phillips W R C. 2009. Waveinduced upper-ocean mixing in a climate model of intermediate complexity. Ocean Modelling, 29(3): 189-197. DOI:10.1016/j.ocemod.2009.04.003
Cho K D. 1982. On the influence of the Yellow Sea Bottom Cold Water on the demersal fishing grounds. Journal of the Korean Society of Fisheries and Ocean Technology, 18(1): 25-33.
Craig P D, Banner M L. 1994. Modeling wave-enhanced turbulence in the ocean surface layer. Journal of Physical Oceanography, 24(12): 2 546-2 559. DOI:10.1175/1520-0485(1994)024<2546:MWETIT>2.0.CO;2
Durski S M, Glenn S M, Haidvogel D B. 2004. Vertical mixing schemes in the coastal ocean:comparison of the level 2.5 Mellor-Yamada scheme with an enhanced version of the K profile parameterization.. Journal of Geophysical Research, 109(C1): C01015.
Editorial Board for Marine Atlas. 1992. Marine Atlas of Bohai Sea, Yellow Sea, East China Sea: Hydrology. China Ocean Press, Beijing, China. p. 89. (in Chinese)
Galperin B, Kantha L H, Hassid S, Rosati A. 1988. A quasiequilibrium turbulent energy model for geophysical flows. Journal of the Atmospheric Sciences, 45(1): 55-62. DOI:10.1175/1520-0469(1988)045<0055:AQETEM>2.0.CO;2
Ghantous M, Babanin A V. 2014. One-dimensional modelling of upper ocean mixing by turbulence due to wave orbital motion. Nonlinear Processes in Geophysics, 21(1): 325-338. DOI:10.5194/npg-21-325-2014
Guan B X. 1963. A preliminary study of the temperature variations and the characteristics of the circulation of the cold water mass of the Yellow Sea. Oceanologia et Limnologia Sinica, 5(4): 255-284. (in Chinese with English abstract)
Han L, Yuan Y L. 2014. An ocean circulation model based on Eulerian forward-backward difference scheme and threedimensional, primitive equations and its application in regional simulations. Journal of Hydrodynamics, 26(1): 37-49. DOI:10.1016/S1001-6058(14)60005-6
Ho C P, Wang Y X, Lei Z Y, Xu S. 1959. A prelimenary study of the formation of Yellow Sea Cold Mass and its properties. Oceanologia et Limnologia Sinica, 2(1): 11-15. (in Chinese with English abstract)
Jochum M, Briegleb B P, Danabasoglu G, Large W G, Norton N J, Jayne S R, Bryan F O. 2013. The impact of oceanic near-inertial waves on climate. Journal of Climate, 26(9): 2 833-2 844. DOI:10.1175/JCLI-D-12-00181.1
Large W G, McWilliams J C, Doney S C. 1994. Oceanic vertical mixing:a review and a model with a nonlocal boundary layer parameterization. Reviews of Geophysics, 32(4): 363-403. DOI:10.1029/94RG01872
Li A, Yu F, Diao X Y, Si G C. 2015. Interannual variability of temperature of the northern Yellow Sea Cold Water Mass. Acta Oceanologica Sinica, 37(1): 30-42. (in Chinese with English abstract)
Li A. 2016. The Study on Interannual Variability of Yellow Sea Cold Water Mass. Institute of Oceanology, University of Chinese Academy of Sciences, Beijing, China. p. 23-26. (in Chinese with English abstract)
Li M, Zhong L J, Boicourt W C. 2005. Simulations of Chesapeake Bay estuary:sensitivity to turbulence mixing parameterizations and comparison with observations. Journal of Geophysical Research, 110(C12): C12004. DOI:10.1029/2004JC002585
Liu J, Ning P. 2011. Species composition and faunal characteristics of fishes in the Yellow Sea. Biodiversity Science, 19(6): 764-769. (in Chinese with English abstract)
Lü X G, Qiao F L, Xia C S, Wang G S, Yuan Y L. 2010. Upwelling and surface cold patches in the Yellow Sea in summer:effects of tidal mixing on the vertical circulation. Continental Shelf Research, 30(6): 620-632. DOI:10.1016/j.csr.2009.09.002
Mellor G L, Yamada T. 1982. Development of a turbulence closure model for geophysical fluid problems. Reviews of Geophysics, 20(4): 851-875. DOI:10.1029/RG020i004p00851
Mellor G. 2003. The three-dimensional current and surface wave equations. Journal of Physical Oceanography, 33(9): 1 978-1 989. DOI:10.1175/1520-0485(2003)033<1978:TTCASW>2.0.CO;2
Miao J B, Liu X Q, Xue Y. 1991. Study on the formational mechanism of the Northern Yellow (Huanghai) Sea cold water mass (I)-solution of the model. Science in China Series B, 34(8): 963-976.
Moon J H, Hirose N, Yoon J H. 2009. Comparison of wind and tidal contributions to seasonal circulation of the Yellow Sea. Journal of Geophysical Research, 114(C8): C08016.
Murphy A H. 1988. Skill scores based on the mean square error and their relationships to the correlation coefficient. Monthly Weather Review, 116(12): 2 417-2 424. DOI:10.1175/1520-0493(1988)116<2417:SSBOTM>2.0.CO;2
Noh Y, Ok H, Lee E, Toyoda T, Hirose N. 2016. Parameterization of Langmuir circulation in the ocean mixed layer model using LES and its application to the OGCM. Journal of Physical Oceanography, 46(1): 57-78. DOI:10.1175/JPO-D-14-0137.1
Qiao F L, Yuan Y L, Ezer T, Xia C S, Yang Y Z, Lü X G, Song Z Y. 2010. A three-dimensional surface wave-ocean circulation coupled model and its initial testing. Ocean Dynamics, 60(5): 1 339-1 355. DOI:10.1007/s10236-010-0326-y
Qiao F L, Yuan Y L, Yang Y Z, Zheng Q A, Xia C S, Ma J. 2004. Wave-induced mixing in the upper ocean:distribution and application to a global ocean circulation model. Geophysical Research Letters, 31(11): L11303.
Ren H J, Zhan J M. 2005. A numerical study on the seasonal variability of the Yellow Sea cold water mass and the related dynamics. Journal of Hydrodynamics, 20(S1): 887-896. (in Chinese with English abstract)
Shchepetkin A F, McWilliams J C. 2003. A method for computing horizontal pressure-gradient force in an oceanic model with a nonaligned vertical coordinate. Journal of Geophysical Research, 108(C3): 3 090. DOI:10.1029/2001JC001047
Shchepetkin A F, McWilliams J C. 2005. The regional oceanic modeling system (ROMS):a split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean Modelling, 9(4): 347-404. DOI:10.1016/j.ocemod.2004.08.002
Smyth W D, Skyllingstad E D, Crawford G B, Wijesekera H. 2002. Nonlocal fluxes and Stokes drift effects in the K-profile parameterization. Ocean Dynamics, 52(3): 104-115. DOI:10.1007/s10236-002-0012-9
Song X, Lin X P, Wang Y. 2009. The preliminary study of long-term variability of the Yellow Sea cold water in summer and its possible reasons. Journal of Guangdong Ocean University, 29(3): 59-63. (in Chinese with English abstract)
Su Y S, Su J. 1996. A preliminary study on the surface cold water in the Bohai and Yellow Sea during summer and its formation mechanisms. Acta Oceanologica Sinica, 18(1): 13-20. (in Chinese)
Umlauf L, Burchard H. 2003. A generic length-scale equation for geophysical turbulence models. Journal of Marine Research, 61(2): 235-265. DOI:10.1357/002224003322005087
Wang B D. 2000. Characteristics of variations and interrelations of biogenic elements in the Huanghai Sea Cold Water Mass. Acta Oceanologica Sinica, 22(6): 47-54. (in Chinese with English abstract)
Warner J C, Sherwood C R, Arango H G, Signell R P. 2005. Performance of four turbulence closure models implemented using a generic length scale method. Ocean Modelling, 8(1-2): 81-113. DOI:10.1016/j.ocemod.2003.12.003
Wijesekera H W, Allen J S, Newberger P A. 2003. Modeling study of turbulent mixing over the continental shelf:comparison of turbulent closure schemes. Journal of Geophysical Research, 108(C3): 3 103. DOI:10.1029/2001JC001234
Xin M, Ma D Y, Wang B D. 2015. Chemicohydrographic characteristics of the Yellow Sea Cold Water Mass. Acta Oceanologica Sinica, 34(6): 5-11. DOI:10.1007/s13131-015-0681-0
Yang H W, Cho Y K, Seo G H, You S H, Seo J W. 2014. Interannual variation of the southern limit in the Yellow Sea Bottom Cold Water and its causes. Journal of Marine Systems, 139: 119-127. DOI:10.1016/j.jmarsys.2014.05.007
Yao Z G, Bao X W, Li N, Li X B, Wan K, Song J. 2012. Seasonal evolution of the northern Yellow Sea Cold Water Mass. Periodical of Ocean University of China, 42(6): 9-15. (in Chinese with English abstract)
Yu F, Zhang Z X, Diao X Y, Guo J S, Tang Y X. 2006. Analysis of evolution of the Huanghai Sea Cold Water Mass and its relationship with adjacent water masses. Acta Oceanologica Sinica, 28(5): 26-34. (in Chinese with English abstract)
Yu X J, Guo X Y. 2018. Intensification of water temperature increase inside the bottom cold water by horizontal heat transport. Continental Shelf Research, 165: 26-36. DOI:10.1016/j.csr.2018.06.006
Zhang Y K, He X M, Gao Y F. 1983. Preliminary analysis on the modified water masses in the north Yellow Sea and the Bohai Sea. Transactions of Oceanology and Limnology, 2: 19-26. (in Chinese with English abstract)
Zhang Y K, Yang Y L. 1996. Analyses of the variational characteristics of the north Huanghai Sea Cold Water Mass. Marine Forecasts, 13(4): 15-21. (in Chinese with English abstract)
Zhao B R. 1985. The fronts of the Huanghai Sea Cold Water Mass induced by tidal mixing. Oceanologia et Limnologia Sinica, 16(6): 451-460. (in Chinese with English abstract)
Zhao B R. 1987. The continental shelf fronts induced by tidal mixing in the Huanghai Sea. Journal of Oceanography of Huanghai & Bohai Seas, 5(2): 16-23. (in Chinese with English abstract)
Zhao B R. 1989. The study on the basic characteristics and formation mechanisms of the intense thermocline in the Bohai Sea, Yellow Sea and the northern East Sea. Acta Oceanologica Sinica, 11(4): 401-410. (in Chinese)
Zhu J Y, Shi J, Guo X Y, Gao H W, Yao X H. 2018. Air-sea heat flux control on the Yellow Sea Cold Water Mass intensity and implications for its prediction. Continental Shelf Research, 152: 14-26. DOI:10.1016/j.csr.2017.10.006
Zieger S, Babanin A V, Erick Rogers W, Young I R. 2015. Observation-based source terms in the third-generation wave model WAVEWATCH. Ocean Modelling, 96: 2-25. DOI:10.1016/j.ocemod.2015.07.014