The Hamburg Shelf Ocean Model (HAMSOM) used in this work is a full-fledged three-dimensional primitive equation model with a free surface. The Cartesian coordinates are used with an Arakawa C-grid and a semi-implicit numerical scheme. The numerical framework of the ocean model was described by Backhaus (1985). The sea-ice model presented here is a two-class (ice and open water) model that takes into account both the ice dynamics and ice thermodynamics properties. It is a modification of the viscous-plastic model suggested by Hibler (1979) with three-layer thermodynamics according to Semtner (1976). The fundamental sea ice equations, thermodynamic properties, and rheology were introduced in detail by Ólason and Harms (2010) and Olason (2016).

The simulated area in this study is the Bohai Sea (see Figure 1 ), covering the latitude and longitude of 37°N~41.25°N and 117.25°E~122.5°E, respectively. The simulation period for an ice season covers from Nov. 1 to Mar. 31 on the next year, with a time step of 60 seconds and a horizontal grid resolution of 1’×1’ (1/60 arc degree, approximately 1.8 km in the meridional direction and 1.5 km in the zonal direction). There are 12 z-level grids used in vertical plane.

Bathymetry for the model domain is interpolated from the high-precision data provided by the Beihai Branch of the Ministry of Natural Resources of China. The model is initialized with monthly climatological salinity data from the OGCM for Earth Simulator (OFES) datasets (Sasaki et al., 2008) and a fused temperature data constructed by the Simple Ocean Data Assimilation (SODA) datasets (Carton et al., 2018) and a high-resolution satellite remote sensing assimilation SST (Group for High Resolution Sea Surface Temperature, GHRSST) (Donlon et al., 2012) using the “vertical projecting” assimilation method (Jia et al., 2022a). Due to the shallow depth of the Bohai Sea and the strong vertical mixing of seawater under the influence of energetic winds in winter, the dynamical steady state can be reached just a few tens of hours after the simulation starts (Liu, 2013), which means that after the simulation starts, the current velocity field in the Bohai Sea quickly adapts to the forcing of surface wind and open-boundary tides, and there is very little dependence on the initial current velocity and the initial elevation fields. Therefore, the initial current velocity and the initial elevation fields are set to zero. The harmonic constants of eight main tidal constituents at the open boundary are provided by The Oregon State University (OSU) Tidal Databases (Erofeeva et al., 2020).

The hourly meteorological data obtained from the global Climate Forecast System Reanalysis (CFSR) produced by the US Center for Environmental Prediction (NCEP) above the sea surface, including the 2 m air temperature, wind speed, cloud cover, precipitation rate, and specific humidity (Saha et al., 2011), are interpolated for the atmospheric forcing. Figure 3 displays the 2 m air temperature at the observation stations Tanggu and Platform-A (the positions of the two stations are labeled in Figure 1 ) during the winter of 2011/2012 and the corresponding values extracted from the NCEP dataset. The mean absolute errors at the two stations are 0.26°C and 0.25°C, respectively. This indicates that NCEP data can realistically reflect the weather system variations in the Bohai Sea during winter when large-scale cold air fronts occur.

Compared with coastal observation, sea survey, radar, and other observation methods, satellite remote sensing has the ability of large-scale, rapid, synchronous, and long-term continuous observation, and has unparalleled advantages in economy. In this study, the ice edge distance (i.e., the distance between the intersection of the central axis of Liaodong Bay and the ice edge line and the axis apex, see Figure 4A ) and the ice area observations from the satellite inversion, as well as the observed information for ice thickness and ice concentration in the three local waters (Erjiegou water (L1), Dabijiashan water (L2) and Xipaotai water (L3), the locations are shown in Figure 1 ) shown in Table 1 are obtained from the North China Sea Marine Forecasting Center of the Ministry of Natural Resources and its annual winter Sea Ice Land and Coastal Monitoring Report in the Bohai Sea. The satellite images that appeared in this paper are extracted from the NASA Worldview application (https://worldview.earthdata.nasa.gov), part of the NASA Earth Observing System Data and Information System (EOSDIS). The inversion data for ice thickness from GOCI ocean color remote sensing satellite shown in Figure 5 is obtained from Yan et al. (2024).

Mixed layer thickness determines the amount of heat content of the mixed layer, which in turn has an important influence on the hydrographic properties of the mixed layer (Wang et al., 2000). In the ice model, mixed layer thickness can influence the amount of sea ice freezing/melting at the ice-ocean interface by participating in the calculation of the variation rates of mixed layer temperature and salinity and the freezing rate of new-formed ice. In existing sea-ice simulations, the mixed layer thickness in the ice thermodynamics module is set to the thickness of the first water layer (Arthun and Schrum, 2010; Ólason and Harms, 2010; Olason, 2016; Lin and Zhao, 2019). This should be based on the premise that there is little difference between the first water layer thickness and the actual mixed layer thickness. Whereas the Bohai Sea is shallow, the vertical mixing of the water column is outright under the influence of high winds in winter, so the hydrographic properties of the water column are in vertical homogeneity from the sea surface to the seabed in the Bohai Sea (Liu, 2013; Yan et al., 2022). In other words, setting the mixed layer thickness to the first water layer thickness would cause a serious underestimation for the mixed layer thickness in the Bohai Sea during the ice season. In order to investigate the sensitivity of the ice simulations in the Bohai Sea to the mixed layer thickness, ice simulations in 2009/2010 ice season and 2011/2012 ice season are modeled using three different mixed layer thickness settings, including the first water layer thickness (default setting, 6 m), half of the water depth and the entire water depth.

The ice area results, defined here as the area of grid cells with at least a 15% sea ice concentration, of the above simulations are shown in Figure 6 , the time-averaged absolute errors of different cases (shown in the legend) are also calculated. It can be found that for both ice seasons, the simulated ice area obtained by setting mixed layer thickness to the entire water depth are overall increased and closer to the satellite observations than the other two cases, and the calculated time-averaged absolute errors are reduced in both ice seasons. It is calculated that for the 2009/2010 ice season, setting mixed layer thickness to entire water depth leads to a 10% improvement in the ice area simulations, compared to the default setting. While for the 2011/2012 ice season, 17% improvement is made.

To further clarify the effect mechanism of the mixed layer thickness to the sea ice variables, we present the time series of the ice coverage-averaged interfacial thermal parameters (including oceanic heat flux, conductive heat flux through ice, freezing/melting rate of the ice base (positive values represent ice freezing rates and negative values represent ice melting rates) and ice concentration) obtained by setting mixed layer thickness to the first water layer thickness and the entire water depth in the 2009/2010 ice season as shown in Figure 7 . It can be noticed that compared to setting mixed layer thickness as the first water layer thickness (blue lines), there is a significant decrease in the oceanic heat flux especially in the ice melting phase when setting it to the entire water depth (red lines), while the conductive heat flux through ice remains essentially unchanged, and the freezing/melting rate of the ice base grows slightly, as well as the ice concentration.

The effect mechanism of the mixed layer thickness on sea ice variables can be preliminary sorted out, based on the variations of the above parameters. For all grids except for those where the water depth equals the first water layer thickness, the mixed layer thickness setting to the entire water depth is deepened compared to the first water layer thickness. The deepening of the mixed layer thickness implies an increase in the total heat content of the mixed layer. This also means that when the mixed layer is changed by the same amount of heat, the magnitude of the change in the mixed layer temperature is reduced, i.e., the variation rate of the mixed layer temperature is reduced. As a result, the warming of the mixed layer due to the heating effect of solar radiation decreases during the ice melting phase, leading to a decrease in the positive deviation of the mixed layer temperature from the ice base temperature ( T m − T f in Equation 1), and a consequent decrease in the oceanic heat flux ( Figure 7A ) as measured by this temperature deviation. The oceanic heat flux and the conductive heat flux through ice jointly determine the heat budget at the ice-ocean interface. Since the conductive heat flux through ice is calculated from the vertical temperature gradient within the ice layer and is not directly constrained by the mixing layer temperature, the effect of the mixing layer thickness on the conductive heat flux through ice is negligible ( Figure 7B ). Therefore, the decrease in oceanic heat flux to the ice base implies a decrease in the amount of heat gained, manifesting net heat loss, resulting in an increase in the freezing rate of the ice base ( Figure 7C ), which ultimately leads to an increase in sea ice variables such as ice concentration ( Figure 7D ), ice thickness (same trend with the ice concentration, picture omitted) and ice area ( Figure 6A ). The corresponding mechanism schematic is shown in Figure 8 .

In summary, the deepening of the mixed layer somewhat corrects the heat content of the mixed layer, and hence the variation rates of the hydrographic properties within the mixed layer, and ultimately brings corrections to the sea ice variables in the Bohai Sea.

The oceanic heat flux is the key factor determining the heat balance at the ice-ocean interface (Wettlaufer, 1991). The oceanic heat flux Q w i can be estimated by the following bulk formula (McPhee, 1992):

where the ρ w is sea water density, c p is the specific heat capacity of sea water, u * is the ice-ocean friction velocity. T m is the mixed layer temperature and T f is the local freezing point (i.e., the temperature of ice base). C h is the HEC at the ice-ocean interface.

Although the well-established formulation scheme (Equation 1) of oceanic heat flux has been introduced in many ice models applied to the Bohai Sea, the HEC at the ice-ocean interface is still set as an empirical constant. This parameter characterizes the efficiency of heat transfer between the ice and the ocean and is related to turbulent and molecular motions within the ice-ocean boundary layer (McPhee, 1992; Sirevaag, 2009; Dansereau et al., 2014; Ha et al., 2016). Ji et al. (2002) surmised that the HEC varies with ice thickness and the hydraulic roughness of the ice under-surface based on the atmospheric, hydrographic and sea-ice observations from point stations in the JZ20-2 Sea Area of Liaodong Bay, 1997-1998. More specific physical processes involved in ice-ocean heat transmission could be ignored if setting the HEC to a constant (McPhee, 1992). Therefore, the relationship between the HEC and ice thickness and the ice-ocean friction velocity based on the parameterization proposed by McPhee (1992) was introduced into our model to obtain the time and space-varying HEC in our previous study (Jia et al., 2022b). The HEC parameterization formulas are as follows:

Here, the ϕ t u r b and ϕ m o l are the non-dimensional temperature changes across the turbulent boundary layer and the molecular transition sublayer embedded within the turbulent boundary layer (see Figure 2 ), respectively. Turbulent motions caused by shear forces of ocean currents under ice dominate in the turbulent boundary layer, while molecular viscosity and diffusivity motions dominate in the molecular transition sublayer. The Von-Karman’s constant k is 0.4. ξ N is a constant, was approximated to 0.05 by McPhee (1992) at the time the parameterization was proposed which applied in the Arctic Ocean. Dansereau et al. (2014) also set it to 0.052 in the ice simulations of the Pine Island ice shelf, West Antarctica. In this paper, it is taken to be 0.05. η * is a buoyancy factor, equal to 1 when the stratification is neutrally stable (McPhee, 1992). f is the Coriolis parameter. z 0 is the hydraulic roughness of the ice under-surface, u * is the ice-ocean friction velocity. The used calculation schemes of both z 0 and u * were proposed by Mellor and Kantha (1989) (see Equations 5, 8 in the following content, respectively). b is an empirical constant, 0.6, taken as the approximation to the laboratory result of Yaglom and Kader (1974), following McPhee (1992). R e is the Reynolds number, R e = u * h m o l ν , where h m o l is the thickness of molecular transition sublayer, h m o l = 30 z 0 (Yaglom and Kader, 1974), v is the kinematic molecular viscosity (v= 1.8 × 10 − 6 ). The molecular Prandtl number P r is 12.9.

Our previous findings indicated that the ice simulation accuracy was effectively improved when using the above parameterized HEC than using constant HEC, especially for the ice melting phase. After HEC parameterization, the time-averaged absolute error of the simulated ice area reduced from 1855 km² (HEC=0.001, default setting) to 1154 km². More detailed descriptions of the difference in the results between using the HEC parameterization and constant values were provided in Jia et al. (2022b).

To further improve the applicability of the above HEC parameterization to the Bohai Sea ice simulations, the new work in this study over our previous work (Jia et al., 2022b) lies in evaluating the calculation schemes for the two core parameters, the hydraulic roughness of the ice under-surface z 0 and the ice-ocean friction velocity u * . For the hydraulic roughness of the ice under-surface, in addition to the previously used Mellor and Kantha (1989) scheme, another two schemes proposed by McPhee (1992) and Dansereau et al. (2014), respectively, are also induced here to detect the model performance of the Bohai Sea ice. The formulas of the three schemes assessed in this study are:

(Mellor and Kantha, 1989),

(McPhee, 1992),

(Dansereau et al., 2014),

where h i is the ice thickness, h 1 is the depth at which the first water layer located. v i → is the ice velocity, and v w → is the mixed layer velocity under the ice.

For the ice-ocean friction velocity, two calculation schemes proposed by Mellor and Kantha (1989) (used previously) and McPhee et al. (2016) are assessed here. The two schemes formulas are:

(Mellor and Kantha, 1989)

(McPhee et al., 2016),

where C d w is the ice-ocean drag coefficient, which is set as 0.0055 following Fujisaki et al. (2011).

The simulated ice area using different hydraulic roughness calculation schemes and different friction velocity calculation schemes are compared respectively in Figure 9 .

As we can see from Figure 9A , the hydraulic roughness calculation scheme proposed by Mellor and Kantha (1989) give the best model performance, the two fluctuations occurred from Feb. 16 to 20 and Feb. 22 to Mar. 1, 2012 are well reflected, however, the other two schemes not captured. And Figure 9B shows that the friction velocity calculation scheme proposed by McPhee et al. (2016) has slightly better model performance than the other scheme. Thus, the hydraulic roughness of the ice under-surface z 0 calculation scheme still use the scheme proposed by Mellor and Kantha (1989), and the ice-ocean friction velocity u * calculation scheme proposed by McPhee et al. (2016) is selected to apply in the HEC parameterization.

In order to check the effectiveness of the updated HEC parameterization with new ice-ocean friction velocity calculation scheme in improving the model performance, the simulated ice area obtained by a constant HEC (0.001, model default setting), using the previously HEC parameterization in Jia et al. (2022b), and using the updated HEC parameterization, respectively, are compared in Figure 9C . Our simulations suggest that the HEC parameterizations can improve the accuracy of ice simulations to a large extent compared with the simulation using constant HEC value (0.001). Significant improvement is seen during the ice melting phase, with two notable fluctuations satisfactorily captured using the parameterized HEC. And the updated HEC parameterization with new ice-ocean friction velocity calculation scheme leads to the best performance, with its time-averaged absolute error decreasing further from the previously HEC parameterization, from 1154 km² to 1132 km². Therefore, the updated HEC parameterization is applied in the subsequent numerical simulations of the Bohai Sea ice in this study.

To comprehensive verify the ice simulations in the 2011/2012 ice season obtained after the adjustments of the two key parameters in the ice thermodynamics module, firstly, the simulated ice thickness and ice concentration in three local waters [Erjiegou water (L1), Dabijiashan water (L2) and Xipaotai water (L3)] on the top of Liaodong bay ( Figure 10 ) are validated by the corresponding observed information during the ice growth, severe, and melting phases ( Table 1 ) collected from the Bohai Sea Ice Land and Coastal Monitoring Report in the 2011/2012 ice season. It can be seen that the ice thickness time series during the whole ice season at the three waters basically experienced multiple fluctuations with varying magnitude. Compared with the observed information, the simulated ice thickness and ice concentration in the three local waters are generally within reasonable limits, however, the ice concentration in Xipaotai water during the ice severe phase was underestimated by about 30%.

Then, three model days (Jan. 05, Jan. 23, and Feb. 26 in 2012) are selected to compare the simulated ice extent and the corresponding satellite images ( Figure 4A ), the calculated ice edge distances are also labeled in the lower left corner of each graph. Since sea ice mainly exists in Liaodong Bay in the 2011/2012 ice season, only the Liaodong Bay area is shown in the figure. Our results indicate that both the simulated ice extent and ice edge distance are basically consistent with the satellite data, the errors of the simulated ice edge distances on all three days did not exceed 2 nautical miles. A comparison of the time series of the simulated ice edge distance ( Figure 4B ) and ice area ( Figure 4C ) with the observations from the satellite inversion, as well as their error analysis ( Table 2 ) also suggest that the temporal trend of ice variables is in good agreement with the satellite inversion observations. Quite significant correlations are found between the simulated ice edge distance (Pearson’s correlation coefficient r=0.92) and ice area (Pearson’s correlation coefficient r=0.98) with the corresponding satellite observations. In summary, the satisfactory model performance in 2011/2012 ice season allows to perform in-depth thermodynamic analysis of the ice-ocean interface.

To interpret the temporal variation characteristics of the heat fluxes at the ice-ocean interface in the Bohai Sea, the time series of the atmospheric conditioning factors (including sea surface air temperature and wind speed), the oceanic conditioning factor (i.e., mixed layer temperature), the ice variables (including ice drift velocity, ice thickness, ice concentration, and the freezing/melting rate of ice base), and the interfacial heat fluxes (including the conductive heat flux through ice, oceanic heat flux, and the latent heat flux) during the 2011/2012 ice season are shown in Figure 5 . The correlation coefficients between the heat fluxes with the atmospheric, oceanic and ice conditioning factors are also calculated and listed in Table 3 .

To investigate the spatial variation characteristics of the heat fluxes at the ice-ocean interface and compare their differences between ice growth state and ice melting state, the distributions of the key interfacial thermal parameters on two representative days, with one day in the ice growth phase and the other in the ice melting phase, are illustrated in Figures 11 , 12 , respectively.

Firstly, the result shows that the oceanic heat flux is always positive both in ice growth state and ice melting state. This is because the oceanic heat flux represents the unidirectional heat transfer from the mixed layer to the ice base and is measured by the mixed layer temperature elevation above freezing (calculated only when the mixing layer temperature is greater than the freezing point). The value of the oceanic heat flux is larger in the marginal ice zone (MIZ) (about 20~80 W/m²) with relatively thin ice than in the inner ice zone (IIZ) (< 20 W/m²) with relatively thick ice, and the value in the nearshore area in the northern especially northeastern part of Liaodong Bay is also relatively large when in the ice-melt phase (shown in Figures 11A , 12A ). This phenomenon can be explained from both dynamic and thermodynamic aspects. From the aspect of dynamics, the mixed layer under the MIZ with lower ice concentration (i.e., more polynyas) is more susceptible to the wind energy input, which strengthens the ice drift velocity (shown in Figures 11F , 12F ) and turbulent mixing, further promoting heat exchange within the ice-ocean boundary layer. From the aspect of thermodynamics, the polynyas in the MIZ make the underlying mixed layer more susceptible to atmospheric thermal driven, and the mixed layer temperature can be heated above freezing by more solar radiation (shown in Figures 11E , 12E ).

We can also see that in the ice growth state, the value of the conductive heat flux through ice in the whole field is mostly positive (shown in Figure 11B ), which means that the heat is transferred upward from the ice base. The magnitude in MIZ (about 10~18 W/m²) is smaller than in the IIZ (about 20~40 W/m²). While in the ice melting state, the value of the conductive heat flux through ice in the whole field is mostly negative (shown in Figure 12B ), which means that heat is being transferred downward. The magnitude in MIZ (about 22~28 W/m²) is larger than in the IIZ (<20 W/m²). The simulated value range here is relatively consistent with that in the Arctic ocean. Fan et al. (2017) computed the conductive heat flux through ice in the Arctic Ocean in the four winter months (November–February) for a long period of 36 years (1979–2014) based on monthly mean Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) sea ice thickness fields, and found that the value range is about 10~30 W/m², most values fall in the 15~25 W/m².

The change of latent heat flux is the result of competition between the conductive heat flux through ice and oceanic heat flux, characterizing the heat gain or loss at the ice-ocean interface. During the growth of the ice, the latent heat flux in the whole field is almost positive (shown in Figure 11C ), indicating that the heat is transmitting upward from the ice base, and the ice-ocean interface is losing heat. The degree of this heat loss in the IIZ is obviously greater than that in the MIZ. During the ice melting phase, the latent heat flux in the whole field is negative (shown in Figure 12C ), indicating that the heat is transmitting downward to the ice base, and the ice-ocean interface is gaining heat. The degree of this heat gain in the MIZ and nearshore area is greater than that in the IIZ.

As for the freezing/melting rate of ice base, its spatial distribution closely resembles the latent heat flux (shown in Figures 11D , 12D ). During ice growing and melting, the whole field is positive (freezing) and negative (melting), respectively. As shown in Figures 11G , 12G , the ice thickness in the MIZ is thinner than the IIZ. The simulated ice coverage and ice edge line are generally confirmed by satellite images. The simulated spatial distribution characteristics of ice thickness in the ice growth state are also consistent with the satellite image, while in the ice melting state, the simulated ice thickness in Liaodong Bay is not well characterized by thick in the east side and thin in the west side. Since the spatial distribution of ice thickness in the Liaodong Bay is highly influenced by its clockwise circulation field (Wu et al., 2005), the ice dynamic process including mixed layer currents, the ice velocity needs to be further studied and improved in the future.

As stated the influence mechanism of mixed layer thickness on ice variables described in section 3.1, the mixed layer thickness can effectively affect the heat budget at the ice-ocean interface by determining the variation rates of hydrographic properties within the mixed layer (Wang et al., 2000). Therefore, it is important to define reasonable values for the mixed layer thickness in the sea ice module that meets the environmental conditions of the local area. For sea ice simulations in the areas with significant mixed layer seasonal or spatial variation, setting mixed layer thickness as the thickness of the first water layer (the default setting) in the ice model would bias the calculation of the variation rates of the hydrographic properties of the mixed layer, especially in z-level height coordinate models. Due to the vertical mixing homogenization of the water column during winter, the mixed layer thickness at each grid determined using the temperature difference threshold method (the difference between the temperature at the depth and the surface temperature is less than 0.2°C) in the Bohai Sea is almost equal to its water depth (figure omitted). So we directly make the mixed layer thickness the entire water depth in the Bohai Sea ice simulations, which is a reflection of the special atmospheric and oceanic environmental conditions of the Bohai Sea. Wang et al. (2000) indicated that the hydrographic properties under ice as well as the oceanic heat flux at the ice-ocean interface are modulated by bathymetry in the Bohai Sea. A consistency of the sea ice edge line with the isobath contour line in the Bohai Sea provides support for this finding (Li et al., 2020; Yan et al., 2022). For other marginal seas such as the Baltic Sea and Sea of Okhotsk, they also experience strong seasonal sea ice with different bathymetry, atmospheric and oceanic environmental conditions. The mixed layer thickness in the Sea of Okhotsk is characterized by significant seasonal and spatial variations, a mixed layer of uniform temperature nearly at the freezing point extending down to a depth of about 300 m was observed in winter of the southwestern of the Okhotsk Sea (Ohshima et al., 2001), which is far deeper than the first water layer thickness in general models. Seasonal variations are also a major feature of mixed layer thickness in the Baltic Sea, about half of the water column in the Baltic Sea is dominated by the seasonal signal directly because mixed layer thickness varies in the range of 20-30 m at a mean water depth of only 52 m (Meier et al., 2003). Considering the unique oceanographic characteristics of these seas, we suggest to use a temperature difference threshold method to determine the mixed layer thickness at each grid in order to more accurately determine the extent of the hydrodynamic mixing of the water column under the ice.

An updated parameterization of the heat exchange efficiency (i.e., HEC) between the ice and the ocean is applied in the Bohai Sea ice simulations in this study. Here, the HEC at the ice-ocean interface is defined as the non-dimensional temperature changes across the turbulent boundary layer and the molecular transition sublayer, considering the impact of the ice-ocean friction velocity and the hydraulic roughness of the ice under-surface. The values of oceanic heat flux modeled by Su et al. (2005) using the constant HEC (2.2×10^-5) in the MIZ of the Liaodong Bay range from 30 W/m² to 100 W/m², which are relatively larger than our values (about 20~80 W/m²) using the parameterized HEC. Using the oceanic heat (about 20 W/m² to 50 W/m²) provided by the Yellow Sea warm current in winter (Wang et al., 1984) as a reference (Wang et al., 2000), the oceanic heat flux modeled in this study has a more reasonable value. Many ice models for simulating ice shelf in polar regions have adopted this HEC parameterization (e.g., Little et al., 2009 [HIM]; Holland et al., 2010 [MICOM/POLAIR]; Makinson et al., 2011 [MICOM]; Timmermann et al., 2012 [FESOM] and Galton-Fenzi et al., 2012 [ROMS]). Nevertheless, constant HEC is still adopted in some ice models. From the simulations shown in Figure 9C , the HEC parameterization scheme has a better performance and is able to capture more detailed ice fluctuations during the ice melting phase. Dansereau et al. (2014) also found that important differences appeared in the ice melting rate distribution under the Pine Island ice shelf, West Antarctica between simulations using constant HEC and parameterized HEC. They indicated that the regions of largest melting coincide with strong outflow plumes and fast mixed layer currents when using the parameterized HEC, which cannot be simulated using constant HEC. The role of the ice-ocean friction velocity is highlighted, they also pointed that the turbulence measured by friction velocity dominate over temperature gradient in setting oceanic heat flux through the ice-ocean boundary layer. Therefore, a careful evaluation of the calculation scheme for ice-ocean friction velocity is necessary. We propose that a dynamic HEC should be used in modeling marginal ice zones where melting occurs frequently.

To our knowledge, the existing studies of ice-ocean interfacial heat fluxes in the Bohai Sea mainly focus on the oceanic heat flux. The main conclusions from existing studies include: 1) the evolution of the Bohai Sea ice is more dominated by ocean heat content (i.e., oceanic heat flux) of the mixed layer than salt (Yan et al., 2022); 2) the oceanic heat flux in the Bohai Sea has large value at the beginning of the ice season and then decreases as the ice season progresses (Ji et al., 2002; Su et al., 2005), and 3) the oceanic heat flux is controlled by the mixed layer thickness and mixed layer temperature (Wang et al., 2000).

These findings are also reflected in our results. Further, we found that the oceanic heat flux in the Bohai Sea rises again at the end of the ice season after reaching a minimum (close to zero) during the middle ice season. Su et al. (2005) obtained the similar result based on their numerical results of the Bohai Sea ice in the 2000/2001 ice season using the POM ice-ocean coupled model. The main difference between our results and theirs, however, is that values of the oceanic heat flux we modeled are larger at the end of the ice season than at the beginning. The stronger rebound of the oceanic heat flux at the end of the ice season is due to the stronger rewarming of air temperatures at this period ( Figure 5A ), causing a greater increase in the mixed layer temperature elevation above freezing ( Figure 5C ) and thus the heat transferred from the ocean to the ice ( Figure 5D ). Yang (2015) obtained an optimal solution for the oceanic heat flux in the Arctic Ocean based on an optimal identification model for the ice base sublayer energy balance system, using an improved genetic mathematical algorithm. Their results show that the oceanic heat flux generally declines from the early stage of ice growth, and increases when entering the ice melt period, until reaching the maximum value when the ice melts completely. This finding is consistent with our simulation results.

In addition, the new outcomes in this study which haven’t been reported in existing studies for the Bohai Sea include:

The model parameter adjustments and simulations analysis in this study only focus on the ice-ocean thermal process. Ice dynamic processes need to be evaluated or improved in the future. In particular, the value/calculation of the ice-ocean drag coefficient (C_dw ) needs to be further evaluated as it affects both ice thermodynamic and dynamic processes. Constant ice-ocean drag coefficient is usually employed in sea ice simulations (Dansereau et al., 2014). Holland and Feltham (2006) set C_dw to 0.0015 at the ice shelf base to account for smoothing effects by melting and ice pumping. Dansereau et al. (2014) found that reasonable cavity-averaged melt rate under PIIS can be matched when using C_dw in the range of 0.005-0.01. In this study, a C_dw value of 0.0055 following Fujisaki et al. (2011) is adopted in the HEC parameterization. The ice-ocean drag coefficient is a highly uncertain parameter, related sensitivity studies for this coefficient will be followed up in the future.