The improvements to the regional South China Sea Operational Oceanography Forecasting System (SCSOFSv1)

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


Introduction
The South China Sea (SCS) is located between 2°30′S～23°30′N and 99°10′E～121°50′E, the largest in area and the deepest in depth, a semi-closed marginal sea in the western Pacific. Its area is about 3.5 million km 2 , and its maximum depth is about 5300 m at the central region. It connects to the East China Sea by the Taiwan Strait to the northeast, to the North Pacific Ocean by the Luzon Strait (LUS) to the 30 east, to the Java Sea by the Karimata Strait to the south. Numerous islands, irregular and complex coastal boundaries, and drastic changes in bottom topography all together contribute to the great complex distribution of topography in the SCS.
The upper-layer basin-scale ocean circulations of the SCS are mainly controlled by the East Asian Monsoon (Hellerman and Rosenstein, 1983), showing a cyclonic gyre in winter and an anti-cyclonic gyre 35 in summer (Mao et al., 1999;Chu and Li, 2000). The multi-scale oceanic circulation dynamical processes of the SCS are affected by various factors, i.e. the Kuroshio intrusion through the LUS (Nan et al., 2015;Farris and Wimbush, 1996;Liu et al., 2019), internal waves Li et al., 2015) or internal solitary waves (Zhang et al., 2018;Zhao and Alford, 2006;Cai et al., 2014) generated in the LUS and propagating in the northern SCS, the SCS throughflow as a branch of the Pacific to Indian Ocean 40 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, it has always been a challenge to simulate or reproduce the ocean circulations, not to mention forecast future oceanic status by Operational Oceanography Forecasting System (OOFS).

Within coordination and leadership of Global Ocean Data Assimilation Experiment (GODAE)
OceanView (GOV, https://www.godae-oceanview.org; Tonani et al., 2015;Dombrowsky et al., 2009), in recent decade or two, several regional OOFSs have been developed and operated based on the state-ofthe-art community numerical ocean models in different regions of the ocean. Tonani et al. (2015) summarized that there were 19 regional systems running operationally in total till 2015. implemented directly according to model experiences or theory knowledges, without standalone evaluation. The performances from a few integrated updates will be shown in Section 4 for different upgrading stages.

110
The SCSOFSv2 is still built based on ROMS, while whose version has been updated from v3.5 (svn trunk revision 648 in 2013) to v3.7 (svn trunk revision 874 in 2017). ROMS 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.
Firstly, we redistributed the land-sea grid mask layout to enable systems mesh land boundary fit the 115 actual coastline better (Fig.1). By comparing with the Fig. 1 from Zhu et al. (2016), a few sea areas had been changed from land to sea or inverse, e.g. along the coast of China mainland, the Vietnam and the

125
The bathymetry is replaced by the General Bathymetric Chart of the Oceans (GEBCO_2014 Grid) global continuous terrain model for ocean and land, which is with 30 arc-second spatial resolution in SCSOFSv2, from ETOPO1 data set in SCSOFSv1, which is with 1 arc-minute grid resolution from U.S. National Geophysical Data Center (NGDC) . It is also merged with the measured topographic data in the coastal areas along China mainland, and adjusted with the tidal range. Then it is smoothed by applying a selective 130 filter 8 times to reduce the isolated seamounts on the deep ocean, so that the "slope parameter" r=Δh/2h is lower than a maximum value r 0 =0.2 for 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, slope parameter (r) and Haney number, change from 0.22 and 9.78 in SCSOFSv1 to 0.17 and 13.80 in SCSOFSv2, respectively. The maximum 135 depth is set to be 6000m still, but the minimum depth changed from 10m in SCSOFSv1 to 5m in SCSOFSv2 (Wang, 1996). The final smoothed bathymetry is shown in Fig.1.
For the vertical terrain-following coordinate, it has been increased from 36 s-coordinate layers in SCSOFSv1 to 50 layers in SCSOFSv2. The transformation equation from the original formulation is also changed to an improved solution (Shchepetkin and McWilliams, 2005). The original vertical stretching 140 function (Song and Haidvogel, 1994) is replaced with an improved double stretching function (Shchepetkin and McWilliams, 2005), to make it preserve a sufficient resolution in the upper 300m in order to resolve the thermocline well. In this case, the thinnest layer changes from 0.16m in SCSOFSv1 to 0.09m in SCSOFSv2 near the surface.

160
The net surface heat flux correction is still following Barnier et al. (1995)'s method in SCSOFSv2, but the parameter (dQ/dSST) of kinematic surface net heat flux sensitivity to sea surface temperature (SST) is calculated using 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 whole domain as in SCSOFSv1. So the parameter dQ/dSST varies temporally and spatially. Meanwhile, we use the infrared 165 Advanced Very High Resolution Radiometer (AVHRR) satellite data in SCSOFSv2, which is an analysis constructed by combining observations from different platforms on a regular grid via optimum interpolation and provided by National Centers for Environmental Information (NCEI), instead of using the merged satellite's infrared sensors and microwave sensor, and in-situ (buoy and ship) data global daily SST (MGDSST) obtained from the Office of Marine Prediction of the Japan Meteorological 170 Agency (JMA) in SCSOFSv1.
The North Equatorial Current (NEC) is an interior Sverdrup steady current in the subtropical NP and located at about 10°N-20°N, and usually bifurcates into two branches after encountering the western boundary along the Philippine coast in the west of 130°E (Qiu and Chen, 2010). However, the NEC is separated into two branches in SCSOFSv1 affected by model eastern lateral boundary setting, its main 175 8 branch located at about 9.5°N-13°N, the other branch located at 14.5°N-17°N (Fig. 2a), which is clearly not in line with the fact. The cause for above result is that the Guam Island (shown in red circle in Fig. 2, located about 13°26′N, 144°43′E) is included in SCSOFSv1, whose location is too close to the eastern lateral boundary. It changes the bathymetry from over 3500m to below 500m suddenly, as a big block to the NEC once flowing into the model domain from eastern lateral boundary. To resolve this problem, the 180 eastern lateral boundary has been moved westward from 145°E to 144°E to narrow the model domain and exclude the Guam Island in SCSOFSv2. It is found that the simulated NEC keeps the form of one main current until 130°E, then bifurcates into the southward-flowing Mindanao Current and the northward-flowing Kuroshio in SCSOFSv2 (Fig. 2b). Also, it is shown that the Kuroshio of eastern  For the advection schemes of momentum, third-order upstream and fourth-order centered schemes are used in both horizontal and vertical, respectively. Harmonic mixing scheme is used for both viscosity for 195 momentum and diffusion for tracers in horizontal. Mellor-Yamada Level-2.5 vertical turbulent mixing closure scheme is used for both momentum and tracers. All of them in SCSOFSv2 are set to be same as in SCSOFSv1. Table 1 summarizes the main different characteristics between SCSOFSv1 and SCSOFSv2 after upgrading. The SCSOFSv2 is run with 5s time step for the external mode, and 150s for the internal mode under all new configurations mentioned above and in Section 3. The reason for modifying time step is related to the change of the discrete schemes, which will be illustrated in Section 3. A 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).

205
The daily mean of model results is archived and used for subsequent evaluation.

Highlights and sensitive updates and their impacts
Most of bias or errors in the operational systems are mainly induced by some major recurring problems, for example sea surface atmosphericexternal forcing, intrinsic deficiencies of numerical model (e.g., discrete schemes, parameterization schemes for sub-grid scale), initial errors, and the assimilation 210 schemes. In this section, we elaborate solutions, that are not mentioned in Sect.2, to such problems applied in SCSOFSv2. All of them have significantly improved the model skill of SCSOFS from different aspects, such SST, three dimensional temperature and salinity structures, comprehensive simulating skill especially for the meso-scale processes.

215
The air-sea interaction is one of the most essential physical processes that affect vertical mixing and thermal structure of the upper-ocean. The air-sea fluxes mainly include momentum flux, fresh water flux and heat flux. SST is an important indicator of ocean circulation, ocean front, upwelling and sea water mixing, whose variation mainly depending on the air-sea interaction, the ocean thermal and dynamical factors (Bao et al., 2002). Thus, for OOFS and ocean numerical modelling, simulation and forecast 220 accuracy of SST is one important metric to evaluate the modelling and forecasting performance.
The accurate input of sea surface atmospheric forcing plays a key role to excel in model simulation of SST. ROMS provides two methods to introduce sea surface atmospheric forcing: one is directly forcing ocean model by providing momentum fluxes (wind stress), net fresh water fluxes, net heat fluxes and shortwave radiation fluxes from atmospheric datasets; the other is employing the COARE3.0 bulk 225 algorithm (Fairall et al., 2003) to calculate air-sea momentum, fresh water and heat turbulent fluxes using the set of atmospheric variables from atmospheric datasets including wind speed at 10m above sea surface, mean sea level air pressure, air temperature at 2m above sea surface, air relative humidity at 2m above sea surface, downward longwave radiation flux, precipitation rate and shortwave radiation fluxes (Large and Yeager, 2009). The calculation for the air-sea fluxes, sensible heat flux, latent heat flux, and 230 longwave radiation can be referenced to Li et al. (2021). Since the SST using in the calculation of those three air-sea fluxes is extracted from ocean model, the increase of SST induces their variations of sensible heat flux, latent heat flux, and longwave radiation as a result, which then lead to increasing loss of ocean heat, and inhibiting further increase of SST, and vice versa. It means that an effective negative feedback mechanism could form between SST and SST-related heat fluxes. In this case, it is much easier to 235 maintain the simulated SST could maintain at a reasonable level. The first method is employed in SCSOFSv1, and the second, bulk algorithm, is employed in SCSOFSv2.
In order to evaluate the performances of different sea surface atmospheric forcing methods, we conduct a special experiment by changing the method based on SCSOFSv1, here named the experiment as BulkFormula. In this experiment, we use the merged satellite SST analysis with a multi-scale optimal 240 interpolation called the Operational SST and Sea Ice Analysis (OSTIA) system, which globally coverage 11 on a daily basis at a horizontal grid resolution of 1/20° (~6 km) and producvided by the Met Office (Donlon et al., 2012), to verify the results of SCSOFS.

SST (lower panels)
Figure3 shows the distributions of monthly mean SST differences in January, April, July, and October, 2014 to stand for Winter, Spring, Summer, and Autumn, respectively. SST differences are calculated with SCSOFSv1, BulkFormula, and SCSOFSv1.32 subtracts OSTIA SST, respectively. It is found that 250 the simulated SST are higher than OSTIA SST for all three sets of results. The difference from SCSOFSv1 is pronouncedly higher than the differences from BulkFormula and SCSOFSv1.32. The maximum differences mainly occur near coast ( Fig.3 upper panels), especially for a few bays embedded into the mainland which is hard to be resolved well with 2-3 horizontal grids at 1/30° resolution and in very shallow water depth in SCSOFSv1. This is because sea surface atmospheric forcing data is not 255 accurate enough near the coast, and provide abnormally more heat to ocean causing the continuously heating up of coastal water. Thus, simulated SST is beyond normal level in SCSOFSv1. This phenomenon can be alleviated significantly by introducing the effective negative feedback mechanism between model's SST and air-sea heat flux by employing the COARE 3.0 bulk algorithm, which is employed in both BulkFormula and SCSOFSv2 ( However, domain averaged RMSE of monthly mean SST differences from SCSOFSv1 is lower than that from BulkFormula in January and February, especially in the shallow region around the Taiwan Island.
It indicates that COARE 3.0 bulk algorithm does not work better than direct forcing at all time or 275 everywhere, even with effective negative feedback mechanism. This may be surface forcing field data dependent, and the accurate dataset of sea surface atmospheric forcing is more effective than the forcing methodology selection . It also may suffer from the complicated air-sea interactions and tidal mixing missed in the model.
Spurious diapycnal mixing is one of traditional errors in state-of-the-art atmospheric and oceanic model, especially for the 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). Marchesiello et al. (2009)

290
For SCSOFSv1, a third-order upstream horizontal advection scheme (hereafter referred to as U3H), a fourth-order centered vertical advection scheme (hereafter referred to as C4V), and the scheme of horizontal mixing on epi-neutral (constant density, hereafter referred to as ISO) surfaces for tracers are selected (Shchepetkin and McWilliams, 2005). We have encountered same problem with Marchesiello et al. (2009) for temperature ( Fig.5b and 5c) and salinity ( Fig.6b and 6c) in deep layer. Figure 5 and 6 300 show the distributions of monthly mean temperature and salinity at 1000m layer in January from GDEMv3 climatological initial fields, and the simulated results from the fifth and the eleventh model years by using the scheme combination of third-order upstream horizontal advectionU3H, fourth-order centered vertical advectionC4V and horizontal mixing on epi-neutral surfacesISO (hereafter referred to as UCI) and the combination of the fourth-order Akima scheme (Shchepetkin and McWilliams, 2005) 305 for both horizontal and vertical advection terms and the scheme of horizontal mixing along Geopotential surfaces (constant Z) for tracers (hereafter referred to as AAG), respectively, and other settings are identical with SCSOFSv2. Figure 7 shows the comparisons of time series of domain averaged monthly mean temperature and salinity at 1000 m layer simulated using the scheme combinations of UCI in SCSOFSv1 and AAG in SCSOFSv2, respectively. In order to save computation costs, we only run the 310 model with scheme combination of UCI for over 16 years till it reaches stable status.
The fourth-order Akima scheme is a little different from the fourth-order centered scheme by replacing the simple mid-point average with harmonic averaging in the calculation of curvature term. Since the time stepping is done independently from spatial discretization in ROMS, the Akima scheme represents its advantage of reducing spurious oscillations, which arises with nonsmoothed advected fields, with 315 respect to the fourth-order centered McWilliams, 2003, 2005).
During the spin-up phase of the model from the initial conditions derived from GDEMV3, the temperature at 1000 m increases from 3.0-12.0℃ in initial (Fig.5a) to 3.0-17.2℃ (Fig.5b), and the domain averaged monthly mean value quickly increases from 4.4 ℃ in initial to 5.1 ℃ (Fig.7a) in indicating a fast increasing speed and strong spurious diapycnal mixing (Fig.7). Those values are even higher in January of the eleventh model year, the ranges (minimum and maximum value) reach to 3.0- 17.3℃ for temperature (Fig.5c) and 34.26-34.73 for salinity (Fig.6c). The domain averaged values are 5.3 ℃ for temperature and 34.56 for salinity (Fig.7), respectively. The areas with increasing temperature and salinity are mainly located at steep slopes and nearby regions, e.g. the central basin of SCS, the 330 Sulawesi Sea and the equatorial Pacific Ocean. To fix this problem, we tested various model settings and compiling options available in ROMS, such as 335 increasing the number of vertical levels, changing the advection and diffusion schemes, horizontal mixing surfaces for tracers, horizontal mixing schemes. Details of how tested model settings effect on the spurious diapycnal mixing are beyond the scope of this paper, which will be discussed in a separate paper. Based on test results, we conclude that the spurious diapycnal mixing problem can be suppressed significantly by employing AAG scheme combination (Fig.5d, e and Fig.6d, e) in SCSOFSv2.

340
The monthly mean temperature at 1000 m layer from SCSOFSv2 varies from 3.0-12.0℃ in initial condition to 3.0-11.5℃ (Fig.5d), and the domain averaged monthly mean value increases slightly from 4.4 ℃ in initial to 4.5 ℃ (Fig.7a) in January of the fifth model year. The salinity at 1000 m varies from 34.26-34.62 in initial condition to 34.24-34.63 (Fig.6d), and the domain averaged monthly mean value only slightly varies from 34.505 in initial to 34.509 (Fig.7b) in January of the fifth model year. Those 345 values show little variation till January of the eleventh 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), respectively. For the increment of domain averaged values, temperature is about 0.2℃ and salinity is about 0.03, yet remaining stable after 20 model years (Fig.7).
It is suggested that spurious diapycnal mixing has been suppressed significantly by AAG scheme 350 combination, which can preserve the characteristics of water masses in deep ocean well. Meantime, the temperature and salinity bias in the subsurface layer have been improved significantly, which will be shown in later of this paper.
In addition, it is found that the model skill for SST has also been improved significantly while the new AAG scheme employed in SCSOFSv1.32 ( Fig.3 and Fig.4)

Data assimilation scheme
As mentioned as Zhu et al. (2016)   we have divided the 7-day time window into 56 3-hour time slots (Fig.8), and archived 57 snapshots with 415 a 3-hour interval while model forecast run following the previous analysis run. Then the innovations can be calculated within each 3-hour time slot by using the observations subtracts the nearest model forecast.
It means that the maximum temporal misfit of the innovations between the observation and model forecast would be decreased from 7 days to 1.5 hours by using FGAT. Meanwhile, the localization is still used with the radius set to be 150 km as in SCSOFSv1.

420
In SCSOFSv1, the analysis increments of sea surface height and three-dimensional temperature, salinity, zonal and meridional velocities produced by each analysis of data assimilation are applied to the model initial fields at one time step. It would induce model significant initial shock and spurious high-frequency oscillation due to the imbalance between the increments and the model physics inevitably (Lellouche et al., 2013;Ourmières et al., 2006), and usually causes a rapid growth of forecast error and even model 425 blow up after a few assimilation cycles or one or two-year period after the intermittent assimilation run.
It is a threat to the stability and robustness of 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.
In our case, we get the tendency term by dividing the increments with the total number of time steps 430 within an assimilation cycle as in most IAU methodologies, in order to make sure the time integral of tendency term equals the analysis increment calculated by EnOI.

435
Once including the FGAT and IAU methods in EnOI scheme, the whole system integral strategy has 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. It means that once physical ocean model finishes 7 days run (does not need to archive snapshot fields) and outputs a restart field, the EnOI data assimilation module starts to calculate the analysis increments at the restart field time and adds it to 440 the restart field, then the physical ocean model makes a hot-start from the updated restart field to run 7 days for next cycle.
However, in SCSOFSv2, two times model integration is needed due to considering the FGAT and IAU methods (Fig.8). It means that physical ocean model needs to be integrated 14 days in each assimilation 445 cycle, to add the tendency term to the model prognostic equations due to the IAU method used during the first 7 days run (referred to as "Analysis Stage"), to output restart field at the end of 7th day for hot starting ocean model in next cycle, and to output 3-hourly snapshots forecast fields during the second 7 days run (referred to as "Forecast Stage") to be used in next cycle by FGAT method. The model outputs from the Analysis Stage are referred to as "Best Estimate", and from the Forecast Stage are referred to as "Forecast". The analysis increments are defined at the 3.5th day, but not at the end of 7th day as in SCSOFSv1, with the observed SLA and Argo T/S vertical profiles data within the 7-day time window and AVHRR SST data on the 4th day used by FGAT method.

Inter-comparison and accuracy assessment
In order to show the improvements of different SCSOFS sub-versions during the upgrading process, the 455 results of inter-comparison and assessment are shown 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 are used for inter-comparison and validation among different global or regional OOFSs or assimilation systems originally Hernandez et al., 2015;Divakaran et al., 2015). It includes four metrics, namely, bias for consistency, RMSE for quality or accuracy, Anomaly Correlation 460 (AC) for pattern of the variability and skill scores for the utility of a forecast. They are calculated according to differences between model values and reference measurements in observations space for each variable over a given period and spatial domain. The physical variables used in Class 4 metrics are SST, SLA, T/S profiles, surface currents and sea ice. Reference measurements, providing ocean "truth", are selected as follow, SST from in-situ drifting BUOY, SLA from AVISO along-track data, temperature 465 and salinity from Argo profiles, respectively. They are assembled by GOV IV-TT participating partners on a daily basis .
It is virtually impossible to test and validate exhaustively performances of all upgrades mentioned in Section 2 and 3. Here, we separate the whole upgrading procedure from SCSOFSv1 to SCSOFSv2 into four stages with three more sub-versions (v1.1, v1.2, v1.3) according to the reality. By respecting to the 470 previous version, the major upgrades in each new version are listed in Table21. In this paper, we use Class 4 metrices and select the first four physical variables, SST, SLA and T/S, to inter-compare and assess the accuracy among different sub-versions of SCSOFS (Table 32). Since all the reference measurements data mentioned above have not been used in SCSOFS for those sub-versions 475 without data assimilation, they are independent reference observation from SCSOFS except for SCSOFSv2. The inter-comparison and validation among those sub-versions without data assimilation are conducted for the model free-run results in 2013, and between v1.3 and v2 are conducted in 2018 to validate the performance of MOOAS.
heat flux correction with more accurate observed SST data can improve accuracy of SST simulation for the whole year (v1.2 in Fig.10b). We also found that OISST data is closer to OSTIA than MGDSST (figure not shown). Due to benefit from those changes, the maximum and minimum value of SST RMSE have decreased from 1.92℃ and 0.71℃ in v1 to 1.52℃ and 0.60℃ in v1.2 for the whole year 2013, 500 respectively. It is worth mentioning that AAG schemes combination not only improves the deep layer temperature, but also contributes to the improvement of SST due to internal baroclinic vertical heat transport. The maximum and minimum value of SST RMSE is 1.21℃ and 0.52℃ in v1.3. For the results with data assimilation in v2, the maximum and minimum value of SST RMSE is only 1.13℃ and 0.32℃, respectively. It is better than the result in v1.3 year-round.

505
For the horizontal distribution of SST RMSE, large values are mainly located at the areas near equator, coast areas and northern lateral boundary, with most of values larger than 1.5℃ and maximum value about 6.67℃ in v1 (Fig.10c). In v1.3, due to the contributions of all the above model updates, the pattern of RMSE is similar with v1 basically without significant variations, but the maximum value decreases to 3.91℃ and most of values are less than 1.2℃. After applying MOOAS in v2 (Fig.10d), only a few large 510 RMSE values are located at the eastern coast of Philippine island with the maximum value of 2.09℃ and most of values lower than 0.8℃. It indicates that the performance of SST in SCSOFSv2 has been improved significantly due to all the updates mentioned above.

SLA
For the whole improving process, the accuracy of SLA is also continuously increasing from version v1 515 to v2, with RMSE decreasing from 21.6cm in v1 to 8.5cm in v2 with PI being 60.6%, for the annual mean of whole model domain averaged in 2013 (or in 2018 for v1.3 and v2) (

540
For the horizontal distribution of SLA RMSE, large values over 20 cm are mainly located in the area of NEC pathway, continental shelf of the northeastern SCS and northeast of LUS, with maximum value of 32.7cm in v1 (Fig.11b). In v1.3 (Fig.11c) northeastern SCS and western Pacific (Fig. 12a), large SLA RMSE in Fig. 11b and Fig. 11c indicating that pure physical ocean model cannot capture meso-scale process well without SLA data assimilated 550 (Fig.12b). However, Fig. 11d shows significant reduction with SLA RMSE, indicating that meso-scale eddies can be represented by SCSOFSv2 with along-track SLA data assimilated and agreement with satellite observations well (Fig. 12c).

T/S profiles
For 3D T/S distribution, by comparing model results with climatology T/S profiles, the results from first three versions show poor correlation with observations ( Fig.13a and Fig.14a) and large RMSE ( Fig.13b and Fig.14b), i.e. 1.44-1.75℃ for T and 0.13-0.14 for S (Table 32), even if they decrease with model 560 updates. Especially, for the vertical distribution, the RMSE can reach to larger than 3℃ for T and 0.3 for S in thermocline and halocline, respectively, and remained larger than 1℃ for T in deep layer and 0.1 for S above 700m depth ( Fig.13d and Fig.14d). This may result from spurious diapycnal mixing due to UCI schemes combination employed. Those updates in v1.1 and v1.2 can only slightly improve 3D T/S, and cannot contribute to their intrinsic improvements, neither for surface forcing nor for lateral boundary 565 conditions, except for surface layer of shallower than 100m.
However, once AAG schemes combination employed in v1.3, the improvements to 3D T/S are obvious with respect to the first three versions (Fig.13a,b and Fig.14a,b). The AC increases to 0.38-0.48 for T and 0.30-0.44 for S, and RMSE decreases to 0.96-1.03℃ for T and 0.10-0.11 for S, respectively (Table   32). For the vertical distribution, the AC remains around 0.4 for both T and S in the whole water column, 570 and over 0.6 for T in the surface layer ( Fig.13c and Fig.14c), RMSEs significantly decrease to less than 2℃ for T in thermocline and 0.25 for S in halocline, and less than 1℃ for T and 0.1 for S in deep layer ( Fig.13d and Fig.14d).
For the horizontal distribution of 3D T/S RMSE, RMSE of T is more likely being more than 1.5℃ with maximum and minimum values being 4.45℃ and 0.49℃ (Fig.13e), and RMSE of S is larger than 0.1, 575 with By employing MOOAS, accuracy of 3D T/S has been improved continuously in v2 compare to v1.3 for all the metrics in 2018 ( Fig.13 and Fig.14). The mean AC has increased from 0.38 to 0.57 for T, and from 0.30 to 0.51 for S. The mean RMSE has decreased from 0.96℃ to 0.67℃ for T, and from 0.11 to 590 0.08 for S (Table 32). For vertical distribution of AC for T, it's over 0.6 in surface, over 0.4 above 600m, and over 0.3 in deep layer (Fig.13c). RMSE of T is less than 1.5℃ for the whole vertical profile, and the maximum value is located at the thermocline similar with other versions, but the error decreases dramatically (Fig.13d). Unlike T, vertical AC of S does not show significant improvement in v2 with respect to v1.3 below 200m, and it shows a little higher than which in v1.3 (Fig.14c) in above 200m. S 595 RMSE is less than 0.25 for the whole vertical profile, with the maximum value located at surface and decreasing with depth, and decrease to less than 0.05 below 600m (Fig.14d). For the horizontal RMSE distribution in v2, most of T RMSE is larger than 0.8℃ with maximum and 600 minimum values being 1.96℃ and 0.03℃ (Fig.13f); and most of S RMSE is greater than 0.1, with maximum and minimum values being 0.35 and 0.01 (Fig.14f), respectively, in 2018.

Conclusions
This study illustrates major updates applied on SCSOFSv1 in both physical model settings, inputs and EnOI data assimilation scheme in recent couple years following the recommendations of Zhu et al. (2016), 605 such as land-water grid mask redistribution, data sources for bathymetry, initial condition, sea surface forcing and open boundary condition changing to higher spatial and temporal resolution, moving the eastern lateral boundary westward, increasing vertical layers of model, and so on.
Three most significant updates are highlighted in this paper. Firstly, sea surface atmospheric forcing method has been changed from direct forcing to BulkFormula to acquire effective negative feedback 610 mechanism of air-sea interaction by using COARE3.0 bulk algorithm. Upgrades lead to more reasonable SST simulation with eliminating of abnormal values, significantly dropping of the maximum value of monthly mean differences between simulated SST and OSTIA, and decreasing of domain averaged Secondly, tracers advection term discrete scheme UCI has been substituted with AAG in order to suppress spurious diapycnal mixing problem. After this substitution, the domain averaged monthly mean temperature at 1000m layer decreases from 5.1℃ to 4.5℃, and which of salinity decreases from 34.54 to 34.509, in January of the fifth model year, respectively. Even after 20 model years, domain averaged 620 values of temperature and salinity increments are about 0.2℃ and 0.03, suggesting that AAG schemes combination can well preserve the characteristics of water masses in deep ocean. In addition, model skill for SST also can benefit from AAG schemes combination with annual mean domain averaged RMSE decreasing from 1.00℃ to 0.77℃, showing 23% improving rate for the performance.
Thirdly, the original EnOI method in SCSOFSv1 has been upgraded to new MOOAS by adding four new 625 functions. The multi-source observation data (SST, SLA, and Argo T/S profiles) can be assimilated simultaneously; Hanning high-pass filter is applied to ensemble members from 10 years free run while calculating the background error covariances to improve the dynamic dependency; FGAT method with 3-hour time slot is used to calculate the innovations; and IAU technique is employed with 7-day time window to apply analysis increment to the model integration in a gradual manner.

630
Moreover, inter-comparison and accuracy assessment among five versions are conducted based on GOV IV-TT Class 4 metrics for four physical variables, SST, SLA, and T/S profiles. The improvement of accuracy of simulated SST mainly attributes to more accurate observed SST data source used for sea surface heat flux correction, BulkFormula method for sea surface atmospheric forcing, AAG temperature advection discrete scheme. The improvement of SLA accuracy mainly benefits from good 635 representations of NEC pattern caused by modification of model eastern lateral boundary, mean sea level air pressure correction, and T/S baroclinic structures improvement due to AAG employed. The improvement of 3D T/S mainly benefits from AAG non-spurious diapycnal mixing schemes combination employed.
application. With respect to v1.3, domain averaged annual mean SST RMSE decreases from 0.66℃ to 0.52℃ with PI being 21.2%, SLA RMSE decreases from 12.9cm to 8.5cm with PI being 34.1%, T profile RMSE decreases from 0.96℃ to 0.67℃ with PI being 30.2%, S profile RMSE decreases from 0.11 to 0.08 with PI being 27.3%, in v2 while using MOOAS.
Although SCSOFSv2 has improved greatly comparing to the previous versions, some biases still exist, 645 such as in surface and the structures of temperature and salinity in subsurface, especially for the thermocline and halocline. We plan to continue to improve the system in both physical model settings and data assimilation scheme for next step, such as sub-grid parameterization scheme for unresolved physical processes, vertical turbulent mixing scheme to consider wave mixing, more accurate input and forcing data source, and assimilating more or new types of observations (glider or mooring T/S, drifting 650 buoys, in-situ velocity from moorings) into the system.
Code and Data availability. The latest version of the source code for EnOI and ROMS trunk used to producethe results in this paper can be accessed via https://doi.org/10.5281/zenodo.5215783.
Author Contributions. XZ performed the physical model improvement and free-run simulations, designed and wrote the paper. XZ and ZZ updated MOOAS and performed the data assimilation 665 simulations. SR and AL analysed and assessed model results. SR, HW and YZ helped in reading and commenting on the paper. MZ helped in polishing the paper.
Competing interests. The authors declare that they have no conflict of interest.