**1. Introduction**

With the fast development of the coastal economy and marine aquaculture, estuaries are suffering a lot from terrestrial and oceanic pollutants. Understanding the characteristics of pollutants and estimating the marine environmental capacity (MEC) in estuaries provides a theoretical foundation for coastal management. The deteriorated environment of bays feeds back to and constrains the economy of coastal cities. Therefore, identifying sources and sinks of pollutants, and estimating the estuarine environmental capacity is quite important for coastal development (Yoon et al., 2020; Halpern et al., 2008; Syvitski et al., 2005) [1–3].

The main reason for the deterioration of the offshore water environment is that a large number of pollutants produced on land are discharged into it, including inorganic nitrogen, active phosphate, and heavy metals, etc. (Chen et al., 2008; Cui et al., 2013) [4,5].

The first step to conduct total discharge control is to assess the MEC of pollutants. (Wu et al., 2005) [6]. MEC refers to the maximum load of pollutants that can be accommodated in a specific sea area under the premise of making full use of the marine selfpurification capacity without causing pollution damage (Linker et al., 2013) [7].

Sanmen Bay is located on the coast of the East China Sea (Figure 1a). It is a macro-tidal turbid estuary, with a maximum tidal range of approximately 2 m at the bay mouth and suspended sediment concentration (SSC) of 1.192 kg/m<sup>3</sup> at the middle of the bay. It has

**Citation:** Yao, Y.; Zhu, J.; Li, L.; Wang, J.; Yuan, J. Marine Environmental Capacity in Sanmen Bay, China. *Water* **2022**, *14*, 2083. https:// doi.org/10.3390/w14132083

Academic Editor: Karl-Erich Lindenschmidt

Received: 24 May 2022 Accepted: 23 June 2022 Published: 29 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

four tributaries. With the fast development of human activities, the pollution problem is becoming increasingly heavy in the bay. has four tributaries. With the fast development of human activities, the pollution problem is becoming increasingly heavy in the bay.

**Figure 1.** (**a**) Model domain for Sanmen Bay, (**b**) Location of observation stations, (**c**) Initial conditions of water exchange model (The concentration distribution in Sanmen Bay is set as 1 unit, red in the figure, the inflow outside the bay and at the boundary is set as 0, dark blue in the figure), (**d**) Water quality control stations in Sanmen Bay. **Figure 1.** (**a**) Model domain for Sanmen Bay, (**b**) Location of observation stations, (**c**) Initial conditions of water exchange model (The concentration distribution in Sanmen Bay is set as 1 unit, red in the figure, the inflow outside the bay and at the boundary is set as 0, dark blue in the figure), (**d**) Water quality control stations in Sanmen Bay.

Identification of the contamination source is a part of the systematical analysis to find the pollution source fundamentally. Because how much a certain water body can hold is on the premise of pollution intensity which also plays the role to decide how to distribute for such a specific object. According to the various forms of pollution being discharged into the water body, Identification of the contamination source is a part of the systematical analysis to find the pollution source fundamentally. Because how much a certain water body can hold is on the premise of pollution intensity which also plays the role to decide how to distribute for such a specific object.

pollution sources can be divided into point source pollution and non-point source pollution. Non-point source pollution is referred to as water body pollution caused by rainfall run-off, this kind of contaminant entered into the soil or underground water body in a wide, microcrystalline, and dispersive way. Since the 1960s, research on non-point source pollution has caught the interest of scientists throughout the world. So far non-point source pollution feature, impact factor, load quantification on pollution output, and mechanism of pollutant migration and transformation have made great achievements. According to the various forms of pollution being discharged into the water body, pollution sources can be divided into point source pollution and non-point source pollution. Non-point source pollution is referred to as water body pollution caused by rainfall run-off, this kind of contaminant entered into the soil or underground water body in a wide, microcrystalline, and dispersive way. Since the 1960s, research on non-point source pollution has caught the interest of scientists throughout the world. So far non-point source pollution feature, impact factor, load quantification on pollution output, and mechanism of pollutant migration and transformation have made great achievements.

A certain amount of contaminant emission poured is permitted. Because natural water body holds a certain amount of environmental capacity for a sort of pollution. Total emission radically depends on assimilative capacity, distribution of total water pollutants based on the limit of water environmental capacity. A certain amount of contaminant emission poured is permitted. Because natural water body holds a certain amount of environmental capacity for a sort of pollution. Total emission radically depends on assimilative capacity, distribution of total water pollutants based on the limit of water environmental capacity.

According to the environmental quality survey results of the Sanmen Bay sea area from 2015 to 2016, and combined with the marine environmental function zoning, the researchers evaluated and classified the current situation of the environmental quality of the sea area (Liang et al., 2021) [8]. A 2-D hydrodynamic and pollutant model of Sanmen Bay is established based on Deift3D. Combined with the current hydrological status of Sanmen Bay, the transport and diffusion laws of COD, TP, and TN in the bay are analyzed (He et al., 2018) [9].

After conducting Clean Water Act in 1972, although industrial and municipal pollution has been under effective control by carrying out NPDES (National Pollutant Discharge Elimination System), the water quality has not been radically improved. Previous research showed that it is non-point source pollution that mainly polluted the rivers, lakes, and surface waters in estuaries. Also, non-point source pollution polluted underground water and degradation of wetland ecosystem. For such a reason, the Clean Water Act includes a guide rule, named TMDL (Total Maximum Daily Load), aiming at controlling both point source pollution and non-point source pollution. The acting emphasis is on figuring out non-point source pollution load and elimination in key water areas. Such a concept (MEC) was raised by Japanese scholars in the environment field in 1968. Japanese researcher Yanowokio (1968) claimed that environmental capacity is determined by environment quality standards, i.e., keeping total contamination loads within a permitted limit [10].

Streeter and Phelps (1952) raised a simple S-P model, which is the earliest form of water quality model [11]. The development of the water quality model can be concluded into 5 periods: 1925–1960, BOD-DO coupled model was raised based on the S-P model; during 1960–1965, spatial variety, physical, kinetic factors, and temperature were introduced into the model as a state variable, in the same the heat exchange between air and the water surface was also considered; during 1965–1970, as the computer started to be applied, people increasingly deeply understand biochemical oxygen consumption; the calculation method was developed from 1-D to 2-D; during 1970–1975, water quality model has developed into mutually non-linearity. In the last 20 years, space dimensionality has been into 3-D, the emphasis of study gradually turned into improving the dependability and evaluation of the model. In this way, dependability on the calculation of environmental capacity has been enhanced [12–16].

The eutrophication degree of seawater, the enrichment degree of heavy metals in sediments and the potential ecological hazard effects were comprehensively analyzed.

In this study, we take Sanmen Bay as an example, to study the characteristics of water exchangeability and marine environmental capacity, using both numerical models and field data.

Sanmen Bay is located in the monsoon subtropical humid climate zone. Affected by the monsoon climate, it has four distinct seasons and a mild climate. The weather changes are complex and disastrous weather is frequent. Disastrous weather of different degrees can be encountered in all seasons. The seasonal variation of temperature is obvious. Sanmen Bay has abundant rainfall, mainly from March to September. The whole year can be roughly divided into two rainy seasons and two relatively dry seasons. The wind direction of Sanmen Bay varies with seasons. Typhoons, rainstorms, and sudden small-scale disastrous weather occur from time to time.

#### **2. Materials and Methods**

#### *2.1. Model Descriptions*

A calibrated two-dimensional hydrodynamic model was used to reproduce the coastal oceanic and environmental characteristics of Sanmen Bay. Delft 3D simulates water surface elevation, velocity, water quality, waves, and morphology. Flow, Hydrodynamics were used in this study. The governing equations of the Delft3D hydrodynamic model in the vertical sigma coordinate system and horizontal curvilinear coordinate system are expressed as follows: namics were used in this study. The governing equations of the Delft3D hydrodynamic model in the vertical sigma coordinate system and horizontal curvilinear coordinate system are expressed as follows:

A calibrated two-dimensional hydrodynamic model was used to reproduce the coastal oceanic and environmental characteristics of Sanmen Bay. Delft 3D simulates water surface elevation, velocity, water quality, waves, and morphology. Flow, Hydrody-

*Water* **2022**, *14*, x FOR PEER REVIEW 4 of 19

$$\frac{\partial \zeta}{\partial t} + \frac{1}{\sqrt{\mathsf{G}\_{\xi\overline{\xi}}\sqrt{\mathsf{G}\_{\eta\eta}}}} \frac{\partial \left( (d+\zeta)\mathcal{U}\sqrt{\mathsf{G}\_{\eta\eta}} \right)}{\partial \xi} + \frac{1}{\sqrt{\mathsf{G}\_{\xi\overline{\xi}}\sqrt{\mathsf{G}\_{\eta\eta}}}} \frac{\partial \left( (d+\zeta)\mathcal{V}\sqrt{\mathsf{G}\_{\xi\overline{\xi}}} \right)}{\partial \eta} = (d+\zeta)\mathcal{Q} \tag{1}$$

where *ξ* is the coordinate direction under the Delft3D curvilinear coordinate system corresponds to the X-axis of the rectangular coordinate system, *η* is the Y-axis, *ζ* is the height of the water surface above the zero scale line of the Z coordinate, d is the depth from the zero scale line of the Z coordinate to the bottom of the water, *U* is the velocity for X-axis, *V* is the velocity for Y-axis, p *Gηη* is conversion coefficient for X-axis and p *Gξξ* is conversion coefficient for Y-axis. where is the coordinate direction under the Delft3D curvilinear coordinate system corresponds to the X-axis of the rectangular coordinate system, is the Y-axis, is the height of the water surface above the zero scale line of the Z coordinate, d is the depth from the zero scale line of the Z coordinate to the bottom of the water, *U* is the velocity for X-axis, *V* is the velocity for Y-axis, √ is conversion coefficient for X-axis and √ is conversion coefficient for Y-axis.

The three-dimensional convection-diffusion equation in the water quality module is as follows: The three-dimensional convection-diffusion equation in the water quality module is as follows:

$$\frac{\text{Cov}\sigma\text{\partial C}}{\partial t} + v\_x \frac{\partial \text{C}}{\partial x} - D\_x \frac{\partial^2 \text{C}}{\partial x^2} + v\_y \frac{\partial \text{C}}{\partial y} - D\_y \frac{\partial^2 \text{C}}{\partial y^2} + v\_z \frac{\partial \text{C}}{\partial z} - D\_z \frac{\partial^2 \text{C}}{\partial z^2} = \text{S} + f\_R(\text{C}, t) \tag{2}$$

where *C* is substance concentration, *D* is diffusion coefficient, *S* is the inflow term, *fR*(*C*, *t*) is the reaction term. where *C* is substance concentration, *D* is diffusion coefficient, *S* is the inflow term, (,) is the reaction term.

#### *2.2. Model Configurations 2.2. Model Configurations*

**2. Materials and Methods** *2.1. Model Descriptions*

The model domain contains Sanmen Bay and its adjacent seas. The open ocean boundary is from 28◦310 N to 29◦260 N, the east part can extend to 122◦270 E. Jiaojiang River runoff is considered at the open boundary. The measurements were conducted in the bay between the Nantian station and the Linhai station (Figure 2a). The model domain contains Sanmen Bay and its adjacent seas. The open ocean boundary is from 28°31′ N to 29°26′ N, the east part can extend to 122°27′ E. Jiaojiang River runoff is considered at the open boundary. The measurements were conducted in the bay between the Nantian station and the Linhai station (Figure 2a).

**Figure 2.** Tidal elevation in (**a**) Jiantiao station and (**b**) Gangdi station (2 December 2009 to 12 December 2009).

According to the calculation region, to generate orthogonal curvilinear grid scatter with the 2-D hydrodynamic model. The Mesh quantity is 629 × 674. In the study region, the minimal edge length is about 50 m, the maximum edge length is about 220 m, and the computational time step is set at 60 s. Charted depth data is gained from the historical chart. ADI method is used to solve the problem. To calculate the environment capacity, we firstly determine water quality goal by dividing functional the water area, then conduct a numerical simulation to consider the quantitative response relation between pollution emissions and water environmental quality.

*Water* **2022**, *14*, x FOR PEER REVIEW 5 of 19

chart. ADI method is used to solve the problem.

**Figure 2.** Tidal elevation in (**a**) Jiantiao station and (**b**) Gangdi station (2 December 2009 to 12 De-

with the 2-D hydrodynamic model. The Mesh quantity is 629 × 674. In the study region, the minimal edge length is about 50 m, the maximum edge length is about 220 m, and the computational time step is set at 60 s. Charted depth data is gained from the historical

According to the calculation region, to generate orthogonal curvilinear grid scatter

The computational domain is about 100 km on the x-axis and 110 km on the y-axis.

The model was run for a further one month for the period from 1 December 2009 to

The grid consisted of 423,946 elements, forming a mesh of orthogonal curvilinear grid with variable cell widths ranging from 220 m in the area of the open sea to 50 m in Sanmen Bay. The bathymetry data were interpolated linearly and were corrected with the satellite

The computational domain is about 100 km on the x-axis and 110 km on the y-axis. The grid consisted of 423,946 elements, forming a mesh of orthogonal curvilinear grid with variable cell widths ranging from 220 m in the area of the open sea to 50 m in Sanmen Bay. The bathymetry data were interpolated linearly and were corrected with the satellite chart (Figure 1a). *2.3. Model Validation* Hourly water elevation from 3 December to 12 December in 2009 is taken from

The model was run for a further one month for the period from 1 December 2009 to 30 December 2009. Jiantiao station and Gangdi station are compared with model elevations in Figure 2. The result of verification basically agreed with the measured data. The relative error of the

To calculate the environment capacity, we firstly determine water quality goal by dividing functional the water area, then conduct a numerical simulation to consider the quantitative response relation between pollution emissions and water environmental quality. whole process could be controlled within 10%. According to spring tides and neap tides verification, the result can basically reflect the tide wave transformation of Sanmen Bay. The current velocity and direction of each verification point coincides with the meas-
