Improvements in the regional South China Sea Operational Oceanography Forecasting System (SCSOFSv2)

The South China Sea Operational Oceanography Forecasting System (SCSOFS), constructed and operated by the National Marine Environmental Forecasting Center of China, has been providing daily updated hydrodynamic forecasting in the SCS for the next five days since 2013. This paper presents recent comprehensive updates to the configurations of the physical model and data assimilation scheme in order to improve the forecasting skill of the SCSOFS. This paper highlights three of the most sensitive 15 updates, including the sea surface atmospheric forcing method, the discrete tracer advection scheme, and modification of the data assimilation scheme. Inter-comparison and accuracy assessment among the five versions were performed during the entire upgrading process using the OceanPredict Inter-comparison and Validation Task Team Class4 metrics. The results indicate that remarkable improvements have been made to the SCSOFSv2 with respect to the original version known as SCSOFSv1. The domain averaged 20 monthly mean root-mean-square errors of the sea surface temperature and sea level anomaly have decreased from 1.21 °C to 0.52 °C and from 21.6 cm to 8.5 cm, respectively.


25
The South China Sea (SCS) is located between 2°30′S～-23°30′N and 99°10′E-～121°50′E., It is the largest in area and the deepest in depth, a semi-closed marginal sea in the western Pacific and has the largest area and deepest depths. Its area is about 3.5 million km 2 , and its maximum depth is about 5300 m inat the central region. It is connecteds to the East China Sea by the Taiwan Strait to the northeast, to the North Pacific Ocean by the Luzon Strait to the east, and to the Java Sea by the Karimata Strait to the 30 south. Numerous islands, irregular and complex coastal boundaries, and drastic changes in bottom topography all together contribute to the extremelygreat complex distribution of the topography in the SCS.
The upper-layer basin-scale ocean circulations in the upper-layer of the SCS are mainly controlled by the East Asian Monsoon (Hellerman and Rosenstein, 1983), resultshowing a cyclonic gyre in winter and an 35 anti-cyclonic gyre in summer (Mao et al., 1999;Chu and Li, 2000). The dynamic multi-scale oceanic circulation dynamical processes inof the SCS are affected by various factors, i.e., the Kuroshio intrusion through the Luzon Strait (Nan et al., 2015;Farris and Wimbush, 1996;Liu et al., 2019), the internal waves Li et al., 2015) andor internal solitary waves (Zhang et al., 2018;Zhao and Alford, 2006;Cai et al., 2014) generated in the Luzon Strait and propagateding westward in the northern SCS, 40 the SCS throughflow as a branch fromof the Pacific Ocean to the Indian Ocean throughflow (Wei et al., 2019;Wang et al., 2011), and energetic mesoscale eddy activities Xu et al., 2019;Zhang et al., 2016;Zheng et al., 2017;Hwang and Chen, 2000;Wang et al., 2020). The multi-scale dynamical mechanisms in the SCS are too complex to understand clearly as yet,. iIt has always been a challenge to simulate or reproduce the ocean circulations, as well asnot to mention forecast the future oceanic status 45 using theby Operational Oceanography Forecasting System (OOFS).
Within the coordination and leadership of the Global Ocean Data Assimilation Experiment OceanView (GOV, https://www.godae-oceanview.org; Tonani et al., 2015;Dombrowsky et al., 2009), in the lastrecent decade or two, several regional OOFSs have been developed and operated based on the stateof-the-art community numerical ocean models forin different regions of the ocean. Tonani et al. (2015) 3 reportedsummarized that a total ofthere were 19 regional systems were running operationally in total untill 2015.  . Thise model consists of a fine-resolution (2 km) coastal model around Japan and an eddy-resolving (10 km) Western North Pacific model with one-way nesting.; tThe Chinese Global operational Oceanography Forecasting System was developed and operated based on the Regional Ocean Modelling System (ROMS, Shchepetkin and McWilliams, 2005) and NEMO by 65 the National Marine Environmental Forecasting Center, covering six6 subdomains from global to polar regions, Indian Ocean, Northwest Pacific, Yellow Sea and East China Sea (Kourafalou et al., 2015), and South China Sea (Zhu et al., 2016), with their horizontal resolutions ranging from 1/12° to 1/30°. It isshould be noted worth noting that there are considerable differences among thoese systems in many aspects, such as the model codes, area coverage, horizontal and/ vertical resolutions, and data assimilation 70 schemes, which are based onand so on, according to the user needs andor regional ocean characteristics.
In order to better satisfy the end users' needs, these OOFSs hasve been upgradeding and improveding constantly since they began operation. In general, most improvements to theof OOFSs weare implemented by increasing the horizontal or vertical grid resolution, changing the data assimilation schemes into a more sophisticated level, assimilating more diverse sources of observation data, and by 75 benefiting from the growth of high-performance computing power and global or regional observation networks. Initially, the MOVE/MRI.COM was developed based on a three-dimensional variational analysis scheme and was implemented in 2008 ., tThen, it was updated to athe fourdimensional variational analysis scheme to provide better representation of mesoscale processes (Usui et al., 2017). The Mercator Ocean International global monitoring and forecasting system had been 80 routinely operated in real time with an intermediate-resolution ofat 1/4° and 50 vertical levels since early 2001. An uUpgrading byof increasing the horizontal resolution was implemented in December 2010, to consisting of a 1/12° nested model over the Atlantic and Mediterranean. Real time daily services with a global 1/12° high-resolution eddy-resolving analysis and forecasting were delivered by an updated system, since 19 October, 2016. Moreover, Mercator Ocean International also continues to implement 85 regularly updates by increasing the system's complexity, such as expanding the geographical coverage, improving the models, and assimilating schemes, and have developinged several versions for the various milestones of the MyOcean project and the Copernicus Marine Environment Monitoring Service (Lellouche et al., 2013(Lellouche et al., , 2018. As mentioned in the literature of Zhu et al. (2016), the regional SCS Operational Oceanography 90 Forecasting System (SCSOFS, here after named it as SCSOFSv1) has been developed and routinely operated in real time since the beginning of 2013. It has continued to be upgraded by modifying the model settings in many aspects, such as the mesh distributions, surface atmospheric field forcing, and open boundary inputs, and so on, andby improving the data assimilation scheme according to the results of comparisongs and validationsg fromconducted by Zhu et al. (2016), in order to provide better services.

95
The primary purpose of this studypaper wa iss to introduceing the updates applied to SCSOFS and to determine which update had the great, but only show the highestest impact on the system. The other results offrom routine system updates andor improvements werewill not determinedbe illustrated or analyseddiscussed in detail. This paper is organized as follows. A detailed description of some general/basic updates applied to the 100 SCSOFS iswill be provided in Section 2. Some highlights and sensitive updates and their impacts ton the performance of the system are describedshown in Section 3. The Rresults of the inter-comparison and assessment ofor the different SCSOFS versions during the upgrading processes based on the 'Class 4 metrics' verification framework (Hernandez et al., 2009) are presentedwill be shown in Section 4. Section 5 contains a summary of the scientific improvements and future plans for the next step.

Physical model description, updates, and input datasets
This section describes severalsome general updates applied to the SCSOFSv1 in the last fewrecent couple years. The newly updated system is referrednamed toas SCSOFSv2 in this paperhere after. In order to isolate the contributions of each modification, different simulations were performed for the respective updates. However, some of the updates werehave been implemented directly according to model 110 experiences or theoreticaly knowledges, without standalone evaluation. The performances offrom a few integrated updates will be shown in Section 4 infor the different upgrading stages.

Figure 1: The model domain and bathymetry of SCSOFSv2
The SCSOFSv2 is still built based on ROMS, whichle whose version has been updated from v3.5 (svn 115 trunk revision 648 in 2013) to v3.7 (svn trunk revision 874 in 2017). In addition to aROMS v3.7 incorporates some changes for the model settings, which facilitating the operational running especially, besides of the major overhaul of the nonlinear, tangent linear, representor, multiple-grid nesting, and adjoint numerical kernels, ROMS v3.7 incorporates several changes to the model settings, which facilitate the operational running. 120 Firstly, we redistributed the land-sea grid mask layout to enable the systems mesh land boundary to fit the actual coastline better (Fig.1). Based on aBy comparisiong with the Fig. 1  can generated significant effects on the sea water volume of sea water transportation in the model domain, and thus, it would contributes to the better simulation of the ocean circulations.
The bathymetry ETOPO1 dataset used in SCSOFSv1, which has a 1 arc-minute grid resolution from the U.S. National Geophysical Data Center, was is replaced by the General Bathymetric Chart of the Oceans (GEBCO_2014 Grid) global continuous terrain model for ocean and land, which has ais with 30 arc-135 second spatial resolution in SCSOFSv2, from ETOPO1 data set in SCSOFSv1, which is with 1 arcminute grid resolution from U.S. National Geophysical Data Center. It wais also merged with the measured topographic data in the coastal areas along China mainland, and was adjusted with the tidal range. Then, it wais smoothed by applying a selective filter 8eight times to reduce the isolated seamounts on the deep ocean, so that the "slope parameter" r=Δh/2h is lower than thea maximum value r0=0.2 for 140 each grid (Beckmann and Haidvogel, 1993;Marchesiello et al., 2009), in order to supress the computational errors of the pressure-gradient (Shchepetkin and McWilliams, 2003). Then, the two grid stiffness ratios parameters, i.e., the slope parameter (r) and the Haney number, were changed from 0.22 and 9.78 in SCSOFSv1 to 0.17 and 13.80 in SCSOFSv2, respectively. The maximum depth was isstill set to be 6000 m still, but the minimum depth was changed from 10 m in SCSOFSv1 to 5 m in SCSOFSv2 145 (Wang, 1996). The final smoothed bathymetry is shown in Fig.1.
For the vertical terrain-following coordinate, it hwas been increased from 36 s-coordinate layers in SCSOFSv1 to 50 layers in SCSOFSv2. The transformation equation ofrom the original formulation wais also changed to an improved solution (Shchepetkin and McWilliams, 2005). The original vertical stretching function (Song and Haidvogel, 1994) wais replaced with an improved double stretching 150 function (Shchepetkin and McWilliams, 2005), to make it preserve a sufficient resolution in the upper 300 m in order to resolve the thermocline well. In this case, the thinnest layer was changeds from 0.  The net surface heat flux correction is still followsing Barnier et al. (1995Barnier et al. ( )'s (1995

method in
SCSOFSv2, but the parameter dQ/dSST, i.e., the of kinematic surface net heat flux sensitivity to sea surface temperature (SST), is calculated using the SST, sea surface atmospheric temperature, atmospheric density, wind speed, and sea level specific humidity, instead of setting a constant number of −-30 W m 2 K -1 for the entirewhole domain as in SCSOFSv1. Therefore, the parameter dQ/dSST varies 175 temporally and spatially. In additionMeanwhile, we use the infrared Advanced Very High Resolution Radiometer (AVHRR) satellite data are used in SCSOFSv2, which is an analysis constructed by combining observations from different platforms on a regular grid via optimum interpolation and is provided by the National Centeres for Environmental Information, instead of using the merged satellite's infrared sensors and microwave sensor, and the in-situ (buoy and ship) data global daily SST (MGDSST) 180 data obtained from the Office of Marine Prediction of the Japan Meteorological Agency used in the SCSOFSv1.

The North Equatorial Current (NEC) is an interior Sverdrup steady current in the subtropical North
Pacific and is located at about 10−°N -20°N., Itand usually bifurcates into two branches after encountering the western boundary along the Philippine coast toin the west of 130°E (Qiu and Chen, 2010). However, 185 the NEC is separated into two branches in the SCSOFSv1 due toaffected by the model's eastern lateral boundary setting,. iIts main branch is located at about 9.5°N −-13°N, and the other branch is located at 14.5−°N -17°N (Fig. 2a), which is clearly not in line with the actual locationsfact. The cause of thefor above result is that the Guam Island (shown in red circle in Fig. 2, located at about 13°26′N, 144°43′E) is included in SCSOFSv1, and itswhose location is too close to the eastern lateral boundary. There is a 190 sudden change in of the bathymetry from over 3500 m to below 500 m, serving as a largebig blockade to the NEC that once floweding into the model domain from the eastern lateral boundary. To resolve this problem, in SCSOFSv2, the eastern lateral boundary whas been moved westward from 145°E to 144°E to narrow the model domain and exclude the Guam Island in SCSOFSv2. It wais found that in SCSOFSv2, the simulated NEC remains as keeps the form of one main current until 130°E ,and then bifurcates into 195 the southward-flowing Mindanao Current and the northward-flowing Kuroshio in SCSOFSv2 (Fig. 2b).
AlsoIn addition, it ihas been shown that the Kuroshio current of eastern Philippine Island and the ocean circulations in theof north-eastern SCS greow stronger whenile the island of Guam wasgot removed.
ThisIt indicates that the location of the lateral open boundary is very important to the results of the model's simulation, and the results areit would be better when it isif being set to far enough away from 200 the island, especially forwhile the islands located in the major ocean circulations. For the advection schemes of the momentum, third-order upstream and fourth-order centered schemes weare used in both the horizontal and vertical directions. A Hharmonic mixing scheme wais used for both the viscosity for momentum and the diffusion for tracers in the horizontal. The Mellor-Yamada Level-2.5 vertical turbulent mixing closure scheme wais used for both the momentum and tracers. In 210 SCSOFSv2, theyAll of them in SCSOFSv2 are all set theto be same as in SCSOFSv1. Table 1 summarizes the main differencest characteristics between SCSOFSv1 and SCSOFSv2 after upgrading. The SCSOFSv2 is run usingwith a 5 s time step for the external mode, and a 150 s time step for the 215 internal mode under all of the new configurations describedmentioned above and those that willto be introduced in Section 3. The reason for the modificationying of the time step is related to the change inof the discrete schemes, which will be illustrated further in Section 3. First, aA 26-years climatology run is conducted for spinning-up at first, and followed by a hindcast run from 2005 to 2018 (Wang et al., 2012). The daily mean of the model results is archived and used for the subsequent evaluation. 220

Highlights, and sensitive updates, and their impacts
Most of the bias andor errors in the operational systems are mainly induced by several some major recurring problems, for example, external forcing, the intrinsic deficiencies of the numerical model (e.g., discrete schemes and sub-grid scale, parameterization schemes for sub-grid scale), initial errors, and the assimilation schemes. In this section, we elaborate upon the solutions to such problems that are applied 225 in SCSOFSv2, which werehas not discussedbeen mentioned in Section 2. All of these solutionsthem have significantly improved the model skills of the SCSOFS from different aspects, such as the SST, the threedimensional temperature and salinity structures, and the comprehensive simulationg skill, especially for the meso-scale processes.

230
The air--sea interactions areis one of the most essential physical processes that affect the vertical mixing and thermal structure of the upper-ocean. The air--sea fluxes mainly include the momentum flux, fresh water flux, and heat flux. The SST is an important indicator of the ocean circulation, ocean front, upwelling, and sea water mixing, and itswhose variations mainly depending on the air--sea interactions, and the ocean's thermal and dynamical factors (Bao et al., 2002). Thus, for the OOFS and ocean numerical modelling, the SST simulation and forecasting accuracy of SST is anone important metric forto evaluatinge the modelling and forecasting performance.
The accurate input of the sea surface atmospheric forcing plays a key role to excel in the performance of the model simulation of the SST. The ROMS provides two methods tof introducinge the sea surface atmospheric forcing: one is directly forcing the ocean model by providing momentum fluxes (wind stress), 240 net fresh water fluxes, net heat fluxes and shortwave radiation fluxes from the atmospheric datasets; the other is employing the COARE3.0 bulk algorithm (Fairall et al., 2003) to calculate the air--sea momentum, fresh water, and heat turbulent fluxes using the set of atmospheric variables from the atmospheric datasets, including the wind speed at 10 m above the sea surface, the mean sea level air pressure, the air temperature at 2 m above the sea surface, the air relative humidity at 2 m above the sea 245 surface, the downward longwave radiation flux, the precipitation rate, and the shortwave radiation fluxes (Large and Yeager, 2009). The calculations ofor the air--sea fluxes, sensible heat flux, latent heat flux, and longwave radiation can be referenced to Li et al. (2021). Since the SST useding in the calculation of thoese three air--sea fluxes is extracted from the ocean model, anthe increase in the of SST induces their variations in these fluxesas a result, which then leads to increaseding loss of ocean heat, and inhibitsing 250 further increases in the of SST;, and vice versa. Thus,It means that an effective negative feedback mechanism canould form between the SST and the SST-related heat fluxes. In this case, it is much easier to maintain the simulated SST at a reasonable level. The first method is employed in SCSOFSv1, and the second method, i.e., the bulk algorithm, is employed in SCSOFSv2.
In order to evaluate the performances of the different sea surface atmospheric forcing methods, we 255 conducted a special experiment by changing the method based on SCSOFSv1, which is referred to as BulkFormula here named the experiment in this paperas BulkFormula. In this experiment, we used the merged satellite SST analysis with a multi-scale optimal interpolation, called the Operational SST and Sea Ice Analysis (OSTIA) system, whicth globally coverage on a daily basis atnd a horizontal grid resolution of 1/20° (~6 km), which is and produced by the Met Office (Donlon et al., 2012), to verify the 260 results of the SCSOFS. October of, 2014, which representto stand for Wwinter, Sspring, sSummer, and Aautumn, respectively.
The SST differences weare calculated usingwith SCSOFSv1, BulkFormula, and SCSOFSv2 minussubtracts the OSTIA, respectively. It iwas found that the simulated SST weare higher than the OSTIA in all three sets of results. The difference from SCSOFSv1 is significantlypronouncedly higher than the differences from the BulkFormula and SCSOFSv2. The maximum differences mainly occur near 270 the coast (Fig.3 upper panels in Fig.3), especially for a few bays embedded into the mainland, which areis nearly impossible to resolve well usingwith 2--3 horizontal grids with at 1/30° resolution and within very shallow water depth in SCSOFSv1. This is because the sea surface atmospheric forcing data areis not accurate enough near the coast, and they provide an abnormally higher amount of heat to the ocean, resulting incausing the continuously heating of the coastal water. Thus, the simulated SST is beyond the 275 normal level in SCSOFSv1. This phenomenon can be significantly alleviated significantly by introducing the effective negative feedback mechanism between the model's SST and the air-sea heat flux usingby employing the COARE 3.0 bulk algorithm, which is employed in both the BulkFormula and SCSOFSv2  The performance of the model's skill for the annual mean SST RMSE iscan be improved by about 21% only by changing the method of sea surface atmospheric forcing method from directly forcing to the COARE 3.0 bulk algorithm due to the effective negative feedback mechanism.
However, the domain averaged RMSE of the monthly mean SST differences forom the SCSOFSv1 is 295 lower than that forom the BulkFormula in January and February, especially in the shallow region around the Taiwan Island. ThisIt indicates that the COARE 3.0 bulk algorithm is not necessarily a panacea, even with an effective negative feedback mechanism. This may be dependent on the surface forcing field data dependent, and the use of an accurate dataset forf the sea surface atmospheric forcing is more importanteffective than the selection of the forcing methodology selection . It also may 300 suffer from the complicated air--sea interactions and tidal mixing missinged in the model.

Discrete tTracers advection term discrete schemes
Spurious diapycnal mixing is one of the traditional errors in state-of-the-art atmospheric and oceanic models, especially for regionalthe terrain--following coordinate regional models, including both the continental slope and deep ocean (Marchesiello et al., 2009;Naughten et al., 2017;Barnier et al., 1998).
305 Marchesiello et al. (2009) from the GDEMv3 climatological initial fields, as well asand the simulated results from the fifth and the 11eleventh model years by using 1) the scheme combinations of the third-order upstream horizontal 325 advection, fourth-order centered vertical advection, and horizontal mixing on epi-neutral surfaces (hereafter referred to as UCI) and 2) the combination of the fourth-order Akima scheme (Shchepetkin and McWilliams, 2005) for both the horizontal and vertical advection terms and the scheme of horizontal mixing along Geopotential surfaces (constant Z) for tracers (hereafter referred to as AAG), respectively., Theand other settings are identical to those ofwith SCSOFSv2. Figure 7 shows the comparisons of the 330 time series of the domain averaged monthly mean temperature and the salinity in theat 1000 m layer simulated using the scheme combinations of the UCI in SCSOFSv1 and the AAG in SCSOFSv2, respectively. In order to lowersave computation costs, we only raun the model with the scheme combination of the UCI for over 16 years untill it reachesd a stable stateus.
The fourth-order Akima scheme is a little different from the fourth-order centered scheme because it 335 replacesby replacing the simple mid-point average with harmonic averaging in the calculation of the curvature term. Since the time stepping is done independently of the from spatial discretization in the ROMS, the Akima scheme has therepresents its advantage of reducing the spurious oscillations, which arises fromwith the non-smoothed advected fields, with respect to the fourth-order centered scheme McWilliams, 2003, 2005).

340
During the spin-up phase of the model from the initial conditions derived from GDEMv3, the temperature at 1000 m increases from the3.0-12.0℃ by initial settings of 3.0 ℃-12.0 ℃ (Fig.5a) to 3.0 ℃-17.2℃ ( Fig.5b), and the domain averaged monthly mean value quickly increases from 4.4 ℃ to 5.1 ℃ (Fig.7a) in 345 Figure 6: The same as Fig. 5, but for salinity.
January of the fifth model year;. tThe salinity at 1000 m increases from the initial settings of 34.26--34.62 by initial settings (Fig.6a) to 34.27--34.68 (Fig.6b), and the domain averaged monthly mean value increases rapidly from 34.50 to 34.54 (Fig.7b) in January of the fifth model year too. In particularEspecially, the increase ing the of domain averaged monthly mean value is almost linearly for 350 both the temperature and salinity in the first 50 months, indicating a fast rate of increase speed and strong spurious diapycnal mixing (Fig.7). Theose values are even higher in January inof the 11eleventh model year, the ranges (minimum and maximum values) reach to 3.0 ℃--17.3 ℃ for temperature (Fig.5c) and 34.26-34.73 for temperature (Fig.5c)for and salinity (Fig.6c), respecitvely. The domain averaged values are 5.3 ℃ for temperature and 34.56 for salinity (Fig. 7), respectively. The areas with increasing 355 temperature and salinity are mainly located on theat steep slopes and nearby regions, e.g., the central basin of the SCS, the Sulawesi Sea, and the equatorial Pacific Ocean. The monthly mean temperature in theat 1000 m layer from for SCSOFSv2 varies from the initial conditions of 3.0 ℃--12.0 ℃ in initial condition to 3.0 ℃--11.5 ℃ (Fig.5d);, and the domain averaged monthly mean value increases slightly from the initial value of 4.4 ℃ in initial to 4.5 ℃ (Fig.7a) in January inof the fifth model year (Fig.7a). The salinity at 1000 m varies from the 34.26-34.62 in initial 370 conditions of 34.26-34.62 to 34.24--34.63 (Fig. 6d);, and the domain averaged monthly mean value only slightly varies slightly from the initial value of 34.505 in initial to 34.509 (Fig.7b) in January inof the fifth model year (Fig. 7b). Theose values exhibitshow little variation untill January of the 11eleven th model year, the ranges are 3.0 ℃--11.3 ℃ for temperature (Fig. 5e) and 34.25--34.63 for salinity (Fig.   6e), and the domain averaged values are 4.6 ℃ for temperature and 34.52 for salinity (Fig. 7), 375 respectively. For tThe increment of the domain averaged value fors, temperature is about 0.2℃ and that for salinity is about 0.03, but theyyet remaining stable after 20 model years (Fig. 7). ThisIt is suggestsed that the spurious diapycnal mixing is significantlyhas been suppressed significantly by the AAG scheme combination, which can preserve the characteristics of the water masses in the deep ocean well. In additionMeantime, the temperature and salinity biases in the subsurface layer arehave been significantly 380 improved significantly, which will be shown in the latter part of this paper.
In addition, it wais found that the model skill for the SST hais also significantly been improved significantly while using the new AAG scheme employed in SCSOFSv2 ( Fig. 3 and Fig.4). The maximum of the monthly mean differences between the simulated SST simulated by SCSOFSv2 and OSTIA is about 3 °C --4 ℃, which is obviously smaller than the results of thefrom BulkFormula.  Secondly, we have introduced the method of computing the anomalies of the ensemble numbers used for 420 constructing the background error covariance following Lellouche et al. (2013). In SCSOFSv1, the anomalies are computed by subtracting a 10-year average from a long-term (typically 10 years) model free run snapshots with a five5-day interval for the ocean state, i.e., the sea surface height and threedimensional temperature, salinity, zonal velocity, and meridional velocity. In addition,And the ensemble is selected within a 60-day window around the target assimilation date from each year, resulting adding 425 in a total ofup to about 130 members in total (Ji et al., 2015;Zhu et al., 2016). However, in SCSOFSv2, a Hanning low-pass filter is employed to create the running mean according to Lellouche et al. (2013)

455
In SCSOFSv1, the analysis increments of the sea surface height and three-dimensional temperature, the salinity, and the zonal and meridional velocities produced by each analysis of the data assimilation weare applied to the model's initial fields at one time step. This inevitablyIt would induced a significant initial shock and spurious high-frequency oscillation into the model due to the imbalance between the increments and the model physics inevitably (Lellouche et al., 2013;Ourmiè res et al., 2006), and it 460 usually resultedcauses in a rapid growth of the forecast error and even lead to the model blowing-up after a few assimilation cycles or one or two-years period after the intermittent assimilation run. This wasIt is a threat to the stability and robustness of the OOFS. Therefore, we introduced the incremental analysis update (IAU) method (Bloom et al., 1996;Ourmiè res et al., 2006) to apply each analysis increment to the model integration as a forcing term in a gradual manner in SCSOFSv2 to diminish the negative impact.

465
In thisour case, we obtainedget the tendency term by dividing the increments bywith the total number of time steps within an assimilation cycle, as in most IAU methodologies, in order to make sure the time integral of the tendency term equalleds the analysis increment calculated by the EnOI. Once including the FGAT and IAU methods were included in the EnOI scheme, the entirewhole system's integral strategy hads to be adjusted by adding one more model integration over the assimilation time window (Lellouche et al., 2013). In SCSOFSv1, only one time model integration is needed only one time.
ThisIt means that once physical ocean model finishes a seven7-day run (does not need to archive snapshot 475 fields) and outputs a restart field., Tthe EnOI data assimilation module starts to calculate the analysis increments at the restart field time and adds it to the restart field,. tThen, the physical ocean model implementsmakes a hot-start from the updated restart field to run the seven7 days ofor the next cycle.
However, in SCSOFSv2, two times model integration areis needed twice due to the use ofconsidering 480 the FGAT and IAU methods (Fig.8). ThisIt means that the physical ocean model needs to be integrated 14 days in each assimilation cycle, to add the tendency term to the model prognostic equations due to the IAU method used during the first seven-7 days run (referred to as the a"Analysis sStage"), to output a restart field at the end of 7th day for hot starting the ocean model in the next cycle, and to output 3threehourly snapshots forecast fields during the second 7seven-days run (referred to as the "Fforecast sStage") 485 to be used in the next cycle by the FGAT method. The model outputs from the Aanalysis sStage are referred to as the "Bbest Eestimate", and those from the Fforecast sStage are referred as "the Fforecast".
The analysis increments are defined at the 3.5th day, but not at the end of the seven7th day as in SCSOFSv1., Twith the observed SLA and Argo vertical profiles data are within the seven7-day time window, and the AVHRR SST data on the 4fourth day are used by the FGAT method.

Inter-comparison and accuracy assessment
In order to demonstrateshow the improvements of the different SCSOFS sub-versions during the upgrading process, the results of the inter-comparison and assessment are presentedshown and discussed in this section, by using the GOV Inter-comparison and Validation Task Team (IV-TT) Class 4 verification framework (Hernandez et al., 2009). Class 4 metrics were originallyare used for inter-comparison and validation among different global or regional OOFSs or assimilation systems originally Hernandez et al., 2015;. TheyIt includes four metrics: the, namely, bias for assessing the consistency, the RMSE for assessing the quality or accuracy, the anomaly correlation for assessing the pattern of the variability, and the skill scores for assessing the utility of a forecast. They are calculated according to differences between the model values and the reference 500 measurements in observations space for each variable over a given period and spatial domain. The physical variables used in the Class 4 metrics are the SST, SLA, Argo profiles, surface currents, and sea ice. The rReference measurements, providing the ocean "truth", are selected as follow:, the SST data from the in-situ drifting BUOY, the SLA data from the AVISO along-track data, and the temperature and salinity data from the Argo profiles, respectively. They are assembled by GOV IV-TT participating 505 partners on a daily basis .
It is virtually impossible to exhaustively test and validate the exhaustively performances of all of the upgrades describedmentioned in Sections 2 and 3. Here, we separate the entirewhole upgrading procedure from SCSOFSv1 to SCSOFSv2 into four stages with three more sub-versions (v1.1, v1.2, and v1.3) according to the reality. By respecting to the previous version, Tthe major upgrades toin each new 510 version with respect to the previous version are listed in Table 2. In this studypaper, we used the Class 4 metrices and selected the first four physical variables, (SST, SLA, and Argo profiles), to inter-compare and assess the accuraciesy of the among different sub-versions of the SCSOFS (Table 3). Since none ofall the reference measurements data describedmentioned above 515 have not been used in these sub-versions of SCSOFS for those sub-versions without data assimilation, they are independent reference observation independent from SCSOFS, except for SCSOFSv2. The intercomparison and validation of theamong those sub-versions without data assimilation weare conducted for the model free-run results in 2013, and the inter-comparison and validation between v1.3 and v2 weare conducted in 2018 to validate the performance of the MOOAS. 520

SST
The accuracy of the SST is continuously increaseding from version v1 to v2, and thewith anomaly correlation increased from 0.52 in v1 to 0.74 in v2, i.e., a with percentage increase being 29.7% 525 improvement., The RMSE is decreaseding from 1.21℃ in v1 to 0.52℃ in v2, i.e., awith percentage increase being 57.0% improvement, for the annual mean of the entirewhole model domain averaged in 2013 (or v1.3 and v2 in 2018) (Table 3). For the versions v1, v1.1, v1.2, and v1.3, their anomaly correlation exhibitedshows significant seasonal variations, with high anomaly correlations in summer and low anomaly correlations in winter. It wais also foundindicated that the accuracy of the SST can be 530 benefited from the sea surface atmospheric forcing method, as well as the usage of more accurate observed SST data for the sea surface heat flux correction, temperature advection discrete scheme, and SST data assimilation. (c) (d)

RMSE in a 1°×1° bin for the versions (c) v1 (c) and (d) v2.(d), tThe calculations wereas performed for yearround in 2013 and 2018, respectively
The improvement of the SST due to sea surface atmospheric forcing method changeding mainly occurred 540 in summer time, exhibitshowing the same pattern as the results ofor the year 2014 in Figs. 3 and 4.
However, Butusing sea surface heat flux correction with more accurate observed SST data for the sea surface heat flux correctioncan improved accuracy of the SST simulation for the whole year-round (v1.2 in Fig. 10b). We also found that the OISST data wereis closer to the OSTIA than the MGDSST ( For the horizontal distribution of the SST RMSE, the large values weare mainly located inat the areas near the equator, coastal areas, and the northern lateral boundary, with most of the values larger than 1.5 ℃ and a maximum value of about 6.67 ℃ forin v1 (Fig. 10c). ForIn v1.3, due to the contributions of 555 all of the above described model updates, the pattern of the RMSE wais similar to that ofwith v1, i.e., basically without significant variations, but the maximum value decreaseds to 3.91 ℃ and most of the values weare less than 1.2 ℃. After applying MOOAS in v2 (Fig. 10d), only a few large RMSE values weare located onat the eastern coast of Philippine Island, with athe maximum value of 2.09 ℃, and most of the values were lower than 0.8 ℃. ThisIt indicates that the performance of the SST in SCSOFSv2 560 whas been improved significantly improved by due to all of the updates describedmentioned above.

SLA
For the entirewhole upgrading process, the accuracy of the SLA is also continuously increaseding from version v1 to v2, with the RMSE decreaseding from 21.6 cm infor v1 to 8.5 cm forin v2, i.e., with percentage increase being a 60.6% improvement, for the annual mean of the entirewhole model domain 565 averaged in 2013 (or in 2018 for v1.3 and v2) ( Table 3) As can be seen fFrom Fig. 11  andwith most of the values were less than 10 cm. It is well known that abundantplenty mesoscale eddies occur oin botheach sides of the Luzon Strait, in the northeastern SCS, and in the western Pacific Ocean (Fig. 12a)., The large SLA RMSEs in Figs. 11b and Fig. 11c indicateing that a pure physical ocean model cannot capture these meso-scale processes well without assimilating SLA data assimilated (Fig.12b).

605
However, Fig. 11d shows a significant reduction in thewith SLA RMSE, indicating that the meso-scale eddies can be represented by SCSOFSv2 due to assimilation of the with along-track SLA data, assimilated and the results are in good agreement with the satellite observations well (Fig. 12c). i.e., 1.44 -°C -1.75 ℃ for temperature and 0.13--0.14 for salinity (Table 3), even thoughif they decrease due to thewith model updates. In particularEspecially, for the vertical distribution, the RMSE can reach to larger than 3℃ for temperature and 0.3 for salinity in the thermocline and halocline, respectively, and it remainsed larger than 1 ℃ for temperature in the deep layer and 0.1 for salinity above a depth of 700 m depth (Figs. 13d and Fig.14d). This may result from the spurious diapycnal mixing caused by thedue 620 to UCI schemes combination schemeemployed. Those updates toin v1.1 and v1.2 can only slightly improved the three-dimensional temperature and salinity, and they did cannot contribute to their intrinsic improvements for, neither for surface forcing nor forthe lateral boundary conditions, with thean exception of the surface layer with depths of lessshallower than 100 m.
However, once the AAG schemes combination scheme was implemented employed in v1.3, the 625 improvements to the three-dimensional temperature and salinity weare significantobvious with respect to the first three versions (Figs. 13a,b and Fig.14a,b). The anomaly correlation increaseds to 0.38--0.48 for temperature and 0.30--0.44 for salinity, and the RMSE decreaseds to 0.96 ℃--1.03 ℃℃ for (c) (b) (a) temperature and 0.10--0.11 for salinity, respectively (Table 3). For the vertical distribution, the anomaly correlation remaineds at around 0.4 for both temperature and salinity in the entirewhole water column, 630 and it was greater than over 0.6 for temperature in the surface layer (Figs.13c and Fig.14c)., The RMSEs significantly decreased to less than 2 ℃ for temperature in the thermocline, and 0.25 for salinity in the halocline, and less than 1℃ for temperature and 0.1 for salinity in the deep layer (Figs.13d and Fig.14d).
For the horizontal distribution of the three-dimensional temperature and salinity RMSEs, the RMSE of the temperature iwas more likely to being more than >1.5 ℃ with maximum and minimum values ofbeing 635 4.45 ℃ and 0.49 ℃ (Fig. 13e), respectively; whike theand RMSE of salinity wasis greatlarger than 0.1,  Fig.14). The mean anomaly correlations has increased from 0.38 to 0.57 for temperature, and from 0.30 to 0.51 for salinity. The mean RMSEs has decreased from 0.96 ℃ to 0.67 ℃ for temperature, and from 0.11 to 0.08 for salinity (Table 3). For the vertical distributions of the anomaly correlation for 655 temperature, it wa's over >0.6 in the surface layer, wasover >0.4 above 600 m, and wasover >0.3 in the deep layer (Fig. 13c). The RMSE of the temperature wais less than 1.5 ℃ for the entirewhole vertical profile, and similar to in other versions, the maximum value iwas located inat the thermocline similar with other versions, but the error decreaseds dramatically (Fig. 13d). In contrast toUnlike temperature, the vertical anomaly correlation of the salinity didoes not show significantly improve below 200 m ment 660 in v2 with respect to v1.3 below 200m, and it was onlyshows a little slightly higher than that ofwhich in v1.3 (Fig. 14c) in above 200 m. TheS salinity RMSE iwas less than 0.25 for the entirewhole vertical profile, with the maximum value located at the surface and decreasing with depth, and decreasinge to less than 0.05 below 600 m (Fig. 14d). For the horizontal RMSE distribution ofin v2, most of the temperature RMSEs wereis greatlarger than 0.8 ℃ with maximum and minimum values ofbeing 1.96 ℃ and 0.03 ℃ (Fig. 13f), respectively; and most of the salinity RMSEs wereis greater than 0.1, with maximum and minimum values ofbeing 0.35 and 0.01 (Fig.14f), respectively, in 2018 (Fig. 14f).

Conclusions
The results of Tthis study illustrate thes major updates applied made to SCSOFSv1 in termsaspects of the physical model settings, inputs, and EnOI data assimilation scheme in the last fewrecent couple years following the recommendations of Zhu et al. (2016), such as redistributions of the land-water grid mask redistribution,; changes in the data sources ofor the bathymetry, the initial conditions, and the sea surface 675 forcing method; changing the and open boundary conditions changing to higher spatial and temporal resolutions;, shifting the eastern lateral boundary westward;, and increasing the vertical layers of the model, and so on.
The Tthree most significant updates weare highlighted in this paper. Firstly, the sea surface atmospheric forcing method hwas been changed from direct forcing to the BulkFormula to acquire an effective 680 negative feedback mechanism ofor the air--sea interactions by using the COARE3.0 bulk algorithm. The Uupgrades lead to more reasonable SST simulations withby eliminating of abnormal values, significantly decreasdropping of the maximum value of the monthly mean differences between the simulated SST and OSTIA, and decreasing theof domain averaged RMSE of the monthly mean SST from 0.99 °C --1.62 ℃ in SCSOFSv1 to 0.87 °C --1.15 ℃ in the BulkFormula run. The annual mean value decreaseds from 685 1.27 ℃ to 1.00 ℃, indicating that the performance of model's skill has improved by about 21%.
Secondly, the AAG scheme was substituted for the tracers advection term discrete scheme UCI has been substituted with AAG in order to suppress the spurious diapycnal mixing problem. After this substitution, the domain averaged monthly mean temperature in theat 1000 m layer decreaseds from 5.1 ℃ to 4.5 ℃, and thatwhich of the salinity decreaseds from 34.54 to 34.509, in January of the fifth model year, 690 respectively. Even after 20 model years, the domain averaged values of the temperature and salinity increments weare only about 0.2 ℃ and 0.03, respectively, suggesting that the AAG combination schemes combination can well preserve the characteristics of the water masses in the deep ocean preserve.
In addition, the model skill for the SST also can benefited from the AAG combination scheme,s and thecombination with annual mean domain averaged RMSE decreaseding from 1.00 ℃ to 0.77 ℃, i.e., 695 showinga 23% improvementing inrate for the performance.
Thirdly, the original EnOI method in SCSOFSv1 hwas been upgraded to thenew MOOAS by adding four new functions. The multi-source observation data (SST, SLA, and Argo profiles) werecan be simultaneously assimilated. simultaneously; The Hanning high-pass filter wais applied to the ensemble members from 10 years of free run while calculating the background error covariances to improve the 700 dynamic dependency.; The FGAT method with a 3three-hour time slot wais used to calculate the innovations; and the IAU technique with ais employed with 7seven-day time window was used to apply analyseis the increment into the model integration in a gradual manner.
Moreover, inter-comparison and accuracy assessment of theamong five versions weare conducted based on the GOV IV-TT Class 4 metrics for four physical variables, i.e., the SST, SLA, and Argo profiles.

705
The improvement in theof accuracy of the simulated SST was mainly dueattributes to the use of more accurate observed SST data source used for the sea surface heat flux correction, the use of the BulkFormula method for the sea surface atmospheric forcing, and the use of the AAG discrete temperature advection discrete scheme. The improvement of theSLA accuracy of the SLA as mainly due to thebenefits from good representations of the NEC pattern obtainedcaused by modifyicationg ofthe 710 model's eastern lateral boundary, the mean sea level air pressure correction, and the improvement of the three-dimensional temperature and salinity baroclinic structures improvement due toby using the AAG schemeemployed. The improvement of the three-dimensional temperature and salinity mainly benefiteds from the use of the AAG non-spurious diapycnal mixing combination schemes combination employed.
FinallyAt last, the remarkable improvements infor all of the above four variables are also benefited from 715 use of the MOOAS application. With respect to v1.3, for the v2 using the MOOAS, the domain averaged annual mean SST RMSE decreaseds from 0.66 ℃ to 0.52 ℃, i.e., awith percentage increase being 21.2% improvement., The SLA RMSE decreaseds from 12.9 cm to 8.5 cm, i.e., a with percentage increase being 34.1% improvement., The temperature profile's RMSE decreaseds from 0.96 ℃ to 0.67 ℃, i.e., awith percentage increase being 30.2% improvement., The salinity profile's RMSE decreaseds from 0.11 to 720 0.08, i.e., a with percentage increase being 27.3% improvement, in v2 while using MOOAS.
Although SCSOFSv2 is greatlyhas improved greatly compareding to the previous versions, some biases still exist, such as the structures of the temperature and salinity profiles in the subsurface, especially infor the thermocline and halocline. We plan to continue to improve the system in terms of both the physical model settings and the data assimilation scheme in thefor next step, including asuch as sub-grid 725 parameterization scheme for the unresolved physical processes, a vertical turbulent mixing scheme to consider wave mixing, a more accurate input and forcing data source, and assimilationg of more or new types of observations (glider or mooring three-dimensional temperature and salinity profiles, drifting buoys, in-situ velocity data from moorings) into the system.