NEMO-Bohai 1.0: a high-resolution ocean and sea ice modelling system for the Bohai Sea, China

Severe ice conditions in the Bohai Sea could cause serious harm to maritime traffic, offshore oil exploitation, aquaculture, and other economic activities in the surrounding regions. In addition to providing sea ice forecasts for disaster prevention and risk mitigation, sea ice numerical models could help explain the sea ice variability within the context of climate change in marine ecosystems, such as spotted seals, which are the only ice-dependent animal that breeds in Chinese waters. Here, we developed NEMO-Bohai, an ocean–ice coupled model based on the Nucleus for European Modelling of the Ocean (NEMO) model version 4.0 and Sea Ice Modelling Integrated Initiative (SI³) (NEMO4.0-SI³) for the Bohai Sea. This study will present the scientific design and technical choices of the parameterizations for the NEMO-Bohai model. The model was calibrated and evaluated with in situ and satellite observations of the ocean and sea ice. The model simulations agree with the observations with respect to sea surface height (SSH), temperature (SST), salinity (SSS), currents, and temperature and salinity stratification. The seasonal variation of the sea ice area is well simulated by the model compared to the satellite remote sensing data for the period of 1996–2017. Overall agreement is found for the occurrence dates of the annual maximum sea ice area. The simulated sea ice thickness and volume are in general agreement with the observations with slight overestimations. NEMO-Bohai can simulate seasonal sea ice evolution and long-term interannual variations. Hence, NEMO-Bohai is a valuable tool for long-term ocean and ice simulations and climate change studies.

1 Introduction

The Bohai Sea is the southernmost seasonal frozen sea in the Northern Hemisphere (Yan et al., 2017). The formation of sea ice in the Bohai Sea mainly depends on the geographical environment and hydrometeorological characteristics (Ding, 1999). Specifically, the Bohai Sea is located in the continental shelf area, and the average water depth is only 18 m (Su and Wang, 2012), which indicates low oceanic heat content in winter. Following a northern continental climate, the Bohai Sea is affected by the cold Siberian air every winter, which causes the sea surface temperature of the Bohai Sea to be significantly lower than that at the same latitude (Zhang et al., 2016; Donlon et al., 2012). In addition, the sea surface salinity of the Bohai Sea is about 30 PSU, which is the lowest in the entire coastal waters in China (Yan et al., 2020). It means that the Bohai seawater freezes before reaching the maximum density. Therefore, it even more easily convects and loses heat just before freezing. The sea ice season is generally from December to March. The sea ice in the Bohai Sea also exhibits significant interannual variability, which could cover half of the sea area during severe winters (Su and Wang, 2012).

With 4 of China's top 10 ports (Tianjin, Tangshan, Dalian, Yingkou), the Bohai Sea is an important economic zone in China (Fu et al., 2017). Under severe sea ice conditions, such as during the 2009/2010 winter, 296 ports in this region were frozen. In particular, the 42 km artificial waterway of Huanghua Port was covered by sea ice; thus, the maritime traffic was severely restricted, and 7157 fishing vessels were damaged (Gu et al., 2013). About 2×10⁵ ha offshore aquatic farms were covered by sea ice, and the aquaculture industry losses reached ∼ CNY 6 billion (∼ USD 878 million) (Zhang et al., 2013). In addition, ice floes drift on the Bohai Sea surface under the drive of winds, tides, and ocean currents. When the ice floes meet offshore structures (oil and gas platforms, navigation lights, and lighthouses) or ships, convergent ice motion seriously threatens marine transportation and oil and gas exploration (Leppäranta, 2011; Li et al., 2011; Yan et al., 2019). Thus, sea ice monitoring and forecasting are essential for assessing sea ice hazard risk and preventing ice disasters.

In addition, spotted seals (Phoca largha) are listed as a least concern species by the International Union for the Conservation of Nature Red List of Threatened Species and as a Grade II state protection species by China's Wild Animal Protection Law (1988). They are distributed in the Northwest Pacific Ocean with Liaodong Bay in the Bohai Sea as their southernmost geographic breeding site. At the end of October each year, the seals travel from the Yellow Sea to Liaodong Bay. However, the population of the spotted seal has experienced several drastic declines, from ∼2300 to ∼1000 individuals, over the past 3 decades due to hunting and environmental pollution (Yan et al., 2018; Yang et al., 2017). They den in coastal areas and on offshore drifting ice. When the sea ice melts, they gradually leave Liaodong Bay and embark on a long journey back to the Pacific Ocean. Therefore, it is crucial to study long-term changes of Bohai Sea ice to study its impacts on spotted seals.

Most research has focused on the polar region, however, sea ice at middle to high latitudes is also sensitive to climate change (Bai et al., 2011; Gong et al., 2007; Yan et al., 2017). Knowledge of regional sea ice variability is vital to studying how climate change will affect regional basins. The traditional view on the Bohai Sea ice is concentrated on extracting sea ice data from satellite imagery (Karvonen et al., 2017; Shi and Wang, 2012; Su et al., 2019). As the use of remote sensing data to retrieve sea ice data has limitations in providing a continuous record for a long time series, numerical modelling is an effective tool to understand the sea ice development from freezing to thawing (Hordoir et al., 2019; Pemberton et al., 2017).

Sea ice models can reasonably simulate sea ice conditions, which can primarily supplement remote Earth observations. To our knowledge, there are only very few studies on developing regional coupled ocean–ice models for the Bohai Sea, and most of them have focused on short-duration simulations, such as 1-week or 1-year case studies. Wang et al. (1984) reported the first sea ice dynamic–thermodynamic model for the Bohai Sea, which simulated the sea ice freezing–melting cycle at a resolution of 20 km × 20 km driven by monthly climatic forcing. Wu (1991) and Wu et al. (1997) tried to simulate sea ice evolution based on a dynamic–thermodynamic model with consideration of an ice thickness distribution function consisting of three types (level ice, drift ice, and open water). Bai and Wu (1998) combined the sea ice dynamic–thermodynamic model, atmospheric boundary layer model, and tide model at a resolution of 8.64 km × 11.11 km. Liu et al. (2003) applied a high-resolution atmosphere model to reproduce sea ice conditions in 2002 using a sea ice model initialized by high-resolution satellite data from HY-1 and MODIS. Su et al. (2004) simulated the winters of 1998 and 2000 with the coupled Princeton Ocean Model (POM) and a viscous–plastic dynamic–thermodynamic sea ice model. Tang et al. (2010) simulate the sea ice freezing–melting process in 2003 with a coupled POM and a sea ice model considering the thickness distribution function. Zhang and Zhang (2013) used a finite-volume, primitive equation community ocean model to simulate sea ice in 2003. Li et al. (2021) simulated sea ice variations in the severe winter of 2009–2010 in the Bohai Sea based on version 3.6 of the Nucleus for European Modelling of the Ocean (NEMO). To the best of our knowledge, there has been no long-term simulation study on Bohai Sea ice so far that aims to understand the interannual or interdecadal variability of the area.

NEMO is a state-of-the-art numerical modelling framework designed for research activities and forecasting services in the ocean and climate sciences (Madec et al., 2016). The NEMO code, including its sea ice component, enables the investigation of ocean–ice dynamics and thermodynamics and their interactions with the atmosphere. NEMO offers a wide range of applications from short-term forecasts and climate projections (Drouard and Cassou, 2019; Obermann-Hellhund et al., 2018; Uotila et al., 2017; Voldoire et al., 2013) to process studies (Courtois et al., 2017; Declerck et al., 2016; Feucher et al., 2019). While global models offer a poor representation of regional ocean processes, regional models have been developed based on the framework of NEMO in recent years for various seas, e.g. the Atlantic marginal basins (Graham et al., 2018; O'Dea et al., 2017), the Baltic and North seas (Hordoir et al., 2019; Pemberton et al., 2017), the northwestern Mediterranean Sea (Declerck et al., 2016), the southern Indian Ocean (Schwarzkopf et al., 2019), the South China Sea (Thompson et al., 2019), and the Black Sea (Gunduz et al., 2020). However, a NEMO-based regional model for the Bohai Sea has not been attempted for long-term climate studies until now.

Long-term sea ice simulations could help improve our understanding of thermodynamic and dynamic sea ice processes in the Bohai Sea, crucial for sea ice disaster prevention, spotted seal habitat studies, and regional climate change studies. The paper aims to report the development of NEMO-Bohai and assess its performance. The paper is organized as follows: Sect. 2 introduces the model and observation dataset. Model comparison and validation are carried out in Sect. 3. The analysis of sea ice variation based on the 22-year hindcast simulations of NEMO-Bohai is presented in Sect. 4. A discussion and summary are provided in Sect. 5.

2 Model setup: NEMO-Bohai

NEMO-Bohai is localized to the Bohai Sea based on the NEMO ocean engine (Madec et al., 2016). We apply the NEMO 4.0 beta and Sea Ice Integrated Modelling Initiative (SI³) model in a regional configuration covering the Bohai Sea. The regional configuration, NEMO-Bohai, is nested into the global configuration using one-way lateral boundary conditions with a flow relaxation scheme algorithm (Engedahl, 1995). The global $\mathrm{1}/\mathrm{4}$^∘ ocean–ice model (ORCA025) maintained by the DRAKKAR group is performed. A detailed description of ORCA025 is provided by Barnier et al. (2006). In our simulation, forcing fields are provided from the Japanese 55-year Reanalysis (Harada et al., 2016; Kobayashi et al., 2015), covering the 55 years from 1958 with $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ spatial resolution. A monthly climatological runoff flux with Antarctic ice shelves and iceberg meltwater flux (Dai et al., 2009; Depoorter et al., 2013) is used. To limit the atmospheric forcing biases from propagating into the upper ocean and avoid salinity and temperature drift, which both impact the overturning circulation, the sea surface temperature (SST) and salinity (SSS) are restored monthly based on the World Ocean Atlas 2013 (Locarnini et al., 2013; Zweng et al., 2013).

The global simulation is spun up in January 1987 with the rest state and ended in 2017. The time step is set to 900 s. The ORCA025 configuration is performed using 240 CPU cores on a Cray XC40 system on the Sisu supercomputer and requires ∼6000 CPU hours per simulated year. The model output frequency is set to every 5 d from 1987 to 1994 for the test run and increases to every day for the period of 1995–2017. The size of the output files per simulated year is approximately 2 TB. The open boundary for NEMO-Bohai is forced by barotropic and baroclinic modes. The sea surface heights (SSHs), temperatures, salinity levels, and velocities are relaxed to outer-model values over a single-point zone at the unstructured open boundary of the model domain. In addition, the use of the XIOS server makes NEMO-Bohai more user-friendly.

2.1 NEMO-Bohai: The ocean model

The NEMO-Bohai domain consists of one central zone and three bays (see Fig. 1). The area covers 37–41^∘ N and 117–122^∘ E. The horizontal resolution is 0.025^∘ in orthogonal curvilinear coordinates, which is equivalent to a resolution of approximately 2 km ($\mathrm{204}\times \mathrm{244}=\mathrm{49}$ 776 grid points). The bathymetry interpolates into the target grids according to the ETOPO1 1 arcmin dataset. We further manually remove the isolated land cells with a number less than 4. The vertical discretization scheme adopts the z^∗ formulation. There are 10 vertical levels, and the top layer is set to 1 m. The partial step is applied to the bottom topography. The time step of the ocean model is set to 90 s, while it calls the sea ice model every three time steps. The split-explicit time-stepping method (ln_dynspg_ts = true.) is applied to compute barotropic and baroclinic modes and their interactions.

Figure 1 The location and bathymetry of the Bohai Sea. Black dots denote the coastal tide gauge stations. Red dots and black triangles indicate station locations with sea surface temperature and sea surface salinity observations, respectively. Black squares denote the oil platforms where the sea ice thickness observations were conducted. Line A–B represents the profile section along which temperature and salinity were presented (Figs. 6 and 7, respectively), while line C–D indicates the Bohai Strait. The straight line in the east denotes an open boundary for the regional ocean model.

The initial conditions of the numerical simulation in the Bohai Sea for 1 July 1995, including SST and SSS, are obtained by interpolating the temperature and salinity fields from the global ocean simulations with the ORCA025 configuration, as shown in Fig. 2. The regional model is forced with lateral ocean boundary conditions, tidal forcing, atmospheric forcing, and river runoff during the study period. Two kinds of boundary conditions are used for the setting of lateral open ocean boundary. A flow relaxation scheme is applied to baroclinic velocities and tracers (Engedahl, 1995), while a Flather boundary condition (Flather, 1976) is used for barotropic dynamics, such as SSH and barotropic velocities (u2d, v2d). NEMO-Bohai has one open boundary located 100 km away from the Bohai Strait, and we set the relaxation zone width to 1 (nn_rimwidth = 1). The lateral open-ocean boundary conditions for 1995 to 2017 are also taken from a series of global ocean simulations with the ORCA025 hindcast. The Oregon State University Tidal Inversion Model is used for the barotropic mode (Egbert and Erofeeva, 2002), including 11 tidal harmonics (M2, S2, N2, K2, K1, O1, P1, Q1, M4, MS4, MN4). We choose the TPXO9.1 global tide model, which is the latest version with a $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ resolution.

Figure 2 Examples of (a) an initial sea surface temperature (SST) field and (b) an initial sea surface salinity (SSS) field interpolated from the global ORCA025 model simulation on 1 July 1995.

The NEMO-Bohai atmospheric forcing is derived from the ERA5 dataset, the ECMWF's latest reanalysis product covering the period from 1979 to the present, which has replaced the widely used ERA-Interim dataset. The data cover the Earth with a horizontal spatial resolution of 30 km and represent the atmosphere using 137 levels from the surface up to a height equaling 0.01 hPa. The forcing files consist of 1-hourly instantaneous fields of the 2 m air temperature, 10 m wind speed, downward short-wave and long-wave radiation, sea level pressure, specific humidity, and precipitation. The equivalent inverse barometer SSH is calculated from the atmospheric pressure (ln_apr_dyn = true.) with the referenced pressure of 101 000 N m⁻². In addition, the river runoff, which provides a significant volume of freshwater to the Bohai Sea and influences stratification, circulation, and primary productivity, is considered in the model. The daily runoff volume data for the rivers flowing into the Bohai Sea, such as the Yellow River, Liao River, and Hai River, are acquired from JRA-55 from 1995 to 2018. The runoff forcing is considered freshwater with a salinity of 0 PSU and the same temperature as the ocean surface.

The bottom roughness impacts the dynamics of the tide, ocean circulation, and storm surges in the Bohai Sea. A constant bottom roughness (rn_z0) is used to calculate the drag coefficient for all bottom grids. To better reflect changes in the tide and circulation, we tune the bottom roughness accordingly. Regarding the configuration in the Gulf of Finland (rn_z0 = $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{3}}$ m) (Westerlund, 2019), we set rn_z0 to $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{1}}$, $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{2}}$, $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{3}}$, $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{4}}$, $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{5}}$, and $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{8}}$ while keeping all the other variables consistent. Results show that reducing the bottom friction can increase the phase difference and amplitude of the tide. In NEMO-Bohai, the smaller the roughness value, the more accurate the tide. However, it also needs to be combined with the actual situation at the bottom of the Bohai Sea. With reference to Zhang and Zhang (2013), the bottom roughness had a magnitude of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ m. Therefore, the bottom roughness in NEMO-Bohai is determined to be $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{3}}$ m. The bottom drag coefficient (rn_Cd0) is set to the default value after a set of experiments. The turbulent kinetic energy (TKE) closure scheme is used for vertical mixing (Blanke and Delecluse, 1993) with the default values of vertical eddy viscosity and diffusivity coefficients, and the flux-corrected transport scheme (Zalesak, 1979) is responsible for the tracer advection. Some critical parameters of the ocean part of NEMO-Bohai are listed in Table 1.

Figure 3 Comparison of sea surface height (SSH) between the observations and NEMO-Bohai simulations at (a) Yantai, (b) Dalian, (c) Jingtanggang, and (d) Bayuquan stations in the Bohai Sea.

Table 1 Key physical parameters in the ocean namelist (namelist_oce) in NEMO-Bohai.

Figure 4 The simulated sea surface temperature (SST) compared with in situ observations at Laohutan, Dalian (a–h), and Zhifudao, Yantai (i–p), for the periods of 2000–2001 and 2010–2015.

2.2 NEMO-Bohai: The sea ice model

SI³ is recommended by the Sea Ice Working Group to reduce duplication and make better use of development resources (Aksenov et al., 2019). SI³ merges the capabilities of three sea ice models previously used in NEMO (CICE, GELATO, and LIM). SI³ is a horizontal dynamic and vertical thermodynamic sea ice model. The thermodynamics and dynamics are separated and rely upon different frameworks and sets of hypotheses: thermodynamics uses the ice thickness distribution and the mushy-layer frameworks, whereas dynamics assumes continuum mechanics, e.g. Leppäranta (2011). Thermodynamics and dynamics interact in two ways: first, advection impacts the ice state variables; second, the horizontal momentum equation depends on, among other things, the ice state. The modified elastic–viscous–plastic rheology is used for sea ice dynamics (Bouillon et al., 2013; Kimmritz et al., 2017). Bohai Sea ice is relatively thin and has low salinity, which is different from polar ice (Gu et al., 2013). Therefore, the parameters in NEMO-Bohai need to be modified accordingly.

Figure 5 Simulated monthly mean current velocities at the surface and 16 m depth in August 2012 and February 2013. The monthly mean current velocities are calculated based on the outputs with hourly intervals. The black lines and arrows represent the streamlines and directions of the current vector field, respectively. The filled contours denote the current speed (in m s⁻¹).

In NEMO-Bohai, we selected and adjusted a series of sea ice model parameters. We increase the number of ice categories (jpl) and the number of ice layers (nlay_i) because Massonnet et al. (2019) suggested that increasing the resolution of the ice thickness distribution results in small but non-negligible variations in the ice extent and volume. The ice thickness is discretized into 10 categories (0.0, 0.0476, 0.0976, 0.150, 0.205, 0.263, 0.324, 0.388, 0.455, 0.526, 0.601 m), and the thermodynamic calculations used five ice layers (nlay_i = 5). The minimum ice thickness used in remapping is set at the minimum value (rn_himin = 0.01 m) to capture new thin ice formations better. The thickness of new ice formed in open waters (rn_hinew) is set to 0.02 m as the value must be larger than rn_hnewice. According to the test run of NEMO-Bohai, the sea ice thickness is generally overestimated. Compared to the referenced value of 2.0 m, we reset the average ice thickness (rn_himean) to 0.20 m according to the observations (Li et al., 2004; Yuan et al., 2012).

The ice initialization is activated (ln_iceini = true), the initial ice salinity (rn_smi_ini_n) is set to 7 PSU, and the initial real snow thickness (rn_hts_ini_n) is set to 0.1 m according to in situ observations (Bai and Wu, 1998; Gu et al., 2014). For the snow and ice albedos, we set them to lower values compared to the referenced ones based on the in situ observations (Zheng et al., 2017). The ice strength parameter (rn_pstar) is defined as 2.75×10⁴ N m⁻² based on previous in situ studies (Li et al., 2011). A ridging scheme is considered in SI³. Thus, the parameter rn_hstar, which adjusts the maximum thickness of ridged ice, is reduced from 100.0 to 10.0 m following the NEMO-NORDIC 1.0 configuration (Pemberton et al., 2017) to reflect the generally thinner Bohai Sea ice. In addition, the free-slip lateral boundary condition is chosen for sea ice dynamics (rn_ishlat = 0), which is synchronized with the ocean model, and a landfast parameterization is set (ln_landfast_L16 = true.) because there is landfast ice in eastern Liaodong Bay.

Figure 6 Comparison of vertical profiles of water temperature (^∘C) along transect A–B (locations shown in Fig. 1) between NEMO-Bohai (a–d) and the marine atlas of the Bohai Sea (Chen, 1992) (e–h) in February, May, August, and November.

3 Model comparison and validation

In this section, we analyse how well NEMO-Bohai reproduces ocean properties (SSH, SST, SSS, temperature and salinity stratification, currents, and volume exchanges with the Yellow Sea) and sea ice properties (area, thickness, and volume). To evaluate the model's performance, we compare our model results to in situ and satellite observations from multiple sources.

Figure 7 Similar to Fig. 6 but for salinity (PSU).

3.1 Observational data

To evaluate the SSH, the modelled tide amplitude and phase were compared to the tide tables with hourly intervals of Yantai (37.550^∘ N, 121.383^∘ E), Dalian (38.867^∘ N, 121.683^∘ E), and Jingtanggang (39.217^∘ N, 119.017^∘ E) from the National Marine Data and Information Service during 1 to 31 October 2012. In addition, we use observed in situ data from the National Marine Environmental Monitoring Center. The sampling station is at Bayuquan (39.804^∘ N, 121.456^∘ E), and the sampling period is from 1 to 31 January 2013. For the SST evaluation, we use data from the National Marine Science Data Center. The coastal stations are located at Laohutan in Dalian (38.90^∘ N, 121.70^∘ E) and Zhifudao in Yantai (37.60^∘ N, 121.40^∘ E). The observed records cover 2000–2001, 2010 (July–December), and 2011–2015. For the sea surface salinity, the observed in situ data for eight ocean stations (Wentuozi, Huludao, Zhimaowan, Qinhuangdao, Tanggu, Longkou, Penglai, Beihuangcheng) in the Bohai Sea are from Yuan et al. (2015). Specifically, the observation period in Huludao, Zhimaowan, Qinhuangdao, Tanggu, Longkou, and Beihuangcheng is from 1960 to 2013, while the period in Wentuozi and Penglai is 1987–2001 and 2012–2013, respectively. The climatological vertical profiles of temperature and salinity are from the marine atlas of the Bohai Sea based on the observed data from the 1950s to 1990s (Chen, 1992).

Figure 8 Comparison of daily sea ice area (DSIA) between NEMO-Bohai simulations (black line) and satellite-derived data (blue circles) from 1996 to 2017.

The sea ice model is evaluated against a series of observational datasets. A satellite-derived dataset covering the winters from 1988 to 2017 was retrieved. Sea ice area was extracted from two datasets: the first was based on the zonal threshold method for Advanced Very High Resolution Radiometer (AVHRR) during 1988–2000, and the second was based on an object-based feature extraction method for Moderate Resolution Imaging Spectroradiometer (MODIS) from 2001 to 2017. A detailed description can be found in Yan et al. (2017). The shapefiles derived by the aforementioned methods are further modified by visual interpretation to build a more accurate sea ice area dataset. For the sea ice thickness evaluation, we use the in situ observations from four offshore oil platforms (see Fig. 1) on 2, 6, 9, 16, and 26 January and 2, 9, 12, and 16 February during 2013 (Zeng et al., 2016; Karvonen et al., 2017).

Figure 9 Comparison of the spatial distribution of sea ice from NEMO-Bohai simulation and remote sensing inversion in freeze, severe freeze, and thaw periods for light (a–f), normal (g–l), and heavy (m–r) ice years.

3.2 Ocean

3.2.1 Sea surface height

Reliable precision of the tidal model is a prerequisite for subsequent simulations by NEMO-Bohai. The tidal model is validated by in situ hourly tidal observations at Bayuquan and tide table data at Yantai, Dalian, and Jingtanggang. Figure 3 provides a robust description of the comparison of SSH time evolution at the four stations. The tidal range in the Bohai Strait (Yantai and Dalian) is larger than that at other tide gauge stations. It is clear that the modelled water elevation at Bayuquan and Jingtanggang stations agrees less with the tide table data than at the other two stations, as shown in Table 2. Nevertheless, the model reproduced the semidiurnal M2 tidal cycle well with exact phase and amplitude at four stations with a mean correlation of 0.92, mean absolute error of 0.17 m, and root-mean-square error of 0.22 m. The model reproduces SSH standard deviation well, and the difference is within 6 cm in Yantai, Dalian, and Jingtanggang (Table 2); the largest deviation at Bayuquan is approximately 11 cm. Simulated SSH depends on the model's bathymetry, open boundary forcing, and freshwater flux (Kärnä et al., 2021). The errors are partly due to the deviation of positions between the tidal gauges and the model grid points, as the model has a horizontal grid size of ∼2 km. In addition, errors can be partly due to the inaccurate bottom topography, especially in the shallow water depth zone. As shown in Fig. 6, there is a clear topographic error between the model and the atlas for the profile section A–B. In general, the reasonable accuracy of the hydrodynamic model can partially assure rather realistic simulations of the sea ice growth and melting, as described in Sect. 3.3.

Figure 10 The comparison of sea ice thickness between the simulated and observation data in the Bohai Sea. The vertical red bars represent the range of the observations, while the blue triangles denote the simulations based on NEMO-Bohai.

Table 2 Sea surface height representation, in terms of standard deviation (SD, metres), correlation, mean absolute error (MAE, metres), and root-mean-square error (RMSE, metres), made by NEMO-Bohai for four Bohai Sea stations.

3.2.2 Sea surface temperature and salinity

As shown in Fig. 4, the modelled SST followed the seasonal cycle well, and the mean absolute error was typically less than 1 ^∘C at Laohutan and Zhifudao stations during the observed period. The modelled variations of SST are smoother than observed, and they particularly fail to capture the high fluctuations of daily variations in summer and autumn at both stations. It is also noticeable that the simulated SST at Zhifudao is generally colder in winter and warmer in summer than the observations.

Figure 11 Comparison of daily sea ice volume between satellite-derived data (circles) and NEMO-Bohai simulations (grey bars) from 1996 to 2017.

Table 3 displays the comparisons of multi-year average sea surface salinities at eight ocean stations in the Bohai Sea between simulations from NEMO-Bohai and observations reported by Yuan et al. (2015). Generally speaking, NEMO-Bohai is able to capture the main variations of the sea surface salinity in the Bohai Sea. Values at six ocean stations agree with observations with a relative deviation of less than 5 %, while the modelled values at Huludao and Tanggu are less salty than observed. This demonstrates that the model faces more significant challenges in the low-salinity areas. Primarily, the freshwater river runoff leads to lower salinity in the coastal regions. In NEMO-Bohai, the river runoff is based on climatological estimates without interannual variability. The river salinity is assumed to be 0 PSU, and the river temperature is set to the same value as the SST at the closest grid point. All these reasons might cause underestimations. In addition, the biases of river runoff and shallow water depth (generally <3 m) also need to be taken into consideration.

Figure 12 Spatial distributions of multi-year monthly average sea surface temperature in the Bohai Sea for the period of 1996–2017.

Table 3 Comparison of multi-year average sea surface salinity between the observations and simulations.

^∗ The observed sea surface salinity data is referenced from Yuan et al. (2015).

3.2.3 Current

The simulated monthly mean current velocities at the surface and 16 m depth in February and August are shown in Fig. 5. The monthly mean current velocities are calculated based on hourly model output during August 2012 and February 2013. Figure 5 shows that both the sea surface and 16 m depth currents are usually less than 0.4 m s⁻¹. Due to the blocking effects of the bays, the currents are weak at the head of the three bays, which is consistent with the observations by Chen et al. (1992). The maximum current velocity zone is located in the northern Bohai Strait, in good agreement with the model simulation result of Ji et al. (2019). The inflow and outflow occur in the northern and southern parts of the Bohai Strait in both seasons, respectively, which is consistent with the observations (Zhang et al., 2018). Specifically, the strongest modelled inflow from the Yellow Sea through the Bohai Strait occurs in a narrow channel in its northern part, namely the Laotieshan Channel, which agrees with the observations of Wan et al. (2015). In addition, Lin et al. (2011) suggested that persistent winds drive a cyclonic coastal current in the northern Yellow Sea and that one branch of the current enters the Bohai Sea at the northern Bohai Strait, which transports warm and saline water from the Yellow Sea.

Figure 13 Spatial distributions of multi-year monthly average sea surface salinity (PSU) in the Bohai Sea during 1996–2017.

The monthly mean water volume exchange at the Bohai Strait (see Fig. 1) based on hourly model simulations during August 2012 and February 2013 are also calculated to evaluate the model's performance. The Bohai Sea water exchange with the Yellow Sea is weak due to its half-closed shape and a relatively independent circulation system. The model results show that the inflow from the Yellow Sea to the Bohai Sea in August reaches 6.9×10⁴ m³ s⁻¹, almost double that in February (3.5×10⁴ m³ s⁻¹), which lies in the range of 5×10³ to 8×10⁴ m³ s⁻¹ as indicated by Bian et al. (2016). Our results show that the outflow exists in both months (August: 8.3×10⁴ m³ s⁻¹; February: 6.5×10⁴ m³ s⁻¹), with a larger amount than the inflow. The net flow appears to be outflow both in winter and summer, which is consistent with other model simulation results (Lin et al., 2002; Ji et al., 2019).

3.2.4 Vertical profile

NEMO-Bohai and observed water temperature and salinity profiles along the transect A–B (see Fig. 1) are shown in Figs. 6 and 7, respectively. Observations are from the atlas by Chen (1992), which is based on data from the 1950s to 1990s. The temporally closest 5-year period from 1995 to 2000 of NEMO-Bohai simulations was selected for model–observation comparisons. Common features are found between the model and observations. The Bohai Sea waters are mixed well vertically in autumn and winter, and they have a remarkable homogeneous vertical distribution for both temperature and salinity. In spring and summer, thermal stratification occurs with a significant cold-water core at depth that is eventually eroded in autumn. As apparent in Figs. 6 and 7, the stratification in shallow coastal waters is generally homogeneous. Similar features were reported by Wang et al. (2008), who analysed the seasonal variations of the vertical profiles in the Bohai Sea.

Figure 14 The annual average sea ice area (a), the annual maximum sea ice area (b), the simulated length of ice period (c), and the annual average sea ice volume (d) in the Bohai Sea from 1996 to 2017. The blue dots present the annual average values, and the error bars stand for one standard deviation. The dotted lines represent 22-year trends obtained through the linear fitting.

The model results, however, show some discrepancies compared to the atlas. Although the model reproduces the summer saline stratification, it is weaker than in the atlas. Nonetheless, Li et al. (2015) reported that the summer salinity stratification in the Bohai Sea is possibly weaker than in the atlas, with an observed top-to-bottom salinity difference of 0.6 PSU. The modelled salinity stratification in summer is weaker compared to the atlas, which is possibly caused by the vertical mixing setting with the used TKE closure scheme, and the high setting of vertical diffusivity in the model. In the northern part of the transect, which corresponds to northern Liaodong Bay, a negative salinity bias is visible compared to the atlas. In addition to the reasons mentioned in Sect. 3.2.2, inaccuracies in the ETOPO1 bathymetry, especially in the low-water-depth region seen from Figs. 6 and 7, may also cause these underestimations.

Figure 15 Seasonal cycles of daily sea ice area (a) and sea ice volume (b) in the Bohai Sea. The lines show the daily averages, and the shading denotes 1 standard deviation. The seasonal cycles are calculated for the periods 1996–2017.

3.3 Sea ice

Sea ice is the focus of our study and is, in this section, directly validated against observations with respect to the area, thickness, volume, and their variations.

3.3.1 Sea ice area

To evaluate the performance of NEMO-Bohai in simulating the characteristics of sea ice area, a series of variables were validated, including the daily sea ice area (DSIA), annual maximum sea ice area (AMSIA), and spatial patterns. Figure 8 shows that the model-estimated DSIA agrees well with the satellite-derived observations (r=0.85, p<0.01). In addition to the good agreement of interannual variability, the occurrence dates of the AMSIA are quite close to the observed dates from 1996 to 2017, as listed in Table 4. The percentage of the occurrence dates of the AMSIA less than 3 d accounted for 72.7 %, and the mean absolute error is 3.5 d. In 2013, the maximum sea ice area was simulated on 9 February 2013, which is only 1 d later than reported (North China Sea Branch of State Oceanic Administration, 2013). However, the date retrieved from MODIS was determined on 17 January as the data on 8 February 2013 was missing due to severe cloud cover (Yan et al., 2017). The largest difference happened in the winter of 2010 with 18 d earlier than the observations. Specifically, DSIA exhibited a bimodal distribution with two remarkable peaks in the winter of 2010. Satellite data showed that the AMSIA happened on 12 February, followed by the second maximum on 26 January, while NEMO-Bohai-simulated AMSIA was on 26 January, followed by the second maximum on 12 February. It is worth noting that the satellite shows that the largest DSIA happened in the winter of 2010, while the NEMO-Bohai-simulated DSIA reached the maximum in 2011. On the contrary, the modelled occurrence date of the AMSIA was 6 d later than the observed date (7 February) in 2008.

Table 4 Comparison of the occurrence dates of the annual maximum sea ice area between the observations and model simulations.

Figure 16 Climatological monthly mean sea ice concentration in the Bohai Sea during 1996–2017, calculated as the average of all daily sea ice concentrations for each month during the ice period.

The comparison of the spatial distribution of sea ice from the NEMO-Bohai simulation and remote sensing inversion during the freeze, severe freeze, and thaw periods in a light ice year (2007), normal ice year (2009), and heavy ice year (2010) is shown in Fig. 9. The simulated spatial distributions reflect general characteristics of sea ice evolution in the Bohai Sea except for the bias at the ice edge zones. Sea ice is mainly located in Liaodong Bay in light and normal ice years with extension to Bohai Bay and Laizhou Bay in heavy ice years. It is worth mentioning that NEMO-Bohai has reproduced the development cycle of the sea ice well but shows a relatively slow thawing process.

3.3.2 Sea ice thickness and volume

Figure 10 shows that the modelled sea ice thickness based on NEMO-Bohai mostly lies in the range of in situ observations with a slight overestimation. The mean relative bias of sea ice thickness between the simulations and observations is 4.6 cm (n=37). Other studies also revealed that sea ice thickness is overestimated by NEMO in regional studies, such as in the Baltic Sea (Tedesco et al., 2017) and the Canadian Arctic Archipelago (Hu et al., 2018). When comparing ice thickness, the discrepancy of spatial scales makes it difficult to directly compare the modelled average of grids (∼2 km) to the exact location in situ measurement as sea ice thickness varies significantly in the Bohai Sea, especially in the floating ice zone.

Figure 17 Similar to Fig. 16 but for sea ice thickness.

Sea ice volume is defined as the total ice over the whole Bohai Sea, calculated through sea ice concentration multiplied by ice thickness in all grids. As demonstrated in Fig. 11, modelled daily sea ice volume is in reasonable agreement with satellite-derived data between 1996 and 2017. The modelled sea ice volume is higher than the observed, with the mean relative bias of 15.1×10⁸ m³ (RMSE = 20.3×10⁸ m³). The correlation coefficient between simulations and satellite-derived data is relatively high (r=0.80, p<0.01, n=432), indicating reasonable agreement at daily sea ice volume variations. In 2010, NEMO-Bohai simulations reproduced the range of 0.1–73.0×10⁸ m³, which agrees with the range of 4.4–63.0×10⁸ m³ reported by Liu et al. (2013) based on AVHRR satellite images. However, it is noticeable that there is a lag in modelling sea ice formation during the early freezing period compared to satellite observations.

Figure 18 Frequency distribution of daily sea ice thickness larger than 10 cm (a) and 25 cm (b) in the Bohai Sea during 1996–2017.

Results of model comparison with in situ data and satellite-derived data confirm the robustness of the developed model, which allows us to use it in a more detailed evaluation of the spatial and temporal changes of Bohai Sea ice and study the continuous processes of freezing and melting at daily scales.

4 Initial results and applications

4.1 Ocean

4.1.1 Sea surface temperature

Figure 12 shows the monthly seasonal cycle of average SST in the Bohai Sea. There is an obvious “warm tongue” starting from the Bohai Strait to the central Bohai Sea in winter (December, January, and February) due to the warm current from the Yellow Sea. It turns into a “cold tongue” in summer due to faster temperature increases in shallow coastal waters as the influence of the warm current from the Yellow Sea also weakens. During January and February, the SST is below 0 ^∘C in all three bays. There are large areas with SSTs below −1 ^∘C in Liaodong Bay, which is the coldest sea in China during the winter (Ju and Xiong, 2013). This extremely cold water provides an ideal platform for sea ice formation.

4.1.2 Sea surface salinity

As shown in Fig. 13, except for estuaries, the salinity isolines of the Bohai Sea are roughly parallel to the coast, as the salinity in the coastal water is generally low due to river runoff. The average salinity of the whole Bohai Sea is low in summer (lowest in August at 29.5 PSU) and high in winter (highest in February at 30.0 PSU). The Bohai Strait in the east, which connects the Bohai Sea with the Yellow Sea, exhibits the highest salinity in the entire Bohai Sea, with an average of about 32.8 PSU.

Figure 19 Similar to Fig. 16 but for sea ice velocity.

4.2 Sea ice

4.2.1 Temporal variation of Bohai sea ice

The variation of annual average and maximum sea ice area in the Bohai Sea based on NEMO-Bohai from 1996 to 2017 are shown in Fig. 14a and b, respectively. The annual average sea ice area (AASIA) is an average of the DSIAs during the ice period for each year. The sea ice area exhibits apparent interannual and decadal variability in the study period. The mean AASIA during 1996–2017 is 1.07×10⁴ km², accounting for 14.3 % of the total sea area. 2010/11 and 2011/12 winters exhibited high AASIA values of 1.71×10⁴ and 1.73×10⁴ km², covering 22.8 % and 23.1 %, respectively. The minimum AASIA reaches 0.54×10⁴ km² (7.2 % coverage) in 2007 winter, only accounting for half of the multi-year average AASIA. Analogously to AASIA, AMSIA is the maximum of the DSIAs during the period for each year. The mean AMSIA during 1996–2017 is 2.07×10⁴ km² (27.6 % coverage), and the AMSIA values in 2001/02, 2010/11, 2011/12 and 2016/17 winters are exceptionally high, larger than 3.00×10⁴ km² (40 % coverage). No significant trend in AASIA was observed, while AMSIA showed a non-significant increasing trend with 0.1×10⁴ km² per decade (r=0.08, p=0.72) from 1996 to 2017. In addition, a strong positive correlation between AASIA and AMSIA (r=0.94, p<0.01) was found.

The length of the ice period is defined as days with a sea ice area greater than 100 km² from December to March. It varies from 60 to 110 d during 1996–2017, as shown in Fig. 14c. The mean length of the ice period is 87 d, with a standard deviation of 12 d, and there is a slight increasing but statistically insignificant trend (r=0.12, p=0.59). The annual average sea ice volume (AASIV, defined as the average daily sea ice volume in the ice period for each year) from 1996 to 2017 is 2.24 billion cubic metres. As shown in Fig. 14d, no significant trends can be found for AASIV during the study period. The highest value of AASIV appeared in 2010, which was 3.76 billion cubic metres, and the lowest value appeared in 2007 with 0.7 billion cubic metres, accounting for only 18.6 % of the ice volume of the highest year. The strong interannual variability of sea ice volume is determined by the strong interannual variability of sea ice area (r=0.96, p<0.01). Interestingly, it can be noticed that the ice period is not highly correlated with sea ice area or thickness, implying complicated processes of sea ice forming, growth, and thaw in the Bohai Sea.

Figure 20 Pearson correlations between annual average sea ice area (AASIA) (a–d) and annual average sea ice volume (AASIV) (e–h) and mean winter air temperature (MWAT), mean winter relative humidity (MWRH), mean winter wind speed (MWWS), and mean winter air pressure (MWAP) through detrended analysis for the period 1996–2017. The climate parameters are calculated based on the measured daily meteorological data at 12 meteorological stations surrounding the Bohai Sea during the winter (DJFM).

The climatological seasonal cycles of the sea ice area and volume in the Bohai Sea averaged over 1996–2017 show unimodal variations (Fig. 15). Ice usually starts to form in mid-December and reaches the maximum in early February. It then starts to melt until the Bohai Sea becomes free of ice by mid-to-late March. The climatological mean of the length of the ice period is about 3.5 months, and the freezing period (∼8 weeks) is longer than the thawing period (∼7 weeks). The asymmetrical process is related to changing rates of air temperature in freezing and thawing periods (Yan et al., 2020). Similarly, the seasonal variation of sea ice volume is an asymmetrical process. It is worth mentioning that the standard deviations of sea ice area and volume during the freezing period are quite large, showing that the Bohai Sea ice has a large interannual variability.

4.2.2 Spatial variation of Bohai sea ice

This section will mainly discuss the spatial features of Bohai Sea ice. As shown in Fig. 16, there are substantial seasonal and spatial differences in sea ice concentrations. To be specific, seawater freezes first in Liaodong Bay, then in Bohai Bay, and finally in Laizhou Bay, showing significant variation with latitude; the melting process happens in exactly the reverse order. Similar to the aforementioned temporal changes in the last section, the sea ice concentration reaches its peak in February. Seawater in Liaodong Bay freezes most severely, and the high concentration zone moves from the north in January to the east in February, where it survives the longest until March.

As shown from Fig. 17, sea ice thickness shows quite similar seasonal and spatial features to sea ice concentration. The monthly mean Bohai Sea ice thickness simulated by NEMO-Bohai usually reaches its maximum in February, with a monthly mean thickness of 16.9 cm, following the second-highest monthly mean thickness of 15.8 cm in January. Sea ice is thicker on the eastern coast of Liaodong Bay than that on the western coast in January and February. The maximum sea ice thickness appears near Bayuquan on the eastern coast of Liaodong Bay, where the thickness reaches up to ∼60 cm in February.

Figure 21 Spatial correlation between daily sea ice concentration in the Bohai Sea and (a) vertically integrated ocean heat content, and (b) vertically integrated ocean salt content from the surface to the bottom of the mixed layer during the ice period (DJFM).

According to the guideline for risk assessment and zoning of sea ice disaster issued by the State Oceanic Administration in 2016, a high-risk level is reached when sea ice thickness becomes greater than 25 cm, while the risk level is low when sea ice thickness is lower than 10 cm. When sea ice thickness is between 10 and 25 cm, the risk level is moderate. Accordingly, in this study, the thresholds for moderate-risk and high-risk levels of sea ice disaster are set at 10 and 25 cm, respectively. In Fig. 18, sea ice risk maps are calculated based on daily sea ice thickness from 1996 to 2017, and they clearly show that the high-risk area is mainly located at Liaodong Bay, with the highest risk area in the eastern Liaodong Bay. Thus, the seas around Yingkou, Bayuquan, and Wafangdian bear high exposure to sea ice disasters.

As the drift ice is a major component in Bohai Sea ice, studying its motion characteristics is also essential. As shown in Fig. 19, the direction is mainly northeast–southwest, and the drift exhibits a high spatial variability. The zone with the highest sea ice velocity (∼0.15 m s⁻¹ and primarily directed to the southwest) is visible at the southeast edge of Liaodong Bay next to the high-concentration and high-thickness area. In February, the velocity in the northern Liaodong Bay is significantly lower than that in the southern part near the ice edge.

4.2.3 Correlations between sea ice and climate factors

Synoptic forcing may play an essential role in the changes of Bohai Sea ice, which is primarily influenced by the Eurasian continental climate. To explore the potential regional climate drivers on the evolution of sea ice, correlations between sea ice and various indices based on detrended air temperature, air pressure, wind speed, and precipitation were examined. These indices were calculated based on daily meteorological measurements during the winter (DJFM) from 1996 to 2017 at 12 weather stations (Dalian, Wafangdian, Xiongyue, Yingkou, Jinzhou, Xingcheng, Suizhong, Qinhuangdao, Laoting, Tanggu, Huanghua, Dongying) surrounding the Bohai Sea. These measurements were obtained from the China Meteorological Data Service Center (http://data.cma.cn/en, last access: 3 March 2021). As shown in Fig. 20a, a strong negative correlation ($r=-\mathrm{0.78}$, p<0.01) between the detrended time series of AASIA and the mean winter air temperature (MWAT) is observed. Similarly, the detrended AASIV exhibits even a stronger negative correlation ($r=-\mathrm{0.87}$, p<0.01, Fig. 20e) with MWAT, suggesting that MWAT is an essential regional climate factor for sea ice changes. The detrended Bohai Sea ice area exhibits a weak but insignificant positive correlation (r=0.28, p=0.21, Fig. 20b) with the mean winter relative humidity (MWRH). The correlation between the detrended AASIV and the detrended MWRH is also found to be weak and insignificant (r=0.36, p=0.10, Fig. 20f). In addition, the correlations between the detrended AASIA and the detrended time series of mean winter wind speed (MWWS, r=0.13, p=0.56, Fig. 20c) and mean winter air pressure (MWAP, r=0.04, p=0.85, Fig. 20d) were found to be very weak and insignificant. Similar very weak and insignificant correlations were also found between AASIV and MWWS/MWAP (see Fig. 20g, h). The above correlation analysis indicates that the local temperature is the controlling factor for the sea ice evolution in the Bohai Sea.

Figure 21 illustrates the spatial correlations between daily sea ice concentration with vertically integrated ocean heat and salt content from the surface to the bottom of the mixed layer during the ice period for each grid point. Ocean heat content strongly correlates with sea ice concentration in Liaodong Bay and moderately correlates near the coastal areas. The spatial pattern corresponds quite nicely with climatological monthly mean sea ice concentration (shown in Fig. 16), with the highest values appearing in the eastern part of Liaodong Bay. On the other hand, a very weak positive or negative correlation between daily sea ice concentration and integrated ocean salt content was found in most sea areas, indicating a more complicated relationship throughout the whole Bohai Sea. Therefore, the interannual variability of sea ice is more dominated by ocean heat content than salt.

5 Conclusion and perspectives for the modelling platform

In this study, we provided a detailed description of NEMO-Bohai, a newly developed setup of an ocean–ice model for the Bohai Sea. The primary intent of our study was to test how NEMO-Bohai represents the ocean characteristics and sea ice properties.

Comparisons with observational data confirm that NEMO-Bohai is able to reproduce the ocean properties reasonably well, including ocean surface information (e.g. SSH, SST, SSS), currents, and temperature and salinity stratification. However, the ongoing development of the NEMO ocean engine, which will soon provide an updated version with the inclusion of wetting and drying processes, could improve the model performance in terms of sea surface height in shallow areas (Hordoir et al., 2019), especially for the Bohai Sea, which has an average depth of only 18 m. A higher-precision bathymetry could also improve the SSH performance. In this study, the monthly climatology of coastal runoff flux from Dai et al. (2009) was recorded only for 1948–2004, which is different from the study period. More importantly, the river runoff is based on climatological estimates without interannual variability, and we assume that the river salinity is set to 0 PSU and that the runoff temperature is the same as the ocean surface. Major future development would be using gauge records with observations of temperature and salinity assimilated into the NEMO-Bohai model, which would reduce uncertainties in temperature and salinity. Furthermore, a possible implementation of multiple embedded methods (Hvatov et al., 2019; Schwarzkopf et al., 2019), such as Global Ocean–West Pacific Ocean–East China Sea–Bohai Sea rather than the current Global Ocean–Bohai Sea direct nesting, should be investigated in the future. In addition, in order to carry out more accurate estimations of vertical mixing, it is worth implementing the experiments of turbulent vertical mixing options (Reffray et al., 2015) for the Bohai Sea for further development of NEMO-Bohai. Moreover, interactive feedback from the Bohai Sea to the global ocean could also be considered in a two-way coupling method compared to current one-way nesting.

For the sea ice component, NEMO-Bohai satisfactorily reproduced the seasonal and interannual variabilities of sea ice area compared to the satellite remote sensing data for the period 1996–2017. The modelled dates of the annual maximum sea ice area were 0.9 d earlier than the observed ones. Spatially, the simulation results realistically reflect the main characteristics of Bohai Sea ice evolution compared to the satellite data. Therefore, NEMO-Bohai can reliably be used to provide ice information during the dates without satellite-derived data. Nevertheless, it is worth mentioning that the simulated sea ice area tends to be somewhat overestimated, which is also reported for other seas in earlier NEMO-related publications (Blockley et al., 2014; Massonnet et al., 2011; Rjazin et al., 2019). In spring, the NEMO-Bohai melting process is delayed, particularly in the land-fast ice zone in eastern Liaodong Bay. It is likely that thick ice melts away slower than thinner observed ice. The applied parameterization of the surface albedo may cause a slow melting process, which does not take realistically into account the relevant physical processes, such as surface melt ponds, as the surface albedo continues to decrease until sea ice completely disappears (Mortin et al., 2016). Thus, regional atmospheric forcing data with higher accuracy and resolution can be used for further development. In addition, the space discrepancy between modelled sea ice thickness and extremely limited in situ observations makes their comparison difficult and introduces significant uncertainties.

In conclusion, NEMO-Bohai can simulate ocean and sea ice properties with reasonable skill in a broad spatiotemporal context, especially in terms of seasonal evolution and long-term interannual variations of sea ice. This finding implies that NEMO-Bohai complements the discontinuous satellite data well in sea ice hazard risk analysis. Therefore, NEMO-Bohai is a valuable tool for long-term ocean and ice simulations and climate change studies.