1 INTRODUCTION

Kuroshio is the strongest western boundary current (WBC) in the Pacific. When passing by the Luzon Strait (LS), most of the Kuroshio water flows out of the LS through a looping path, while some of it intrudes into the South China Sea (SCS; Yuan et al., 2006; Chen et al., 2011a; Nan et al., 2011a; Shu et al., 2018a), and leads to complex circulations (Qu et al., 2004; Nan et al., 2011a; Wu and Hsin 2012; Shu et al., 2016a, 2018b) and energetic eddy activities in the northeastern SCS (Fig. 1; Wang et al., 2008a; Chen et al., 2011b; Nan et al., 2011b; Liu et al., 2012b; Shu et al., 2016b; Chu et al., 2017). The spatial patterns and temporal variability of the circulations and mesoscale processes are crucial to local marine environment and air-sea interactions, thus becoming one of the hot research topics (Xian et al., 2012; Zhang et al., 2013; Guan et al., 2014; Huang et al., 2014). However, our knowledge of the variability in the northeastern SCS is still limited due to the lack of observations. Recently, a large number of observations have been collected in the northeastern SCS by different researchers (Tian et al., 2006; Chen et al., 2010, 2011a, c; Yang et al., 2010a, 2015); however, continuous observations are rare. It is worth noting that some researchers have carried out continuous observations in the northeastern SCS. Ocean University of China designed and placed 18 mooring arrays in the northeastern SCS, and parts of station data were used in several studies (Fig. 2; Zhang et al., 2013; Guan et al., 2014; Huang et al., 2014). The mooring array has successfully documented a pair of deep-penetrating mesoscale eddies that generated southwest of Taiwan, China, and could greatly influence the deep circulation in the northeastern SCS (Zhang et al., 2013). Continuous observations always have two goals: one is to capture the potential physical phenomena in the ocean; another is to assimilate the observation results to improve the model simulation. Before assimilating the observation results to the model, designing a suitable observing system is necessary. However, the design of such an observation system as shown in Fig. 2 is based typically on empirical knowledge, rather than quantitative assessments. For the high costs involved in building and maintaining the observing system, designing and optimizing an observing system are necessary; and an observation system simulation experiment (OSSE) is a possible solution to this problem. An OSSE can be used to determine the potential benefit of future observing systems using an existing monitoring system; it help define optimal characteristics of future instruments and evaluate the present sampling distribution, and optimize construction design, and shall play a key role in developing an ocean observing system (Smith, 1993; Sakov and Oke, 2008).

Based on the data assimilation theory, the OSSE, which was originally used in meteorology (Charney et al., 1969) had been widely employed to design observing systems for the world's oceans (Hackert et al., 1998; Bishop et al., 2001; Hirschi et al., 2003; Schiller et al., 2004). Especially, in the tropical Indian Ocean (TIO), there were many applications of the OSSE to assess the ability of the monitoring system (Ballabrera-Poy et al., 2007; Oke and Schiller, 2007; Vecchi and Harrison, 2007; Sakov and Oke, 2008). Some applied the OSSE to place and optimize stations in coastal waters (Frolov et al., 2008; Yildirim et al., 2009; Xue et al., 2011, 2012). Others focused on designing the optimal monitoring sites for biological oceanography research (Lin et al., 2010; Xue et al., 2011, 2012). Varieties of methods have been introduced to the framework of OSSE, which extended the efficiency of OSSE in assessing the ability of observation systems. Using an empirical orthogonal function (EOF) and data assimilation method, a mooring array had been designed and evaluated to monitor the intraseasonal and interannual variability in the TIO (Ballabrera-Poy et al., 2007; Oke and Schiller, 2007; Vecchi and Harrison, 2007; Sakov and Oke, 2008). EOF-based approach was also employed by Yang et al. (2010b) to select the optimal sensor locations for noisy ocean measurements in the Nantucket Sound. A variance quadtree (VQT) algorithm for optimizing survey design was developed by McBratney et al. (1999) and improved by Minasny et al. (2007). Lin et al. (2010) used the VQT algorithm to have optimized plankton survey strategies in the Gulf of Maine. Based on a high-resolution ocean model, some researchers carried out the OSSE using the ensemble Kalman Filter (EnKF) data assimilation approach (Evensen, 2003, 2004; Chen et al., 2009; Xue et al., 2011, 2012; Liu et al., 2018a). According to the Kalman Filter theory, a simplified and computationally efficient, ensemble-based method for OSSE (hereafter EnOSSE) was introduced by Bishop et al. (2001), and had been used intensively for optimal array design in the last few years (Khare and Anderson, 2006; Sakov and Oke, 2008; Ye and Wang, 2011; Wang and Ye, 2012; Liu et al., 2018b). This method is used to design an optimal array based on an ensemble system state, which may come from model results or from long time series of field observations. Sakov and Oke (2008) used the long time series of sea level anomaly (SLA) data from satellite altimeter to design an optimal array to monitor the intraseasonal to interannual variability in the TIO.

There also exists intraseasonal to interannual variations in the northeastern SCS. The Kuroshio intrusion has significant intraseasonal variability (10–30 days and 60 days), the 10–30 days variability is caused by the high-frequency variability of local winds and the baroclinic instability, whereas the 60 days cycle is mainly for the intrusion of eddies (Zhang et al., 2015). The Kuroshio path in the LS also has seasonal variability for the mesoscale eddies and monsoons (Shaw, 1991; Yuan and Li, 2008; Nan et al., 2011a; Yuan and Wang, 2011). The interannual variability of Kuroshio intrusion into the SCS is mainly caused by the El Niño-Southern Oscillation (ENSO; Qu et al., 2004; Liu et al., 2006, 2012a). Some researchers focused on the intraseasonal variability in upwelling intensity in the northern SCS (Shu et al., 2011; Wang et al., 2012; Gan et al., 2015); the upwelling in the northern SCS also has strong interannual variability (Jing et al., 2011; Shu et al., 2018b). About the mesoscale eddies, Chen et al. (2011b) pointed out that the number, radius and lifetime of the eddies have seasonal and interannual variability in the northern SCS. The Luzon cyclonic gyre (LCG) also has strong seasonal variability for the monsoons (Fig. 1; Fang et al., 1998; Gan et al., 2006). Controlled by the monsoons, Kuroshio intrusion and other factors, the circulation and mesoscale eddies in the northeastern SCS present a complex intraseasonal to interannual variations (Su, 2005). Therefore, it is necessary to design an array to monitor the intraseasonal to interannual variations in the northeastern SCS.

Although OSSEs have been used widely for promoting monitoring efficiency in the world's oceans, few applications have been conducted to design and assess a stationary observation network in the northeastern SCS. Therefore, detailed studies for quantitatively examining the present stations' distribution and constructing an optimal observation system are necessary. Furthermore, previous studies usually focused on optimal station distribution in horizontal direction, while the optimal vertical distribution of observation instruments, i.e., which layers play more key roles to the observing system, was often ignored. However, when building and maintaining a stationary mooring array to monitor the current, we always need to synchronously place dozens of observation instruments and monitor units (such as conductivity-temperature-depth (CTD) machine) in different layers. To promote monitoring efficiency and to be more economical, studying optimal vertical layers is also required.

In this study, we first focus on the design of an optimal array and quantitative evaluation of the present in-situ stations using the ensemble-based method for OSSE in the northeastern SCS. Then, we quantitatively examine the performance of observations in different layers and find the optimal design in the vertical. This paper is organized as follows. In Section 2, we briefly describe the method. In Section 3, we introduce the data to be analyzed. Our results are presented in Section 4, and the summary is given in Section 5.

2 METHOD

The EnOSSE (Bishop et al., 2001; Sakov and Oke, 2008) is applied in the present study. Assuming that the uncertainty of a state vector X^n×l can be represented by a background error covariance matrix P^b, which is expressed as follows:

where A is the anomaly field of the state vector X^n×l and can be defined as , of which X is the ensemble mean of state vector X^n×l, n is the number of variables in a state vector and m is the number of samples.

According to the Kalman Filter theory (Evensen, 2003, 2004), after assimilating a set of observations, the background error covariance matrix P^b is updated, which can be represented by

where P^a is the updated P^b, K is the Kalman Gain defined as follows:

R in Eq.3 is the observational error covariance. H is an observation matrix that projects the model data onto the observational points, if assimilating q observations, H is given as H^q×n, when q=1, , k is the site of the observational point to be assimilated. The diagonal elements of matrix P^b represent the variance of ensemble matrix X, so the trace of P^b denotes the sum of the variance of X. After assimilating a set of observations, the optimal set of observation matrix H^opt is expected to minimize the trace of P^a or to maximize the trace of KHP^b in Eq.2.

According to the relation trace(AB)=trace(BA), the solution of Eq.4 can be described by

When searching for the first optimal observational site, H in Eq.5 becomes 1×n vector, replacing P^b in Eq.5 with that in Eq.1. Then, we get the first optimal observational site k, which can be written as follows:

where A_i is the i^th row of A. After finding the first optimal observation, we update P^b with Eq.2, and then calculate the next optimal observation; and P^b is updated again, until all q optimal observations are found. When n≫m, directly manipulating the matrix P^b with a dimension n×n in Eq.2 is very difficult, so a computationally effective solution is given by Bishop et al. (2001), replacing directly updating P^b in Eq.2 by updating the background ensemble A^b as follows:

In the above equations, I is the identity matrix, is symmetric and can be solved by singular value decomposition (SVD):

where U is an orthogonal matrix and Λ is a diagonal matrix. Then, T can be written as

The above method tells us how to find an optimal array, or how to get the optimal set of observation matrix H^opt. When evaluating the present in-situ stations (Fig. 2), the observation matrix H is known, and then the background ensemble A^b can be updated using Eqs.7–10, which leads to the updating of P^b in Eq.1. According to the Kalman Filter theory, the background error covariance matrix P^b characterizes the uncertainty of the initial ensemble, so we can assess the ability of the present observing system according to the updating of P^b.

3 DATA

The present in-situ stations were designed to monitor the current and eddy activities in the northeastern SCS (Zhang et al., 2013; Guan et al., 2014; Huang et al., 2014). In this study, a long time series of sea surface height (SSH) from satellite altimeter were employed to construct the ensemble to represent the variability of the current and eddy activities. The weekly SSH data from the AVISO (Archiving, Validation and Interpretation of Satellite Oceanographic data) is from Oct.-14-1992 to Dec.-9- 2012. The AVISO product combines data from multiple satellites, including TOPEX/Poseidon, Jason-1, the European Remote Sensing (ERS) Satellites 1, 2 and their successor (Ducet et al., 2000). It has a resolution of 0.25°, the detailed information of the data used in this paper can be found in Table 1. SSH samples can reflect the intraseasonal to interannual variability of the geostrophic current of the study domain. To better understand the variability of the current, we compare the geostrophic currents from the AVISO and the surface currents from the Ocean Surface Current Analysis-Real Time (OSCAR). The OSCAR product is available beginning from October 1999, with a horizontal resolution of 1/3° and 5-day interval (Table 1; Bonjean and Lagerloef, 2002; Johnson et al., 2007), and represents the ocean current (both geostrophic and Ekman components) of the upper 30 m. The temporal resolution of OSCAR and AVISO data are 5-day and 7-day, respectively, to compare the two data, we linearly interpolate the OSCAR and AVISO data into 1-day interval. Figure 3 is the time series of the surface current from OSCAR and the surface geostrophic current from AVISO at station H in Figs. 4 & 5a where the SSH has the largest root-meansquare , , where x_ij is the SSH data at the i^th variable and j^th sample (total n variables in the study domain and m samples)). Because the spatial resolution of OSCAR and AVISO data is different, we need to bilinearly interpolate the two datasets into station H. As shown in Fig. 3, the two datasets are well consistent in terms of both amplitude and phase. The u and v components of the geostrophic current from the AVISO agree very well with those of the OSCAR (the u and v components mean the velocity in the x and y direction, respectively), with correlation coefficients of 0.87 (Fig. 3a) and 0.86 (Fig. 3b), respectively. Furthermore, Fig. 4 depicts the distribution of correlation coefficients between the time series of u and v components from OSCAR and AVISO during January 1999–December 2012. The results show that the small correlation coefficient only appears at a small area around Taiwan, China, and Luzon Island (LI), while the large value appears at most of the study domain (Fig. 4). The mean correlation coefficients at the u and v components in the study area are 0.84 and 0.86, respectively, so the results in Fig. 4 indicate that the u and v components from AVISO are well consistent with that of OSCAR. Therefore, we conclude that the variability of the current in the study domain is primarily from the geostrophic component, and the contribution from the Ekman flow is much smaller. Consequently, selecting SSH data from the AVISO as the initial ensemble to represent the intraseasonal to interannual variability of the current and eddy activities is reasonable.

Due to the lack of large number of threedimensional (3D) observations in the northeastern SCS, 3D weekly temperature results from the HYCOM global 1/12° reanalysis products (Table 1; Experiments 19.0 & 19.1; Fox et al., 2002; Cummings, 2005; Cummings and Smedstad, 2013) are used to evaluate the optimal sampling deployment in the vertical. The HYCOM has a high resolution of 1/12° in the horizontal, and contains 32 hybrid layers in the vertical. The surface forcing of the model is from hourly National Centers for Environmental Prediction (NCEP) Climate Forecast System Reanalysis (CFCR) products, including wind stress, heat fluxes, and precipitation. To improve the forecast capability of both temperature and salinity, a 3D variational algorithm is used to assimilate available satellite altimeter observations, in-situ sea surface temperature (SST) and available in-situ vertical temperature and salinity profiles from XBTs, Argo floats and moored buoys. The HYCOM results exhibit a flow pattern similar to the observations in the LS (Zhang et al., 2010; Zeng et al., 2018), it also reproduces the variation in the Kuroshio intrusion reasonably well on both seasonal and interannual scales (Nan et al., 2011a; Yuan et al., 2014). By comparing the monthly T-S diagram of the HYCOM results with that of the observations, Wang et al. (2015) suggested that the model reproduces the monthly temperature and salinity profiles well in the northeastern SCS. Using the in-situ temperature profiles with different sources in the northeastern SCS, Lian et al. (2019) proved the reliability of the HYCOM results in the northeastern SCS. The weekly temperature results used in this study are from Oct.-2-1992 to Dec.-31-2012. To promote computational efficiency of determining the optimal sampling deployment in the vertical, we only consider the sampling deployment in the surface and middle layers of the ocean, so the 3D temperature data of 19 standard layers (0, 10, 20, 30, 50, 75, 100, 125, 150, 200, 250, 300, 400, 500, 600, 700, 800, 900, 1 000 m) are employed in this study.

To validate the representation skill of the HYCOM reanalysis products, we compare the SST from HYCOM with that of the Advanced Very High Resolution Radiometer (AVHRR). The AVHRR product has the longest record (from late 1981 to the present) of SST measurements from a single sensor design (Table 1). Infrared instruments, like AVHRR, can make observations at relatively high resolution but cannot see through clouds. The daily SST from AVHRR used here is an analysis product that has a spatial resolution of 0.25° (Banzon et al., 2016). Figure 6 depicts the time series of SST from HYCOM and SST from AVHRR at station T where SST from HYCOM has the largest RMS (Fig. 5b). The results show that the SST from HYCOM at station T agrees very well with that of AVHRR (the correlation coefficient is up to 0.98). Figure 7 is the distribution of correlation coefficients between the SST time series from HYCOM and AVHRR from Oct.-2-1992 to Dec.-31-2012, we can see that the SST from HYCOM is also very well consistent with that of AVHRR, the mean correlation coefficient in the study area is up to 0.96. Therefore, we can conclude that the temperature data from HYCOM is suitable to represent the intraseasonal to interannual variation in the northeastern SCS.

Because geostrophic balance is dominant in this region (Figs. 3–4), the surface current is proportional to the SSH gradient which reflects the vertical structures. To further validate the representation skill of the HYCOM reanalysis products, it is necessary to validate the surface current of HYCOM (Table 1). Figure 8 is the time series of the surface current from OSCAR and HYCOM at station H. Because we focus on the weekly temperature products from HYCOM, the high-frequency variability of the surface current from HYCOM within 7 days is filtered. As shown in Fig. 8, although the surface current from HYCOM is larger than that of OSCAR, the two datasets are consistent in terms of phase. The correlation coefficients of the u and v components between the two datasets are 0.69 (Fig. 8a) and 0.57 (Fig. 8b), respectively. The distribution of correlation coefficients between the time series of u and v components from OSCAR and HYCOM also indicates that the small correlation coefficient only appears at a small area, while the large value appears at a large area (Fig. 9), which again proves the reliability of the HYCOM results.

To improve the accuracy of the optimal stations' locations, we remap SSH data to the high-resolution grid of the HYCOM results (SSH data from AVISO is bilinearly interpolated into the HYCOM grid). The SSH data within water depth of 150 m are discarded due to the bad quality of offshore altimetry data. This is also applied to the 3D temperature ensemble to maintain consistency of the study domain, which is shown in Fig. 2. The number of variables (n)/samples (m) for the SSH and 3D temperature datasets are 5 131/1 052 and 92 066/1 050 in this study, respectively.

4 RESULT 4.1 SSH ensemble results

SSH is selected as the ensemble to design an optimal array and assess the present in-situ stations in the northeastern SCS in this study. We assume that SSH has an error variance R of 4, 9, and 16 cm², respectively, to discuss the relation between R and the distribution of optimal stations. The EnOSSE method is used to design the optimal array. When assimilating each optimal station, the initial P^b is updated, which leads to decreasing of the variance of the initial ensemble. A decrease in variance means less uncertainty in the initial sample.

As shown in Fig. 10, the mean SSH (MSSH; the MSSH is calculated from the study period in this paper (from Oct.-14-1992 to Dec.-9-2012)) has large values along the main axis of the Kuroshio in the LS and small values in the southwest of the study domain. The RMS of SSH has two large cores. The largest core is located at the west of LI where the LCG and Luzon cyclonic eddy (LCE) are energetic (Fig. 1; Fang et al., 1998; Gan et al., 2006); another large core is located at the northwest of the LS, where a large number of eddies are generated locally (Yuan et al., 2007; Wang et al., 2008b) or shed from the Kuroshio (Caruso et al., 2006; Yuan et al., 2006; Wang et al., 2008a; Nan et al., 2011a). For the three different error variances R, the first five optimal stations are almost located at the same locations. However, the locations of the subsequent optimal stations vary greatly with different R. The first optimal station is located at the northwest of LI, near the core of the largest RMS, where the energetic LCG occurs (Fig. 10). In winter, the LCG develops very well under the northeasterly monsoon, and dominates the flow pattern in the northeastern SCS, but it is greatly reduced by the southwesterly monsoon in summer (Fig. 1; Fang et al., 1998; Gan et al., 2006). The large seasonal variation of the LCG leads to a large RMS there, which needs more attention when we design an observing system; so, the first optimal station is reasonable. The 2^nd and 4^th optimal stations are at the edge of the LCG. The 3^th optimal station is located at the northwest of the LS, where the Kuroshio path changes strongly, and a large number of mesoscale eddies can be found, leading to the second largest RMS of SSH (Fig. 1; Yuan et al., 2006; Zhuang et al., 2010; Nan et al., 2011a; Zu et al., 2013). The 5^th optimal station is located at the southwest of Taiwan, China, which is energetic due to both Kuroshio loop current and eddy activities, especially the anticyclonic eddies (Nan et al., 2011a; Zu et al., 2013). We notice that the first five optimal stations all appear west of the LS due to the strong variation of currents and eddy activities. In contrast, the Kuroshio is stable east of the LS, only the 12^th optimal station when R=4 and 9 cm² and the 11^th optimal station when R=16 cm² can be found there.

After assimilating all the 18 optimal stations, the basin-averaged RMS of the sample ensemble drops sharply to 3.09, 3.52, and 3.97 cm when R=4, 9, and 16 cm², respectively (Fig. 11a), of which the first optimal station contributes most to the reduction, up to 2.71, 2.61, and 2.49 cm, respectively (Fig. 11b). It means that the northwest of LI is the key region. The first five stations reduce the basin-averaged RMS by 4.66, 4.47, and 4.24 cm when R=4, 9, and 16 cm², respectively, accounting for about 65%, 67%, and 68% of the total reduction (Fig. 11a), which shows that the first five stations reduce the RMS considerably and all subsequent stations result in only very small further reduction. So, the first five optimal stations are most valuable, and should receive more attention when we design an observing system for the domain. Sakov and Oke (2008) pointed out that the locations of the optimal stations were determined not only by the variance of the ensemble, but also by the error variance R; however, they did not further discuss the relation between the error variance R and the RMS of the system. Figure 11a suggests that the RMS decreases faster if the error variance R is smaller. A smaller R represents a more accurate sample ensemble; so, the results show that if a more accurate sample ensemble is used to design an observing array, the same number of optimal stations can yield better observation results.

We now set the observation matrix H to evaluate the present in-situ observation stations shown in Fig. 2. We set H^1×n for different stations and assimilate each station for one cycle, and then assess the performance of each station. First, the results again indicate that each station will reduce the RMS much more; in other words, perform better if the error variance R is smaller (Fig. 12). Second, based on the SSH ensemble, the 18^th station is most close to the core of the largest RMS of SSH (Fig. 10) and performs the best, with the basin-averaged RMS dropping by 2.62, 2.51, and 2.38 cm for R=4, 9, and 16 cm², respectively, after assimilating this station for one cycle. The 4^th station is also near the core of the largest RMS, so it also performs well. The 7^th station, which can only reduce the basin-averaged RMS by less than 0.8 cm (Fig. 12), performs the worst. The results imply that the 18^th and 4^th present stations may have a different spatial correlation scale with the other insitu stations, so a correlation analysis is used next to address this problem.

We define C_i as the mean correlation coefficient (MCC) at the i^th variable to the entire computational domain as follows,

where n is the total variables in the study domain, Corr(i, k) is the correlation coefficient of SSH for the i^th and k^th variables.

Figure 13 is the distribution of the MCC calculated using the SSH ensemble. It shows that the MCC of SSH is large in the northwest of LI and has a small value in the southwest of Taiwan, China, where the 7^th station is located. The 18^th station is located in the area where the MCC of SSH has a relative large value of 0.57 (Fig. 14); it means that the SSH there is well correlated with that of the entire domain. Therefore, the 18^th station performs the best among all the 18 present stations based on the SSH ensemble. The SSH of the 7^th station is relatively uncorrelated with that of the entire domain (Fig. 14), which explains why this station performs the worst. The 12^th to 16^th stations are located at the northwest of the LS, and are all hardly correlated with that of the entire domain, which explains why these stations receive the same small level of assimilation effect (Figs. 13–14). Therefore, we conclude that the present observing system may over sample at the northwest of the LS. Note that the tendency of the MCC for SSH ensemble in Fig. 14 is consistent with that of the reduction of RMS in Fig. 12. This implies that when evaluating an observing system based on certain sample ensemble, if a station is well (badly) correlated with the entire domain, it will reduce the RMS of the initial ensemble much more (less) after assimilating, leading to a good (bad) performance of this station.

To further evaluate the performance of the 18 present stations, we compare the capability of the present and optimal stations to reduce the RMS of the ensemble. In this case, we set R=9 cm², and the results are shown in Fig. 15. The basin-averaged RMS of the initial SSH ensemble decreases gradually after assimilating no matter the optimal stations or present stations, and the optimal stations always perform better than the present stations. Assimilation of the first optimal station and the first present station, the basin-averaged RMS sharply drops to 7.64 and 8.94 cm, respectively. Later, the RMS drops more slowly after assimilating the subsequent stations. The RMS decreases to 3.52 and 5.10 cm after assimilating all the 18 optimal and present stations, respectively. It implies that the future observation stations of the same number, if they are in the key regions as revealed above, can receive better observation effects. We notice that the basin-averaged RMS of SSH ensemble decreases slightly after assimilating the present stations from the 12^th to 16^th stations, it also implies that the present observing system may over sample from the 12^th to 16^th stations.

Now, we investigate the horizontal distribution of the RMS after assimilating all the 18 stations (R=9 cm²). The results are shown in Fig. 16. We can see that small RMS only appears at a small area around the present stations, while big RMS appears within a large area of the southwest of the study domain (Fig. 16a). This suggests that the present stations only perform well in a small surrounding area, while performing badly at the southwest of the study domain, implying that the present observing system may under sample there. On the contrary, for the optimal stations, small RMS occurs in most of the domain (Fig. 16b), indicating that the optimal stations perform well in most of the domain. Figure 17 is the proportion of the area less than different RMSs after assimilating all the 18 optimal and present stations. The area ratio of the RMS < 3 cm is no more than 22% for assimilating of either the 18 optimal or present stations, implying that the optimal stations perform almost the same as the present stations at this level. It shows that the optimal and present stations all perform best within a small region. With the increase of the RMS, the difference of the area ratio at the optimal and present stations becomes larger and larger, and reaches its maximum, up to 47% at RMS=4.5 cm (the optimal and present stations could reduce the RMS of about 94% and 47% of the domain to less than 4.5 cm, respectively). This indicates that the optimal stations perform better for a larger area.

4.2 Three-dimensional temperature ensemble results

Previous research often focused on designing an optimal observing array in the horizontal to monitor the current; few concerns were made about the vertical distribution of observation instruments (such as the CTD machine). To promote observing efficiency and be more economical, the study of determining optimal layers is of the same importance. The 3D weekly temperature at 19 standard levels from the HYCOM global 1/12° reanalysis products are used as the sample ensemble here. We treat all the temperature data in different layers as the samples. Unlike the 2D SSH ensemble, which has only n computational variables (Eq.11; the SSH ensemble has n=5 131 variables), the temperature ensemble used here has close to 19×n computational variables, which requires more computational cost (the water depth is no more than 1 000 m in some regions, leading to less than 19 temperature data in the vertical there).

Based on the EnOSSE method, three different error variances (see below) are set to determine the optimal position in the vertical. Note that in this subsection we mainly focus on the vertical distribution of the first 18 optimal layers, so that the horizontal distribution of the optimal stations is not shown or discussed. As shown in Fig. 18a, when R=0.04, 0.25, and 0.64 (℃)², the first optimal layer is 20 m, 20 m and surface, respectively. Furthermore, for different R, all the 18 optimal layers are located within the upper 200 m, meaning that placing observation instruments and monitoring factors within the upper 200-m depth is more important for the observing system. After assimilating all the 18 optimal layers, the basinaveraged RMS is decreased sharply to 0.48, 0.52, and 0.57℃ when R=0.04, 0.25, and 0.64 (℃)², respectively, of which the first optimal layer contributes most to the reduction, up to 0.37, 0.34, and 0.31℃, respectively (Fig. 18b). This indicates the surface and subsurface layers reduce the RMS considerably and thus are most valuable.

5 SUMMARY

The ensemble-based method for OSSE is employed to design the optimal stations and assess the present stations in the northeastern SCS. To monitor the intraseasonal to interannual variability in the northeastern SCS, the SSH from the satellite altimeter is treated as the initial sample ensemble to discuss the horizontal distribution of the optimal stations.

Based on the SSH ensemble, our results show that the most key region is located at the northwest of LI, where SSH has a very large RMS and the energetic LCG occurs. Other key regions are near the edge of the LCG, at the northwest of the LS and the southwest of Taiwan, China. The first five optimal stations account for more than 65% of the total reduction of the RMS, of which the first optimal station contributes most to the reduction; it shows that the first five stations are most valuable. We further discuss the relation between the error variance R and the RMS of the ensemble. These results suggest that the basinaveraged RMS decreases faster if the error variance R is smaller. A smaller initial R means a more accurate initial sample ensemble, so these results imply that if the optimal stations are designed based on a more accurate and representative initial sample ensemble, the same number of stations can yield better observation effects.

For the present in-situ stations, our results indicate that the 18^th (7^th) station performs the best (worst) (Fig. 12). The 12^th to 16^th stations perform badly, and we conclude that the present observing system may over sample at the northwest of the LS. Only the 4^th and 18^th stations are close to the key region where the energetic LCG occurs (Figs. 10 & 13); so, we conclude that the present observing system may under sample at the northwest of LI. To explain these results, correlation analysis is employed. We show that the 18^th station is located in the area where the MCC of SSH is very large (Fig. 13), which means that the SSH there is well correlated with that of the entire domain. Therefore, the 18^th station performs the best among all 18 present stations. The MCC of SSH has a very small value at the southwest of Taiwan, China, where the 7^th station is located (Fig. 13); it means that the SSH of this station is relatively uncorrelated with that of the entire domain, which leads to the worst performance of the 7^th station.

We further compare the capability of present and optimal stations to reduce the basin-averaged RMS of the ensemble. Our results suggest that the optimal stations always perform better than the present stations; the RMS will decrease to 3.52 and 5.10 cm after assimilating all 18 optimal and present stations, respectively (Fig. 15). It implies that the future observation stations of the same number, if they are in the key regions as revealed above, can have better observation effects. Moreover, the present stations only perform well in a surrounding area, while the optimal stations perform better in a larger area (Fig. 16).

To further improve observing efficiency and be more economical, the study of determining optimal layers in the vertical is conducted. Based on the 3D temperature samples, our results show that the key layer is located within the upper 200-m depth, of which the surface and subsurface layers reduce the RMS considerably and are most valuable for the observing system. This implies that when designing an observing array, we should place more observation instruments in the surface and subsurface layers to promote observing efficiency.

6 DATA AVAILABILITY STATEMENT

The authors duly acknowledge the various data sources for the freely available data: SSH (AVISO) data from https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdTAssh1day.html, OSCAR data from https://coastwatch.pfeg.noaa.gov/erddap/griddap/jplOscar.html, HYCOM reanalysis products from http://tds.hycom.org/thredds/catalogs/GLBu0.08/expt_19.0.html and http://tds.hycom.org/thredds/catalogs/GLBu0.08/expt_19.1.html. SST (AVHRR) data from https://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst2Agg.html.

7 ACKNOWLEDGMENT

We are grateful for the two anonymous reviewers' helpful comments to improve our manuscript. The data analysis is supported by the HPCC, Dr. ZHOU Wei, and Ms. SUI Dandan of the South China Sea Institute of Oceanology, Chinese Academy of Sciences (SCSIO).