Marine Environmental Regionalization for the Beibu Gulf Based on a Physical-Biological Model

: A physical–biological ocean model was employed to investigate characteristics of the Beibu Gulf in the northwest South China Sea (SCS) from 2011 to 2015. We adopted the spatially constrained multivariate clustering method to determine the reﬁned marine environmental regionalization using 10 variables from the model output, and compared regionalization differences in ENSO (El Niño–Southern Oscillation) years. The simulated physical and biochemical variables display a wide spectrum of patterns in space and time. The regionalization maps indicated that the Qiongzhou Strait and its adjacent area can be classiﬁed as a separate region, characterized by the rich presence of nutrients, phytoplankton, zooplankton, and detritus, owing to the water invasion from the western Guangdong estuary. As a result of the invasive progression of the SCS, the northern and southern gulf show distinct features over a boundary near 20 ◦ N. In the La Niña year (2011), the classiﬁed boundary of the Qiongzhou Strait-northeastern gulf moved southwards due to enhanced phytoplankton growth. In the El Niño year (2015), the current collision from the northern gulf and SCS resulted in the boundary of the northern and southern gulf moving to approximately 19 ◦ N. These results provide useful guidance on subregional marine management and subregional studies for the gulf.


Introduction
Biogeographic classification and regionalization are critical for developing ecologically representative systems of protected areas. Recent studies have highlighted the need for a more objective division based on ecological parameters, in order to provide an accurate representation of ecological complexity [1][2][3][4][5]. A knowledge of marine regionalization gives insights into climate change research, marine ecosystem protection policies, management of fisheries, and analyses of global food and environmental security [3,4,6]. Current global and regional oceanic regionalization includes the large marine ecosystems (LME) [7], the marine ecoregions of the world (MEOW) [6], the biogeochemical provinces (BGCP) [8,9], biomes and provinces for ecological patterns [10], global ecological marine units (EMUs) [4,5], biogeochemical regions (BGCR) [3], etc. To date, most marine regionalizations depend on representative oceanographic and biogeochemical features, such as sea surface temperature (SST), ocean circulation, bathymetry, and sea surface chlorophyll-a (Chl-a) [7,11,12]. The initial input features for the regionalization cover in-situ data, satellite observations, and model outputs [1,3,8,[11][12][13][14]. Compared to global oceanic regionalization, regional marine regionalization can integrate physical and biochemical factors with a finer resolution, thus providing detailed information for marine ecological policies and management. In particular, regional marine regionalization is commonly implemented on the scale of gulfs [3]. with a finer resolution, thus providing detailed information for marine ecological policies and management. In particular, regional marine regionalization is commonly implemented on the scale of gulfs [3].
The Beibu Gulf (16°00′-21°30′ N, 105°40′-111°00′ E, also named the Gulf of Tonkin) is a semi-enclosed gulf in the northwestern South China Sea (SCS) (Figure 1). The depth of the Beibu Gulf is less than 100 m, with a mean depth of 50 m [15][16][17]. The gulf is rich in natural resources, especially fisheries for both China and Vietnam. However, the gulf environment is at increasing risk, since the coastal rim of the Beibu Gulf has been undergoing rapid economic growth in recent years [18,19]. Thus, it is essential to conduct studies to inform management for the environmental protection of the gulf. Previous studies indicated that the distribution of environmental biochemical variables (e.g., Chl-a, nutrients, zooplankton, fish species, benthos species) within the gulf is strongly influenced by monsoons and river inputs [20][21][22][23][24]. Therefore, the environmental conditions of the gulf are likely diverse in time and space. In addition, physical conditions, such as temperature, salinity, and current, also have non-negligible impacts on the temporal and spatial distributions of environmental variables [16,[25][26][27][28]. Clustering both physical and biochemical variables using an unsupervised machine learning algorithm for the Beibu Gulf could help to define the impacted geographic range of specific dynamic and environmental conditions. In particular, a series of papers [21,29,30] have explored the water mass types in the Beibu Gulf. The temperature and salinity measurements collected by the Chinese Offshore Investigation and Assessment (COIA) cruises in 2006-2007 were used to classify the water masses in the eastern Beibu Gulf. They found that the water masses in summer can be subdivided into five categories: upper water, northern cold water, subsurface water, intermediate water, and lower water; the mixed water in winter can be subdivided into low salinity water and high salinity water. The chemical and biological parameters, including dissolved oxygen, alkalinity, nitrite, active silicate, and Chl-a, were further used to define five water masses: diluted water, mixed water masses, surface water, subsurface water masses, and bottom water. More recently, Wang et al. [31] found that the water of the In particular, a series of papers [21,29,30] have explored the water mass types in the Beibu Gulf. The temperature and salinity measurements collected by the Chinese Offshore Investigation and Assessment (COIA) cruises in 2006-2007 were used to classify the water masses in the eastern Beibu Gulf. They found that the water masses in summer can be subdivided into five categories: upper water, northern cold water, subsurface water, intermediate water, and lower water; the mixed water in winter can be subdivided into low salinity water and high salinity water. The chemical and biological parameters, including dissolved oxygen, alkalinity, nitrite, active silicate, and Chl-a, were further used to define five water masses: diluted water, mixed water masses, surface water, subsurface water masses, and bottom water. More recently, Wang et al. [31] found that the water of the Beibu Gulf can be divided into nine groups using a fuzzy c-means clustering method based on year-long sea surface temperature (SST) and Chl-a data derived from the 4 km MODIS Aqua. However, the clustering approaches mentioned above are limited to specific seasons, areas, and parameters. In addition, the La Niña and El Niño effects may result in variations in ecological systems; comparing the regionalization under different climate conditions may extract the explicit boundaries of the ecological impacts in the gulf.
To better understand the subregional characteristics under different climate conditions and to promote environmental management in the Beibu Gulf, here we used 10 variables from the output of a three-dimensional physical-biological numerical ocean model developed by the Marine Environment Committee (MEC) of Japan [32], as parameters to regionalize the study region. First, we used the mean variables from 2011 to 2015 at monthly intervals to carry out the regionalization. Next, we compared the regionalization using the annual mean outputs of 2011, 2013, and 2015, which represent the La Niña, normal, and El Niño years, respectively, according to the observation data of NOAA [33].

Numerical Model
In the present study, we mainly analyzed the output of the physical-biological model in the area of the Beibu Gulf. In the MEC ocean model, the vertical eddy viscosity, as well as diffusivity, and the horizontal diffusion coefficient, were parameterized using the Mellor-Yamada 2.5-level turbulence closure scheme [34], and the Smagorinsky turbulence closure scheme [35], respectively. The model uses the hydrostatic balance with the Boussinesq approximation for convective flows to solve the primitive equations. The study domain (  Table 1. The detailed equations of the hydrodynamic model were provided by the literatures [36][37][38]. The physical model outputs temperature (TEM), salinity (SAL), and current field (CF). We used a biological model developed by Nakata [42]. It is coupled off-line to the MEC ocean model and uses the TEM, SAL, and CF as forcing. In the biological model, physical, biological, and chemical processes govern the ocean ecological processes via: here B represents the biochemical variables; ∂B ∂t is the partial derivative of the biochemical variables relative to time; u, v, and w are the flow velocities in the x, y, and z directions; A c and K c represent the diffusion coefficients in the horizontal and vertical directions; S(B) is the biogeochemical source/sink term; and the equation for each variable is given in Document S1. As a brief summary, S(B) represents the effect of the biochemical processes on the biochemical variables relative to time, including the processes that control the dynamics of phytoplankton (PHY), cell quota nitrogen (Q N ) of PHY, cell quota phosphorus (Q P ) of PHY, zooplankton (ZOO), dissolved inorganic nitrogen (DIN), dissolved inorganic phosphorus (DIP), non-living particulate organic matter (POM, also known as detritus), dissolved organic matter (DOM), and dissolved oxygen (DO). The biological model schematic is shown in Figure 2.
here represents the biochemical variables; is the partial derivative of the biochemical variables relative to time; , , and are the flow velocities in the , , and directions; and represent the diffusion coefficients in the horizontal and vertical directions; ( ) is the biogeochemical source/sink term; and the equation for each variable is given in Document S1. As a brief summary, ( ) represents the effect of the biochemical processes on the biochemical variables relative to time, including the processes that control the dynamics of phytoplankton (PHY), cell quota nitrogen (QN) of PHY, cell quota phosphorus (QP) of PHY, zooplankton (ZOO), dissolved inorganic nitrogen (DIN), dissolved inorganic phosphorus (DIP), non-living particulate organic matter (POM, also known as detritus), dissolved organic matter (DOM), and dissolved oxygen (DO). The biological model schematic is shown in Figure 2. The biological model shared the same grid as the physical model, with a time step of 60 s. The concrete model settings of parameters are in Document S1. We exported the model simulated values at monthly intervals for further analysis.
The modeled TEM, SAL, PHY, DIN, and DIP in the year 2007 were evaluated against COIA data in varied depths and seasons in our previous study [43]. Note that the COIA observations represent the most comprehensive in-situ investigation available in the Chinese territorial waters over the recent decade.
Here, we further used MODIS Aqua data for sea surface temperature (SST) and Chla (represents PHY) for the validation in the years of 2011-2015. The Taylor diagram [44] of the model validation against MODIS Aqua annual mean data of 2011-2015 is shown in Figure 3. The correlation coefficients are above 0.7 for Chl-a and above 0.8 for SST, suggesting that the model and remote sensing data have a high correlation for all study years. The root mean square error (RMSE) shows the model results and the remote sensing data in close distances. The normalized standard deviations (STDs) of the model are within 1.0, meaning that the model results had a smaller STD than that of the observations. Hence, the overall variation between the model results and remote sensing data for the Chl-a and SST are within acceptable error ranges [45]. Therefore, the comparison against remote sensing data also assesses the accuracy of the model results of 2011-2015. The biological model shared the same grid as the physical model, with a time step of 60 s. The concrete model settings of parameters are in Document S1. We exported the model simulated values at monthly intervals for further analysis.
The modeled TEM, SAL, PHY, DIN, and DIP in the year 2007 were evaluated against COIA data in varied depths and seasons in our previous study [43]. Note that the COIA observations represent the most comprehensive in-situ investigation available in the Chinese territorial waters over the recent decade.
Here, we further used MODIS Aqua data for sea surface temperature (SST) and Chl-a (represents PHY) for the validation in the years of 2011-2015. The Taylor diagram [44] of the model validation against MODIS Aqua annual mean data of 2011-2015 is shown in Figure 3. The correlation coefficients are above 0.7 for Chl-a and above 0.8 for SST, suggesting that the model and remote sensing data have a high correlation for all study years. The root mean square error (RMSE) shows the model results and the remote sensing data in close distances. The normalized standard deviations (STDs) of the model are within 1.0, meaning that the model results had a smaller STD than that of the observations. Hence, the overall variation between the model results and remote sensing data for the Chl-a and SST are within acceptable error ranges [45]. Therefore, the comparison against remote sensing data also assesses the accuracy of the model results of 2011-2015.

Spatially Constrained Multivariate Clustering
First, in order to eliminate the scale difference of the original data value, all ten variable were normalized between 0 and 1 by the min-max method (Equation (1)).
We performed regionalization and mapping processes in ArcGIS Pro 2.5 (Esri, Redlands, CA, USA). Spatially constrained multivariate clustering (SCMC) was used for unsupervised clustering analysis. This method is based on geographic location clustering to ensure that each clustering result has a geographical proximity of the geographic location, and clusters of non-adjacent polygons will not be marked as the same cluster [46,47]. To avoid tiny region areas that would have no ecological meaning, we set the parameter of minimum points of a cluster to be 50. We chose the trimmed Delaunay triangulation (TDT) option as the spatial constraint, because it is appropriate for point features and will ensure that a feature will only be included in a cluster if at least one other cluster member is a natural neighbor (Equations (2)-(4)) [48].
where is the total square sum and is the regression square sum. reflects the degree of retention of changes in the original data after grouping. Therefore, the larger

Spatially Constrained Multivariate Clustering
First, in order to eliminate the scale difference of the original data value, all ten variable were normalized between 0 and 1 by the min-max method (Equation (1)).
We performed regionalization and mapping processes in ArcGIS Pro 2.5 (Esri, Redlands, CA, USA). Spatially constrained multivariate clustering (SCMC) was used for unsupervised clustering analysis. This method is based on geographic location clustering to ensure that each clustering result has a geographical proximity of the geographic location, and clusters of non-adjacent polygons will not be marked as the same cluster [46,47]. To avoid tiny region areas that would have no ecological meaning, we set the parameter of minimum points of a cluster to be 50. We chose the trimmed Delaunay triangulation (TDT) option as the spatial constraint, because it is appropriate for point features and will ensure that a feature will only be included in a cluster if at least one other cluster member is a natural neighbor (Equations (2)-(4)) [48].
where TSS is the total square sum and RSS is the regression square sum. R 2 reflects the degree of retention of changes in the original data after grouping. Therefore, the larger R 2 value of a particular variable, the better the variable can be distinguished. n c is the number of clusters, n i is the number of features in the cluster i, n v is the number of variables, V k ij is the kth variable value of jth features in ith cluster, V k is the average value of the kth variable, and V k i is the average value of the kth variable in the cluster i.
To determine the grouping number, we tested various options ranging from 2 to 15 and calculated the Calinski-Harabasz pseudo F-statistic, which can be used to determine the largest ratio of between-group difference to within-group similarity at the cluster number of f (Equation (5)) [4,49]. Then a normalized average slope function (ASF) was used to determine the optimum number of clusters: after the point of f , there would be a rapid decrease; if there are three or more consecutive slower decreases or inflexion (within 1% of each slope) after the point of k, the "elbow" point k would be the threshold of acceptable flatness (TAF), where an additional clustering number does not significantly reduce the pseudo F-statistic. Hence, k is applied as the final optimum number of cluster [3,12]. The in-region probabilities were also calculated by the evidence accumulation of 100 parallel clustering runs.
In Equation (5), n is the number of features, and n c is the number of clusters.

Regionalization of 2011-2015 Means
The results of the regionalization of 2011-2015 means are shown in Figure 5a. As seen, the inner Beibu Gulf can be clustered into three to five regions with some similarities and differences among depths, and the in-region membership (Figure 5b) is high (>90%) in most regions. The box plots of the normalized six representative variables for each depth, including the maximum, first quartile, medium, third quartile, and minimum values (excluding outliers) in each region are presented in Figure 5c (all variables are in Figure S4), and the average and deviation of the original values of the representative variables in each region is summarized in Table 2 (all variables are in Table S1). In the following results and discussion, we focus on the regionalizations in the inner gulf, while the regionalizations in the outer gulf will be the references.

Regionalization of 2011-2015 Means
The results of the regionalization of 2011-2015 means are shown in Figure 5a. As seen, the inner Beibu Gulf can be clustered into three to five regions with some similarities and differences among depths, and the in-region membership (Figure 5b) is high (>90%) in most regions. The box plots of the normalized six representative variables for each depth, including the maximum, first quartile, medium, third quartile, and minimum values (excluding outliers) in each region are presented in Figure 5c (all variables are in Figure S4), and the average and deviation of the original values of the representative variables in each region is summarized in Table 2 (all variables are in Table S1). In the following This region has the characteristics of higher average PHY, ZOO, DIN, DIP, DOM, and POM than the other regions in the inner gulf. In the central gulf, there is a boundary separating the northern and the southern part of the inner gulf near 20° N. The northern region has higher PHY, ZOO, DIN, DIP, DOM, and POM, with a lower SAL than the southern region in general. At the 55.0-m depth layer, there is a region in the southern mouth of the gulf (Region 8), which has in-between values of PHY, ZOO, DIN, and DIP, compared with its adjacent regions (Region 1 and Region 5).

Comparison for the Regionalization in ENSO and Normal Climate Years
Besides regionalization for the 2011-2015 means, we compared the difference of the marine environmental regionalization under three different climate years, in which 2011 was a La Niña year, 2013 was a normal year, and 2015 was an El Niño year, according to NOAA [33]. The optimum regionalization of the years 2011, 2013, and 2015 remained generally similar to the results of the 2011-2015 means ( Figure 6). However, there were some differences between them. The most significant difference is that there is an individual region (Region 4) on the estuary of the Guangxi coast for the years 2011 and 2015 at the 2.5-m depth layer, and the in-region probabilities in this area is high (nearly 100%) (Figure 7).

Discussion
Since the marine environmental region of a specific marine area can delineate the environmental characteristics of the area, it provides an objective framework for upcoming defined management, policy, and research in specific zones [2,3,6]. In marine areas lacking frequent field investigations such as the Beibu Gulf, using multiple physical and biochemical variables in a high-resolution distribution derived from a physical-biological model to generate multivariate regionalization can continuously monitor and forecast the ecosystem and dynamics, giving insights into dynamic management for different areas [14]. In our study, we used a model simulated on 2011-2015 means and annual means of the representative climate years, and applied to regionalize the Beibu Gulf, in order to find the regional characteristics under different climates in the gulf.

Variable Distribution
As discussed in our previous paper [43], the ecosystem of the Beibu Gulf is strongly influenced by interface processes (monsoon, water flux, and river input). In winter and autumn, the study area is controlled by the northeast monsoon, while in summer it is controlled by the southwest monsoon [20][21][22][23][24]43,50]. In the 2011-2015 means distribution, the effect of northeast monsoon exceeded the southwest monsoon, so the current direction in the Beibu Gulf was westward or cyclonic (Figure 4c). However, the southwestward intrusion of the SCS water from the southern gulf was still significant, which can be seen in the SAL distribution (Figure 4b). Although the TEMs vary across the gulf in different seasons at the sea surface, which was lower in the northern gulf (<24 °C) and higher (>25 °C) in the southern gulf in winter, while higher in both northern and southern gulf (>30 °C) in summer [43,51], the 2011-2015 means did not show significant variations among the gulf at the sea surface. The distribution of SAL reflects the influence of the SCS through the Qiongzhou Strait in the northern and southern gulf mouth. The physical variables of TEM, SAL, and CF form a background for the distribution characteristics of nutrients

Discussion
Since the marine environmental region of a specific marine area can delineate the environmental characteristics of the area, it provides an objective framework for upcoming defined management, policy, and research in specific zones [2,3,6]. In marine areas lacking frequent field investigations such as the Beibu Gulf, using multiple physical and biochemical variables in a high-resolution distribution derived from a physical-biological model to generate multivariate regionalization can continuously monitor and forecast the ecosystem and dynamics, giving insights into dynamic management for different areas [14]. In our study, we used a model simulated on 2011-2015 means and annual means of the representative climate years, and applied to regionalize the Beibu Gulf, in order to find the regional characteristics under different climates in the gulf.

Variable Distribution
As discussed in our previous paper [43], the ecosystem of the Beibu Gulf is strongly influenced by interface processes (monsoon, water flux, and river input). In winter and autumn, the study area is controlled by the northeast monsoon, while in summer it is controlled by the southwest monsoon [20][21][22][23][24]43,50]. In the 2011-2015 means distribution, the effect of northeast monsoon exceeded the southwest monsoon, so the current direction in the Beibu Gulf was westward or cyclonic (Figure 4c). However, the southwestward intrusion of the SCS water from the southern gulf was still significant, which can be seen in the SAL distribution (Figure 4b). Although the TEMs vary across the gulf in different seasons at the sea surface, which was lower in the northern gulf (<24 • C) and higher (>25 • C) in the southern gulf in winter, while higher in both northern and southern gulf (>30 • C) in summer [43,51], the 2011-2015 means did not show significant variations among the gulf at the sea surface. The distribution of SAL reflects the influence of the SCS through the Qiongzhou Strait in the northern and southern gulf mouth. The physical variables of TEM, SAL, and CF form a background for the distribution characteristics of nutrients (DIN and DIP), followed by PHY and ZOO, which decrease from north to south in general, reflecting the seasonal patterns of the water exchange driven by the monsoon. In winter and autumn, water with large amount of nutrients from the western Guangdong estuaries flows through the Qiongzhou Strait into the gulf, making the northeastern gulf a highvalue area for nutrients, and thus promoting the growth of PHY and ZOO, and increases of POM and DOM [52]. Studies have shown that, in summer, lager amount of nutrients are carried by the increase of the river input flow into the gulf, making the estuaries have high concentrations in PHY, ZOO, DOM, and POM, especially on the Guangxi coast [43,[53][54][55]. The model simulated 2011-2015 mean distributions shown in Figure 4 keep the seasonal features mentioned above. The distributions show gradual changes in each environmental variable, which may visually indicate a coarse ecological boundary in the Beibu Gulf.

Marine Environmental Regionalization
Marine environmental regionalization can refine ecological boundaries under the comprehensive interactions of all variables. From the model output distributions, the northern and southern gulf showed varying patterns for most variables.

The Northern Gulf
The area in and to the west of the Qiongzhou Strait can be identified as an individual region in this study. In the study of Wang [31], a cluster could also be identified in the area near the Qiongzhou Strait and to the west of the Leizhou Peninsula. Due to the transport effect of the Qiongzhou Strait, nutrients have relatively high values in this area [23]. This region reflects the influence of a westward current driven by the northeast monsoon [31,51,56]. However, in some other studies, the area to the west of the Leizhou Peninsula was classified as the same region as the area in the northern gulf, as the water mixed [21,29,30,57]. Moreover, in a study of fish community spatial patterns, the area to the west of the Leizhou Peninsula had a special fish community [58], which also verifies the distinct characteristic here.
The estuary on the Guangxi coast could be identified as an individual region at the 2.5-m depth layer in the year of 2011 and 2015. In the regionalization of 2011-2015 means and 2013 it was not identified as an individual region. However, the in-region probabilities of 2011-2015 means was lower here (~50%) (Figure 5b), meaning that this area may also have an alternative chance of being an individual region. Most variables here, especially DOM, PHY, and SAL, have significant differences from the adjacent area ( Figure 4). In the model results and field data, these variables have obvious high-value areas on the surface of the estuary area from spring, and reach a maximum in the summer, reflecting the obvious influence of the river inputs [22,23,31]. According to the study of Chen [21], the water at the surface layer of the estuaries on the Guangxi coast was classified as diluted water. However, this area was smaller for the 2011-2015 means and 2013, so it was not identified as an individual region. In other studies by Chen [29,30], when using TEM and SAL as clustering parameters, this cluster did not exist, indicating that it is the considerable amount of nutrients carried by the rivers and the nutrient-related variables that make this region independent.

The Central Gulf
The differences between the northern and southern parts of the Beibu Gulf are reflected at the three layers, so there is usually a zonal boundary near 20 • N (while in 2015 it was near 19 • N). The boundary current of the SCS usually intrudes into the gulf near 19-20 • N, and then flows southwards along the Vietnam coast, carrying the high-salinity SCS water to the southern gulf (Figure 4b,c). In previous studies, the boundary of the mixed water and the surface water [21] and the boundary of different water mass patterns [29,30] also existed here. The northern gulf is affected by diluted water from the Guangxi estuaries and diluted water from the western Guangdong estuaries through the northern SCS and the Qiongzhou Strait, while the southern gulf is mainly affected by the SCS current from the southern gulf mouth. The regionalization results can determine the specific ecological boundaries of the northern and southern gulf at different depths.

The Southern Gulf
The regionalization result in the southern gulf is relatively complex. The regions in the southern gulf located gradually from the SCS to the inner Beibu Gulf. The more a region penetrates into the gulf, the lower SAL value it has; so it reflects the invasive stages of the SCS. From the perspective of theCF, the hydrodynamic characteristics in the southern gulf are more complicated, and vary by season [16,23,26,59]. The regionalization results reveal that the invasion characteristics of the SCS impact the Beibu Gulf spatially. In other studies of water mass clustering, except in summer [30], this area did not show such complex patterns [21,29]. However, from the results of our model, the horizontal distribution of each environmental variable is closely related to the current and salinity patterns (Figure 4b,c and Figure 9). Therefore, the regionalization of this study identified invasive stages of the SCS current.

The Climate Effect
According to a remote sensing study in the SCS (including the Beibu Gulf), the average SST dropped in the La Niña year of 2011, but did not have a significant increase in the El Niño year of 2015 [60]. In the model simulated results, we also found the TEM distribution in the Beibu Gulf had this trend in general ( Figure 10).  The physiological activities of PHY are significantly affected by TEM, hydrodynamic, and environmental conditions [32,42]. Therefore, the distribution of PHY in 2011 reflected that the range and the value of high PHY areas in the northern gulf was greater than in 2011 and 2015, indicating that the lower sea temperature in the La Niña year was more conducive to PHY growth in the Beibu Gulf [61]. In the La Niña phase, the wind strengthened, resulting in an increase in CV and seawater mixing, so the concentration of PHY increased [62,63]. From the results of the model simulated outputs (Figures 9-11), we found that the La Niña year of 2011 had lower TEM in general, especially in the northern gulf, and the area of high PHY was larger than that of 2013 and 2015 in the northeastern gulf, but that it moved to south. As a result, the region of the Qiongzhou Strait-northeastern gulf in 2011 moved further southwards than that of 2013 and 2015 (Region 5). In addition, in the mouth of the southern gulf, the CV in the southwest direction was greater than in 2013 and 2015, thus the high salinity SCS water intruded more into the inner gulf ( Figure 9). The region in the southern gulf (Region 6 in the upper and middle layers, Region 8 in the bottom layer) shifted further into the inner gulf.

Conclusions
In this study, we used 10 physical and biochemical variables of the physical-biological model output as parameters, and applied an unsupervised machine learning clustering method to regionalize the upper, middle, and bottom layers of the Beibu Gulf.
The model simulated distribution of the variables indicates that, driven by interface processes (such as monsoon, water flux, and river inputs), the northern and the southern gulf show different patterns, in which the northern gulf is mainly influenced by the highnutrient water from the western Guangdong estuaries through the Qiongzhou Strait and river inputs from the Guangxi coast, while the southern gulf is mainly influenced by the invasion of the SCS water.
The regionalization results show that the Beibu Gulf has unique characteristics in the Qiongzhou Strait and the northeastern gulf. This area can be identified as an individual region, reflecting the impact of the water from the western Guangdong estuaries as it advects through the Qiongzhou Strait to the Beibu Gulf. The estuaries on the Guangxi coast showed specific characteristics at the upper layers in 2011 and 2015. Therefore, this can be identified as an individual region, which represents the influence of the river inputs. In the central Beibu Gulf, there is a zonal boundary near 20° N, indicating the difference between the northern and the southern parts of the Beibu Gulf, which can help to identify the ecological boundary of the Beibu Gulf at different depths, because it can be regarded as the boundary between the influence of water in the northern gulf and SCS invasion in In 2015, the SST did not have a significant increase in the Beibu Gulf compared to the normal year of 2013. However, the regionalization boundary of the northern and the southern gulf in the upper and middle layers moved to south (19 • N) compared to 2011 and 2013 (20 • N). In our model output, the northward current from the southern gulf and southward current from the northern gulf confronted near 19 • N in 2015. From the distribution of the SAL (Figure 9), the invasion of high salinity water from the SCS extended to the gulf, to near 19 • N, and had a significant difference with the lower salinity water from the northern gulf. The distribution of the SAL was strongly affected by the current [43], and had most significant difference here, thus it resulted in the boundary of the northern and the southern gulf moving to here. The CF reflects that the total impact of the northeast monsoon on the central gulf was weaker than the southwest monsoon in 2015 ( Figure 9).

Conclusions
In this study, we used 10 physical and biochemical variables of the physical-biological model output as parameters, and applied an unsupervised machine learning clustering method to regionalize the upper, middle, and bottom layers of the Beibu Gulf.
The model simulated distribution of the variables indicates that, driven by interface processes (such as monsoon, water flux, and river inputs), the northern and the southern gulf show different patterns, in which the northern gulf is mainly influenced by the high-nutrient water from the western Guangdong estuaries through the Qiongzhou Strait and river inputs from the Guangxi coast, while the southern gulf is mainly influenced by the invasion of the SCS water.
The regionalization results show that the Beibu Gulf has unique characteristics in the Qiongzhou Strait and the northeastern gulf. This area can be identified as an individual region, reflecting the impact of the water from the western Guangdong estuaries as it advects through the Qiongzhou Strait to the Beibu Gulf. The estuaries on the Guangxi coast showed specific characteristics at the upper layers in 2011 and 2015. Therefore, this can be identified as an individual region, which represents the influence of the river inputs. In the central Beibu Gulf, there is a zonal boundary near 20 • N, indicating the difference between the northern and the southern parts of the Beibu Gulf, which can help to identify the ecological boundary of the Beibu Gulf at different depths, because it can be regarded as the boundary between the influence of water in the northern gulf and SCS invasion in the southern gulf. The regions in the southern gulf reside gradually from the SCS to the inner Beibu Gulf, reflecting the invasive stages of the high-salinity SCS water at different depths.
In the La Niña year (2011), the region of the Qiongzhou Strait-northeastern gulf in 2011 moved southwards, reflecting that lower TEM can enhance the growth of PHY near the Qiongzhou Strait, thus this region became distinct from the adjacent area. In the El Niño year (2015), although there was no significant TEM increase in the gulf, the boundary of the northern and southern gulf moved to near 19 • N, due to the enhanced invasion of SCS.
Marine environmental regionalization can provide insights for refined environmental management for different purposes, and can provide insights into biogeographic subregionalization and the ENSO effects in further studies of the Beibu Gulf.
Supplementary Materials: The following supplemental materials that are available online at https: //www.mdpi.com/2077-1312/9/2/187/s1, Document S1: The equations and parameter settings for the biological model, Figure   Data Availability Statement: The data are not publicly available because some of them are classified as confidential.

Conflicts of Interest:
The authors declare no conflict of interest.