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 15 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 (SI3) (NEMO4.0-SI3) 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 20 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 over-estimations. NEMO-Bohai can simulate seasonal sea 25 ice evolution and long-term interannual variations. Hence, Nemo-Bohai is a valuable tool for long-term ocean and ice simulations as well as climate change studies.


Introduction
The Bohai Sea is the southernmost seasonal frozen sea in the Northern Hemisphere (Yan et al., 30 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 35 lower than that at the same latitude 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 40 could cover half of the sea area during severe winters (Su and Wang, 2012).
With four of China's top ten 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 45 vessels were damaged . About 2×10 5 hectares offshore aquatic farms were covered by sea ice, and the aquaculture industry losses reached ~6 billion yuan (~878 million US dollars) . 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 50 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 the 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), which are distributed in the Northwest Pacific 55 Ocean with the Liaodong Bay of the Bohai Sea as the southernmost geographic breeding site. At the end of October each year, the seals traveled from the Yellow Sea to Liaodong Bay. However, the population 3 of the spotted seal has experienced several drastic declines from ~2300 to ~1000 individuals over the past three 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 60 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 in mid-high latitude 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 65 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). 70 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 one-week or one-year case studies. Wang et al. (1984) reported the first sea ice dynamicthermodynamic model for the Bohai Sea, which simulated the sea ice freezing-melting cycle at a 75 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 the ice thickness distribution function consisting of three types (level ice, drift ice, and open water). Bai and Wu (1998)  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 aiming to understand the interannual or interdecadal variability.
NEMO is a state-of-the-art numerical modelling framework designed for research activities and 90 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/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).

95
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 South Indian Ocean (Schwarzkopf et al., 2019), the South China Sea (Thompson et al., 2019), and the Black Sea 100 (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 105 and assess its performance. The paper is organized as follows: Sec. 2 introduces the model and observation dataset. Model comparison and validation are carried out in Sec. 3. The analysis of sea ice variation based on the 22-year hindcast simulations of NEMO-Bohai is presented in Sec. 4. Discussion and a summary are provided in Sec. 5.

110
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 Initiative (SI 3 ) 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 the flow relaxation scheme algorithm (Engedahl, 1995). 5 The global 1/4° ocean-ice model (ORCA025) maintained by the DRAKKAR group is performed. A 115 detailed description of ORCA025 is provided by BernardBarnier 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 0.5°×0.5° 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 120 temperature drift, which impacts 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 125 system on the Sisu supercomputer and requires ~6000 CPU hours per simulated year. The model output frequency is set to every five days from 1987 to 1994 for the test run and increases to every day for the

NEMO-Bohai: The ocean model
The NEMO-Bohai domain consists of one central zone and three bays (see Fig. 1). The area covers 117°-122° E and 37°-41° N. The horizontal resolution is 0.025° in orthogonal curvilinear coordinates, 135 which is equivalent to a resolution of approximately 2 km (204 × 244 = 49 776 grid points). The bathymetry is interpolating into the target grids according to the ETOPO1 1 arc-minute 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 ten 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 140 the sea ice model every three time steps. The split-explicit time stepping method (ln_dynspg_ts = .true.) is applied to compute barotropic/baroclinic mode and their interactions.  (Engedahl, 1995), while Flather boundary 155 condition (Flather, 1976) is used for barotropic dynamics, such as SSH and barotropic velocities (u2d,

165
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 field of the 2 m air temperature, 10 m wind speed, downward short/long 170 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 −2 . 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 175 the Yellow River, Liao River, and Hai River, are acquired from JRA-55 from 1995 to 2018. The runoff forcing is considered as 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.

180
Regarding the configuration in the Gulf of Finland (rn_z0 = 5.0×10 -3 m) (Westerlund, 2019), we set rn_z0 to 5.0×10 -1 , 5.0×10 -2 , 5.0×10 -3 , 5.0×10 -4 , 5.0×10 -5 , and 5.0×10 -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 is, the more accurate the tide is. 8 However, it also needs to be combined with the actual situation at the bottom of the Bohai Sea. With 185 reference to , the bottom roughness had a magnitude of 1×10 -3 m. Therefore, the bottom roughness in NEMO-Bohai is determined to be 5.0 × 10 -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, 190 1979) is responsible for the tracer advection. Some critical parameters of the ocean part of NEMO-Bohai are listed in Table 1.  (Aksenov et al., 2019). SI 3 merges the capabilities of three sea ice models previously used in NEMO (CICE, GELATO, and LIM). SI 3 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 200 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 Formatted Table   Formatted Table   9 relatively thin and has low salinity, different from polar ice . Therefore, the parameters in NEMO-Bohai are required to be modified accordingly.

205
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) since Massonnet et al. (2019) suggested that increasing the resolution of the ice thickness distribution results in small but nonnegligible variations in the ice extent and volume. The ice thickness is discretized into ten categories ( The ice initialization is activated (ln_iceini = .true.), and 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 220 al., 2017). The ice strength parameter (rn_pstar) is defined as 2.75×10 4 N m −2 based on previous insitu studies (Li et al., 2011). A ridging scheme is considered in SI 3 . Thus, the parameter rn_hstar, which adjusts the maximum thickness of ridged ice, is reduced from 100.0 m to 10.0 m following the NEMO-NORDIC 1.0 configuration (Pemberton et al., 2017) to reflect the general thinner Bohai Sea ice.
In addition, the free-slip lateral boundary condition is chosen for sea ice dynamics (rn_ishlat = 0.), 225 which is synchronized with the ocean model, and a landfast parameterization is set (ln_landfast_L16 = .true.) because there exists landfast ice in the eastern Liaodong Bay.

Model comparison and validation
In this section, we analyze how well NEMO-Bohai reproduces ocean properties (SSH, SST, SSS, 230 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 multiple sources in-situ and satellite observations.

Observational data
To evaluate the SSH, the modeled tide amplitude and phase were compared to the tide tables with we use the in-situ observations from four offshore oil platforms (see Fig. 1) on January 2, 6, 9, 16, 26, February 2, 9, 12, 16 during 2013(Zeng et al., 2016Karvonen et al., 2017).

Sea surface temperature and salinity
As shown in Fig. 4, the modeled SST followed well the seasonal cycle, and the mean absolute error 285 was typically less than 1 °C at Laohutan and Zhifudao stations during the observed period. The modeled variations of SST are smoother than observed, and they particularly fail to capture the high fluctuations of daily variations in summer/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.   Yuan et al. (2015).
Generally speaking, NEMO-Bohai is able to capture the main variations of the sea surface salinity in the 295 Bohai Sea. Values at six ocean stations agree with observations with the relative deviation less than 5%, while the modeled 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 300 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 considerations. 14

Sea currentCurrent
The circulation in the Bohai Sea is mainly barotropic and results from the combined effects of tides and winds .

Vertical profile
NEMO-Bohai and observed water temperature and salinity profiles along the transect AB (see Fig.1) are shown in Fig. 6 and Fig. 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 345 NEMO-Bohai simulations was selected for model-observation comparisons. Common features are found both in the model and observations. The Bohai Sea waters are vertically well-mixed 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, eventually eroded in autumn. As apparent in Fig. 6 and Fig. 7, the stratification in shallow coastal waters is generally 350 homogeneous. Similar features were reported by Wang et al. (2008), who analyzed the seasonal variations of the vertical profiles in the Bohai Sea.
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 355 an observed top-to-bottom salinity difference of 0.6 PSU. The modeled salinity distribution along with transect AB duringstratification in summer is weaker compared to the atlas, which is possibly affectedcaused by the vertical mixing setting with the used TKE closure scheme, and the high setting of vertical diffusivity in the model. In the north part of the transect, which corresponds to the northern Liaodong Bay, a negative salinity bias is visible compared to the atlas. In addition to the reasons 360 mentioned in section 3.2.2, inaccuracies in the ETOPO1 bathymetry, especially in the low water depth region seen from Fig. 6 and Fig. 7, may also cause these underestimations.

Sea ice
Sea ice is the focus of our study and is, in this section, directly validated against observations with 370 respect to the area, thickness, volume, and their variations.   (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 well reproduced the 400 development cycle of the sea ice but with a relatively slow thawing process.  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)

415
The red vertical bars represent the range of the observations, while the blue triangles denote the simulations based on NEMO-Bohai.
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, modeled daily sea ice volume is in reasonable agreement with satellite-derived data between 1996 and 2017. The modeled sea 420 ice volume is higher than the observed, with the mean relative bias of 15.  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 440 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, the coldest sea in China during the winter (Ju and Xiong, 2013). This extremely cold water provides an ideal platform for sea ice formation.

475
The length of the ice period is defined as days with a sea ice area greater than 100 km 2 from December to March. It varies from 60 to 110 days during 1996-2017, as shown in Fig. 14c. The mean length of the ice period is 87 days, with a standard deviation of 12 days, 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 480 billion m 3 . 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 m 3 , and the lowest value appeared in 2007 with 0.7 billion m 3 , 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 485 with sea ice area or thickness, implying complicated processes of sea ice forming, growth, and thaw in the Bohai Sea.
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 25 the maximum in early February. Then it starts to melt until the Bohai Sea becomes ice-free by mid-to-490 late March. The climatological mean of the length of the ice period is about three and a half 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 pretty large, showing that 495 the Bohai Sea ice has a large inter-annual variability.  According to the guideline for risk assessment and zoning of sea ice disaster issued by the State 520 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 cm 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 cm 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 525 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-southwesterly, and the drift exhibits a high spatial variability. The zone with the highest sea ice velocity (~0.15 m s -1 and primarily directed to the southwest) is visible at the southeast edge of Liaodong Bay, next to the high concentration and 535 thickness area. In February, the velocity in the northern Liaodong Bay is significantly lower than that in the southern part near the ice edge.    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.

585
Comparisons with observational data confirm that NEMO-Bohai is able to reproduce reasonably well the ocean properties, including ocean surface information (e.g., SSH, SST, SSS), currents, and temperature and salinity stratification. However, the ongoing development of the NEMO ocean engine, soon providing 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 590 the Bohai Sea with an average depth of only 18 m. In addition, aA 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, different from the study period. More importantly, the river runoff is based on climatological estimates without interannual variability. In addition,, and we assume that the river salinity is set to 0 PSU and the runoff temperature is the same as the ocean surface.

595
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 current Global Ocean-Bohai Sea direct nesting, should be investigated in the future. In addition, in 600 order to carry out more accurate estimation 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 reproduced satisfactorily the seasonal and interannual 605 variabilities of sea ice area compared to the satellite remote sensing data for the period 1996-2017. The modeled dates of the annual maximum sea ice area were 0.9 days 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 610 Formatted: Font: 11 pt 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 the Eastern Liaodong Bay. It is likely that thick ice melts slower away 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 615 relevant physical processes, such as surface melt ponds, as the surface albedo continues to decrease until sea ice is completely disappeared (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 modeled sea ice thickness and extremely limited in-situ observations makes their comparison difficult and introduces significant uncertainties.

620
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 well complements the discontinuous satellite data in sea ice hazard risk analysis. Therefore, NEMO-Bohai is a valuable tool for long-term ocean and ice simulations and climate change studies.

625
Code and data availability The NEMO-Bohai is built upon the standard NEMO code (NEMO 4.0 beta, revision 10226