Advancing Three-Dimensional Coupled Water Quality Model of Marine Ranches: Model Development, Global Sensitivity Analysis, and Optimization Based on Observation System

Abstract

Marine ranching is a stock enhancement project that has been an important part of aquaculture in China. Due to the lack of scientific management, disasters have occurred, resulting in millions of economic losses. Based on the observation system of marine ranches, a three-dimensional water quality model of marine ranches was developed to study the spatio-temporal variation of key ecological indicators, including the concentrations of chlorophyll-a, zooplankton, detritus, nutrients, and dissolved oxygen (DO). The model is coupled offline with the residual current, temperature, and salinity simulated by a regional oceanic modeling system (ROMS). The conservative characteristic finite difference (C-CFD) scheme is introduced to solve the equations, which guarantees model stability and mass conservation and allows for a larger time step compared to traditional difference schemes. In state-of-the-art water quality models, the biogeochemical processes are parameterized. Due to the complexity of the water quality model, a combination of global sensitivity analysis (GSA) and the adjoint method is introduced as the methodology to optimize the model parameters. Morris’ sampling method is implemented as the GSA method to find out the key factors of the water quality model. The optimization of sensitive parameters with the adjoint method significantly improves the model precision, while the other parameters can be set as empirical values. The results indicate that the combination of GSA and the adjoint method is efficient in the parameter optimization of the water quality model. The model is applicable in marine ranches.

1. Introduction

Marine ranching has been an important part of aquaculture in China [1,2]. Marine ranching is an artificial ecosystem; however, due to the lack of scientific management, it is difficult to detect the potential ecological risks of marine ranches in time. Disasters have occurred, resulting in millions of economic losses. With the development of observation systems in marine ranches, a three-dimensional water quality model with optimization that combines GSA and the adjoint method is essential for the sustainable development of marine ranches, to forecast and analyze the spatio-temporal variations of the environment.

Water quality modeling is an important methodology for studying the ecological dynamics of both rivers and seas [3,4]. Based on the observation system deployed in marine ranches, we were able to develop a water quality model of marine ranches. In state-of-the-art water quality models, the biogeochemical processes are parameterized. Due to the uncertainty of the model parameters, uncertainty exists in these models. Therefore, parameter optimization is essential for improving the simulation. However, the complex solution space of the parameter optimization places high demands on both the amount of observation data and the reliability of the optimization method, especially in complex models [5]. Determining the sensitive parameters and reducing the solution space is essential in parameter optimization. It was pointed out that a change made to one parameter could interact with another highly correlated one [6]. Several studies have expounded that the impact of parameters with weak inter-correlations on models cannot be neglected [7,8]. According to the interactions among the parameters, the sensitivity analysis methods can be divided into two types: local sensitivity analysis (LSA) and GSA [9]. GSA is a kind of sensitivity analysis method that tests the impacts of multiple parameters simultaneously and is able to capture the interactions among the parameters, while LSA analyzes the parameters one by one. In this paper, Morris’ sampling method [10] is implemented as the GSA method due to its high efficiency. In this study, we only optimized the sensitive parameters and kept the insensitive parameters as the empirical values to reduce computation.

Model parameter optimization is the methodology that assimilates observations and optimizes the model parameters to improve the model precision. Among the assimilation methods, 4D-Var has been proven to be one of the most effective and powerful approaches in ecological models [5,11,12,13] and is able to optimize the parameters and improve the simulation precision simultaneously.

This paper is organized as follows: the methods, model settings, and data description are introduced in Section 2; the GSA and parameter optimization results and discussions are presented in Section 3; the conclusions are presented in Section 4.

2. Methods

2.1. Water Quality Model

Based on previous studies of water quality models [4,14,15], an offline-coupled three-dimensional hydrodynamic water quality model was developed. This model is different from the previous models in the following ways: although the water quality model is a nitrogen-based model (which means all the state variables except DO are calculated using the unit mmol N/l), the phosphate and dissolved organic phosphorus (DOP) are taken into consideration. In addition, the nitrification process and the state variable DO are included. The state variables of this model include phytoplankton, zooplankton, detritus, dissolved organic nitrogen (DON), DOP, phosphate (PO₄), ammonium (NH₄), nitrate (NO₃), and DO. The schematic structure of the water quality model is presented in Figure 1.

In the water quality model, the phytoplankton growth rate is described as a function that depends on the water temperature, nutrient limitation, and photosynthetically active irradiance. The temperature effect on phytoplankton growth is described as an exponential function [15,16], and the light control of photosynthesis is described as Steele’s function [17]. Michaelis–Menten functions are used to describe the nutrient limitation [14,18]. The phytoplankton biomass budget includes the photosynthesis process, mortality, grazing by zooplankton, and metabolization. The metabolization process includes respiration and exudation and is parameterized as a function related to the phytoplankton growth rate. More details about the model parameters and governing equations can be found in

Table 1. Values and units of the parameters of the water quality model.

Symbol	Description	Value	Confidence Interval	Unit
ρ	Photosynthetically active irradiance fraction of the total solar irradiance	0.43	0.301–0.559	—
Iopt	Optimal light irradiance	72.5	50.750–94.250	W/m2
kPPT_G	Maximum phytoplankton growth rate	0.8	0.560–1.040	1/day
kPPT_D	Mortality rate of phytoplankton	0.05	0.035–0.065	1/day
kPPT_Z	Zooplankton grazing rate on phytoplankton	0.4	0.280–0.520	1/day
kZPT_D	Mortality rate of zooplankton	0.05	0.035–0.065	1/day
kZPT_N	Excretion rate of zooplankton	0.2	0.140–0.260	1/day
kZPT_F	Grazed rate of zooplankton by fish	0.1	0.070–0.130	1/day
kZPT_R	Respiration rate of ZPT at 0 °C	0.03	0.021–0.039	1/day
kDPT_Z	Zooplankton grazing rate on detritus	0.6	0.420–0.780	1/day
kDPT_B	Detritus remineralization rate	0.05	0.035–0.065	1/day
kDON_NH4	The remineralization rate of DON at 0 °C	0.027	0.019–0.035	1/day
kDOP_B	Dissolved organic phosphorus remineralization rate	0.04	0.028–0.052	1/day
kNH4	Half-saturation concentration for ammonium	0.5	0.35–0.65	mmolN/m3
kNO3	Half-saturation concentration for nitrate	0.5	0.35–0.65	mmolN/m3
kPO4	Half-saturation concentration for phosphorus	0.03	0.021–0.039	mmolP/m3
kNH4_NO3	Oxidation rate of NH4-N at 0 °C	0.053	0.037–0.069	1/day
kDPT_DON	Degradation rate of PN at 0 °C	0.01	0.007–0.013	1/day
KSNO3	Half-saturation constant for nitrate	0.5	0.350–0.650	mmolN/m3
KSNH4	Half-saturation concentration for ammonium	0.5	0.350–0.650	mmolN/m3
KSPO4	Half-saturation constant for phosphorus	0.03	0.021–0.039	mmolP/m3
KSDPT	Half-saturation constant for detritus limitation	0.7	0.490–0.910	mmolN/m3
KSPPT	Half-saturation constant for ingestion	0.6	0.420–0.780	mmolN/m3
Pthre	Threshold for overgrazing on phytoplankton	0.12	0.084–0.156	mmolN/m3
ГPPT_G	Temperature coefficient for phytoplankton growth	0.065	0.046–0.085	1/°C
ГPPT_D	Temperature coefficient for phytoplankton mortality	0.065	0.046–0.085	1/°C
ГZPT_N	Temperature coefficient for zooplankton excretion	0.027	0.019–0.035	1/°C
ГZPT_R	Temperature coefficient for the respiration of ZPT	0.061	0.043–0.0793	1/°C
ГZPT_D	Temperature coefficient for zooplankton mortality	0.05	0.350–0.650	1/°C
ГDON_B	Temperature coefficient for dissolved organic nutrient remineralization	0.065	0.046–0.085	1/°C
ГDPT_B	Temperature coefficient for detritus remineralization	0.05	0.350–0.650	1/°C
ГDON_NH4	The temperature coefficient for the remineralization of DON	0.056	0.039–0.073	°C
ГDPT_DON	Temperature coefficient for the degradation of PN	0.063	0.044–0.082	°C
ГNH4_NO3	Temperature coefficient for the oxidation of NH4-N	0.062	0.043–0.081	°C
δPPT_Z	Assimilation efficiency for phytoplankton	0.8	0.560–1.040	—
δDPT_Z	Assimilation efficiency for detritus	0.7	0.490–0.910	—
DOSNH4	Half-saturation constant of NH4 oxidation consumption dissolved oxygen	0.5	0.350–0.650	mg/L
DOSDON	Half-saturation constant of NH4 transform consumption dissolved oxygen	1.0	0.700–1.300	mg/L
DOSDPT	Half-saturation constant of DPT transform consumption dissolved oxygen	1.0	0.700–1.300	mg/L
rZPT_N	Inorganic nutrient fraction of the excretion of zooplankton	0.75	0.525–0.975
Q0	Solar constant	1368	-	W/m2
R	Albedo of sea surface	0.378	-	—
lk0	Coefficient of light extinction	0.8	-	1/m
lk1	Coefficient of light extinction	0.0088	-	1/m (mg/Chla)
lk2	Coefficient of light extinction	0.054	-	1/m (mg/Chla)2/3
Kext	Light attenuation coefficient	0.1	0.07–0.13	m−1
rChl/PN	Chla/N ratio	1.6	--	mg Chla/(mmol N)
rN_P	N/P ratio in phytoplankton and zooplankton	16	--	mmol N/(mmol P)
μSDO	DO exchange in bottom water	0.8	0.560–1.040	mg/m2/day
ГZPT_R	Temperature coefficient for the respiration of ZPT	0.061	0.043–0.079	°C

and the Appendix A.

Regarding the budget of the DO concentration, the dissolved oxygen is mainly supplied by the photosynthesis process and exchanged at the sea surface. The respiration process of phytoplankton and zooplankton, nitrification, remineralization, and exchange with sediment all consume oxygen.

To improve the calculation efficiency of the water quality model, the conservative characteristic finite difference (C-CFD) scheme was introduced as the difference scheme of this model for numeric calculation [19,20]. The C-CFD scheme guarantees model stability and mass conservation and allows for a larger time step compared to traditional difference schemes. More details about the C-CFD scheme can be found in the Appendix B.

2.2. Global Sensitivity Analysis

Based on the Morris sampling method, the frame of global sensitivity analysis was established. The process can be generalized into three steps—discretization, sampling, and sensitivity analysis. The discretization method is used to normalize and discretize each parameter into the same range. Then the sampling method is implemented to choose the parameter randomly. According to the randomly chosen parameter, the elemental effect (EE) is calculated and the process is continuously repeated until the mean value of EE is convergent. The global sensitivity is defined as the mean value of EE.

2.2.1. Discretization

Firstly, the range of each parameter is normalized into a range of [0, 1]. Then, the parameters are discretized into $\left\{0, \frac{1}{p-1}, \frac{2}{p-1},\dots ,1\right\}$. The simulation results of the water quality model are denoted as a function, $y=y\left(x_1,x_2,\dots x_m\right)$, where $x_1,x_2,\dots x_m$ represents the parameters, and m is the number of parameters. The EE of a certain parameter x_i is defined as shown in Equation (1):

$$E{E}_i=\left|\left[y\left(x_1,x_2,\dots ,x_{i-1},x_i+\delta ,x_{i+1},\dots ,x_m\right)-y\left(x\right)\right]/\delta \right|$$

(1)

$$\delta =\frac{m}{p-1},m=1,2,\dots ,p-1$$

(2)

2.2.2. Sampling

According to Morris’ sampling method, δ is random. Each time a random parameter is chosen, the EE is calculated. The process is repeated until the EE of each parameter is calculated and the mean value of EE is convergent. To achieve this goal, a sampling method is required.

A k-dimensional random vector x* is defined, for which all the elements are randomly chosen from $\left\{0, \frac{1}{p-1}, \frac{2}{p-1},\dots ,1\right\}$.

A (k + 1) × k dimensional down triangle matrix B is defined as follows:

$$B=\left(\begin{array}{c}\begin{array}{lll}0\hfill & 0\hfill & 0\hfill \end{array}\begin{array}{cc} \dots & 0\end{array}\\ \begin{array}{lll} 1\hfill & 0\hfill & 0\hfill \end{array}\begin{array}{cc} \dots & 0\end{array} \\ \begin{array}{lll}1\hfill & 1\hfill & 0\hfill \end{array}\begin{array}{cc} \dots & 0\end{array}\\ \begin{array}{lll}1\hfill & 1\hfill & 1\hfill \end{array}\begin{array}{cc} \dots & 0\end{array}\\ \begin{array}{lll}⋮ \hfill & ⋮\hfill & ⋮\hfill \end{array}\begin{array}{cc} \dots & ⋮\end{array}\\ \begin{array}{lll}1\hfill & 1\hfill & 1\hfill \end{array}\begin{array}{cc} \dots & 1\end{array}\end{array}\right)$$

(3)

A k-dimensional diagonal matrix D* is defined, of which all the elements are randomly set as −1 or 1. A (k + 1)-dimensional vector J_{k + 1,1} is defined, of which all the elements are 1.

A (k + 1) × k dimensional matrix J* is defined as follows:

$$J^{\ast }=\left[\left(2B-J\right)D^{\ast }+J\right]/2$$

(4)

Then, the sampling matrix B* is defined as follows:

$$B^{\ast }=\left(J_{k+1,1}{x}^{\ast }+\delta J^{\ast }\right)P^{\ast }$$

(5)

From the matrices given above, it can be deduced that B* obtains its value randomly. In every two adjacent rows of B*, there is only one different parameter.

2.2.3. Elemental Effect

Considering the two adjacent rows of B*:

$$B^{\ast }\left(j\right)=\left(\begin{array}{cc}{x}_1& \begin{array}{cc}\dots & x_{j1}\end{array} \begin{array}{cc}\dots & x_m\end{array}\\ x_1& \begin{array}{cc}\dots & x_{j2}\end{array} \begin{array}{cc}\dots & x_m\end{array}\end{array}\right)$$

(6)

$$x_{j1}-x_{j2}=\delta$$

(7)

Taking $B^{\ast }\left(j\right)$ as the parameters of the water quality model, the EE of parameter j is shown as follows:

$$E{E}_j=\frac{y\left(x_1,x_2,\dots ,x_{j1}\dots ,x_m\right)-y\left(x_1,x_2,\dots ,x_{j2}\dots ,x_m\right)}{\delta y\left(x_1,x_2,\dots ,x_{j2}\dots ,x_m\right)}$$

(8)

where the function $y\left(x_1,x_2,\dots ,x_j\dots ,x_m\right)$ represents the model output. The process is repeated until the EE of each parameter is calculated and the mean value of EE is convergent. Then, the mean value of the elemental effect is defined as sensitivity and the standard deviation is defined as the correlation between each parameter.

2.3. Parameter Optimization

The parameter optimization method used in this water quality model is the 4D-var method, of which the flowchart is shown in Figure 2. To start the optimization process, the parameters are initialized with prior values. Then the water quality model is run, which is also denoted as the forward model. The simulation results of the state variables are used to calculate the cost function. If the value of the cost function is lower than the criterion, the optimization is terminated; otherwise, the adjoint model is run to adjust the model parameters, and the water quality model is recalculated until the value of the cost function is lower than the criterion.

The cost function is defined as follows:

$$J=\frac{1}{2}{{\displaystyle \sum }}^{ }\left[K{\left(P-\overline{P}\right)}^2\right]$$

(9)

where P represents the simulated concentration of phytoplankton, and $\overline{P}$ represents the observed concentration of phytoplankton. K is the weighting matrix. The elements of the weighting matrix are 1 where the observation of the corresponding grid point is available, and 0 otherwise. The cost function serves as an indicator of the model precision. To optimize the parameters and assimilate the water quality model, the cost function should be minimized. Therefore, the optimization problem is defined as follows:

$$\min :\mathit{J}\left(\mathrm{x}\right)$$

(10)

where $\mathit{J}\left(\mathrm{x}\right)$ is already defined in Equation (9) and x represents the parameters to be optimized. ‘min’ represents the minimum value [5]. The Lagrangian multiplier method is introduced [5] to calculate the gradient. The Lagrangian function is defined as follows:

$$L=J+\underset{\Omega }{{\displaystyle \int }}[P^{\ast }G\left(P\right)]d\Omega$$

(11)

where Ω denotes the model domain; $P^{\ast }$ represents the adjoint variable of P, and G(P) represents the Equation (B1). More details about the equations can be found in Appendix B.

According to the Lagrangian multiplier method, the first-order derivative of Equation (11), with respect to the state variables P and parameters, should be zero. Taking the parameter k_{PPT_G} and I_opt as examples, the equations are shown as follows:

$$\frac{\partial L}{\partial P}=0, \frac{\partial L}{\partial k_{PPT\_G}}=0, \frac{\partial L}{\partial I_{opt}}=0$$

(12)

By solving the equations above, the gradient of the parameters can be determined. The dimensions of the parameters and gradient vector of the cost function are denoted as N and $\nabla J\left(x_i^k\right)$, respectively, with respect to the i-th control variable ($1\le i\le N$) at the k-th iteration step. The following iteration is expressed to solve Equation (10):

$$x_i^{k+1}=x_i^k+{\alpha }_k{d}_k, k=1,2\dots$$

(13)

where d_k represents the gradient of x_k, and α_k represents the adjustment coefficient, which satisfies the strong Wolfe conditions [21,22]. The α_k in Equation (13) is defined as shown in Equations (14) and (15):

$$\left(x_i^k+{\alpha }_k{d}_k\right)\le J\left(x_i^k\right)-c_1{\alpha }_k{\left(d_k\right)}^T\left(d_k\right)$$

(14)

$$\left|\nabla J{\left(x_i^k+{\alpha }_k{d}_k\right)}^T{d}_k\right|\le c_2\left|{\left(d_k\right)}^T{d}_k\right|$$

(15)

where 0 < c₁ < c₂, and T denotes the matrix transposition.

Then, the steepest descent method is implemented to adjust the parameters, and it is expressed as follows:

$$d_k=-\nabla J\left(x_i^k\right)$$

(16)

With Equation (16), Equation (13) can be solved. Then, the optimization problem in Equation (10) can be solved.

2.4. Numerical Experiment Design

2.4.1. Model Settings

The hydrodynamic background is calculated by the three-dimensional regional ocean modeling system (ROMS) [23]. The study site was the Bohai Sea and the Yellow Sea, as shown in Figure 3 (117.5° E–127° E, 34° N–41° N). Horizontal orthogonal grids were applied to calculate the model domain, with a resolution of 3’ × 3’. The time step was set as 20 s. There were 6 uniform layers in the vertical direction. The surface atmospheric forcing fields and shortwave radiation were derived from the Comprehensive Ocean-Atmosphere Data Set (COADS) [24]. The salinity and temperature were obtained from the climatological dataset of the World Ocean Atlas 2013 (WOA13; available online at https://www.nodc.noaa.gov/OC5/woa13, accessed on 1 January 2017). The open boundary was set at 34°N. The sea surface height variations induced by astronomical tides at the open boundary were calculated with the tidal harmonic constants of four dominating tidal constituents, including M₂, S₂, K₁, and O₁, which were obtained from Oregon State University Tidal Inversion Software [25]. The simulation of ROMS started one year before the simulation of the water quality model to make the simulated background forcing stable enough for the water quality model. The ROMS started on 1 January 2015 and terminated on 31 December 2016. The calculated sea surface height, residual current, temperature, and salinity were exported as the background forcing for the water quality model with a time range from 1 January 2016 to 31 December 2016.

The sea surface height calculated by the ROMS in the time range from 1 January 2016 to 31 December 2016 was processed using harmonic analysis to extract the tidal constitutes. The cotidal chart is shown in Figure 4 and is consistent with previous studies on tidal models [26,27,28,29]. Then, the water quality model was coupled offline with the hydrodynamic field.

The horizontal grid and vertical layers of the water quality model are consistent with the hydrodynamic model. With the help of the C-CFD scheme, the time step could be set as 1 h, and the scheme guaranteed model stability and mass conservation. The water quality model was run within the period from 1 January to 31 December 2016.

2.4.2. Data Descriptions

Daily observations of 8 marine ranches, including the chlorophyll-a concentration, water temperature, and salinity at full water depth, with a time range from May to August 2016, were obtained from the observation systems of marine ranches. The locations of the ranches are shown in Figure 3. The Pauta criterion was implemented to remove the spurious values of observations. The chlorophyll-a concentration was observed by SBE WQM, and the temperature and salinity were observed by SEB 63 CTD.

To initialize and calibrate the water quality model, we used the observed temperature, salinity, and open-source datasets to set the initial values of the model. The observed chlorophyll-a concentration was used to calibrate the simulated concentration of phytoplankton and optimize the model parameters.

The open-source datasets included the CSIRO Atlas of Regional Seas (CARS2009; available online at http://www.cmar.csiro.au/cars, accessed on 1 January 2017) and the Sea-Viewing Wide Field-of-View Sensor (SeaWIFS; available online at https://seawifs.gsfc.nasa.gov/, accessed on 1 January 2017). The chlorophyll-a concentration from SeaWIFS was used to set the initial background value of the phytoplankton concentration by the coefficient conversion [15]. The nitrate, phosphorus, and dissolved oxygen from CARS2009 were used to set the initial values of the corresponding state variables. The initial values of the other state variables were set as the empirical values from [15]. Before being applied to the water quality model, both datasets were interpolated into the model grid with the modified Cressman interpolation method described in [30].

3. Results and Discussion

3.1. Simulation Results

The simulated phytoplankton concentration was converted to the chlorophyll-a concentration. As shown in Figure 5, with all the aforementioned observations of eight marine ranches, we compared and calibrated the simulation results and observations from June–August 2016. The water quality model captured an increasing trend of the chlorophyll-a concentration during the summer time, which agrees well with the observations. The mean square errors (MAEs) between the model results and observations are listed in

Table 2. The mean absolute error (MAE) between simulated and observed chlorophyll-a concentration in different ranches.

Station	MAE (mmol/L)
Beihai Ranch	0.25
Liugongdao Ranch	0.24
Sanggouwan Ranch	0.25
Sanggouwanxinan Ranch	0.27
Tianerhu Ranch	0.22
Xiaoshidao Ranch	0.27
Xixiakou Ranch	0.25
Yandunjiao Ranch	0.18

.

To further improve the model, GSA and parameter optimization was implemented as follows:

(1)

GSA was applied to the parameters to determine the sensitive parameters and summarize their effects on different state variables;

(2)

The sensitive parameters were optimized and the effects of optimization on the model results were analyzed.

3.2. GSA Results

The EE was standardized to [0, 1] and ranked as shown in Figure 6. The sensitivity of different parameters and their impact on different state variables were different. The maximum phytoplankton growth rate (k_{PPT_G}) and optimal light irradiance (I_opt) had a significant effect on the phytoplankton, ammonium, nitrate, detritus, phosphorus, DOP, and DO, as shown in Figure 6, because both parameters directly influenced the photosynthesis process, and thus, influenced the growth rate of phytoplankton, which indirectly influenced the nutrient uptake process and DO concentration. The zooplankton biomass was sensitive to k_{ZPT_N} (zooplankton excretion rate), k_{PPT_Z} (zooplankton grazing rate on phytoplankton), k_{DPT_Z} (zooplankton grazing rate on detritus), k_{DPT_Z} (assimilation efficiency for phytoplankton), and δ_{PPT_Z} (assimilation efficiency for detritus) since all of these parameters influenced the grazing and excretion process of zooplankton. The DON concentration was sensitive to k_{DON_NN4} (remineralization rate of DON at 0 °C) and Г_{DON_NH4} (temperature coefficient for the remineralization of DON), which means that the remineralization process was a key factor for the DON concentration. To sum up, among all the parameters, k_{PPT_G} and I_opt were the two most sensitive parameters.

3.3. Parameter Optimization

Since most of the state variables were sensitive to k_{PPT_G} and I_opt, we mainly optimized these two parameters. In the actual simulation, the spatial patterns of the parameters were unknown. Generally, all the parameters were set as the constant values listed in Appendix A. During the optimization, we treated the parameters as constant or spatially varying. Taking the joint parameter k_{PPT_G}/I_opt as an example, three numerical experiments were implemented, including simulation without optimization, simulation with the optimization of constant parameters, and simulation with the optimization of spatially varying parameters. The results of the three numerical experiments, as well as the observed chlorophyll-a concentration, are shown in Figure 6. For further analysis of the change in the model precision, the normalized cost function versus the iterations is shown in Figure 7.

In Figure 8, the Taylor diagram shows that the simulation results with parameter optimization are closer to the observation than that without optimization. The results with the optimization of the spatially varying parameters are better than the constant ones, indicating that the spatial variations of parameters cannot be neglected in water quality modeling, while our method is capable of optimizing the spatially varying parameters. By optimizing the parameters, the correlation between the simulated and observation results improved from −0.1 to 0.8, and the deviation was significantly reduced. In Figure 7b, the cost function with the optimization of the spatially varying joint parameters is smaller than the other two and converges to 0.654.

As described in Section 3.1, the chlorophyll-a concentrations in different marine ranches kept increasing and reached their highest values in August, so in this section, the monthly mean chlorophyll-a concentration in August is shown in Figure 9. According to the spatial distribution of the chlorophyll-a concentration, the highest values of the chlorophyll-a concentration are mainly distributed in the three bays of the Bohai Sea, which matches the spatial distribution of pollutants in the Bohai Sea due to terrestrial discharge [5]. It was pointed out that phytoplankton increases during seasons in which nutrient loadings are high [15]. The high concentration of pollutants provides enough nutrients for the growth of phytoplankton. The chlorophyll-a concentrations in marine ranches are also higher than in adjacent regions.

With respect to the model results, it is possible to evaluate the environmental capacity around the regions of marine ranches and conduct reasonable management to avoid overstocking. The impact of marine ranches on coastal regions is nonnegligible. Compared to previous studies on marine ranching [31,32], we have taken a step towards systematic research on marine ranching.

Based on previous studies of water quality models [4,14,15], we developed, calibrated, and optimized a new coupled three-dimensional water quality model, which was proven to be a powerful tool. Compared to other published work on water quality models [4,5,6,11,12,13,14,15], we applied the C-CFD scheme to a water quality model for the first time, which allowed the model to be more computationally efficient, improving the model precision and reducing the calculation cost at the same time. Furthermore, the new combination of GSA and parameter optimization shows its advantage over the single-parameter optimization method and the models without sensitivity analysis in previous studies [5,11], with higher efficiency and precision.

4. Conclusions

In this study, a three-dimensional water quality model was developed to study the spatio-temporal variation of the concentration of chlorophyll-a, zooplankton, nutrients, detritus, and DO in marine ranches and adjacent regions. The water quality model is coupled offline with the hydrodynamic background field simulated by a ROMS. With the help of the C-CFD scheme, the time step can be set as 1 h, and the scheme guarantees model stability and mass conservation. Due to the limitation of the observations in the marine ranches, we only calibrated the model results of the chlorophyll-a concentrations. The calibration results indicate that the model can make a reasonable estimate of chlorophyll-a concentrations.

Many factors can influence the dynamic process of a complex water quality model. To improve the efficiency of parameter optimization, the GSA method based on Morris’ sampling method was implemented. The sensitivity results show that most of the state variables were sensitive to k_{PPT_G} and I_opt. The parameter optimization results indicate that the optimization of the sensitive parameters can significantly improve the model’s precision, while the other parameters can be set as empirical values.

The water quality model is a powerful tool for marine ranch management. To further utilize the model, we will introduce more observations to calibrate the results of other state variables, and further optimize the model parameters to study the spatio-temporal variations of nutrients and DO.