A Three-Dimensional Coupled Hydrodynamic-Ecological Modeling to Assess the Planktonic Biomass in a Subalpine Lake

Abstract

In this study, a coupled three-dimensional hydrodynamic-ecological model was developed to comprehensively understand the interaction between the hydrodynamics and ecological status of a lake. The coupled model was utilized to explore the hydrodynamics, water quality, and ecological status in an ecologically rich subalpine lake (i.e., Tsuei-Feng Lake (TFL), located in north-central Taiwan). The measured data of water depth, water temperature, water quality, and planktonic biomass were gathered to validate the coupled model. The simulated results with a three-dimensional hydrodynamic and water quality-ecological model reasonably reproduced the variations in observed water depth, water temperature, water quality, and phytoplankton and zooplankton biomass. Sensitivity analysis was implemented to determine the most influential parameter affecting the planktonic biomass. The results of sensitivity analysis indicated that the predation rate on phytoplankton (PRP) significantly affects the phytoplankton biomass, while the basal metabolism rate of zooplankton (BMZ) importantly affects the zooplankton biomass. Furthermore, inflow discharge was the most important environmental factor dominating the phytoplankton and zooplankton biomass of TFL. This implies that the runoff in the catchment area caused by rainfall and the heavy rainfall induced by climate change may affect the planktonic biomass of the lake.

1. Introduction

Freshwater lakes provide important and valuable resources such as water supply, environmental and ecological services, recreation, and landscape appreciation. Due to rapid economic development, population expansion, and land cultivation, watershed areas have been overly developed, resulting in excessive nutrient input into lakes during rainfall, causing serious damage to water quality and eutrophication and even resulting in an imbalance of ecosystems. Therefore, reducing external nutrient loads into lakes to prevent deterioration of water quality and eutrophication for sustainable lake management has become an extremely important task [1,2,3,4,5,6].

To fully understand the water quality and ecological health of lakes, on-site observation is an indispensable measure. However, no matter how dense the number of measurement stations is, the water quality and ecological conditions in spatial and temporal distributions are still insufficient, and the cost of on-site observation is too expensive. Another alternative is to utilize mathematical modeling, which can be employed to describe spatiotemporal variations in hydrodynamics, water quality, and ecological variables and can be applied to predict the water quality and ecological conditions due to human activities and after adjustment strategies [7,8,9]. Moreover, coupled hydrodynamic-ecological models have become useful and effective tools for exploring the interactions among physical, biogeochemical, and ecological processes and for further water quality and ecological management in surface waters [10,11,12,13,14]. These models have been widely developed and applied to different water bodies, such as lakes, reservoirs, wetlands, rivers, estuaries, bays, and oceans [15,16,17,18,19,20,21,22].

Recently, different kinds of one-dimensional, two-dimensional, and three-dimensional hydrodynamic-ecological coupled models have been developed and applied to investigate water quality and ecosystems in lakes [8,16,23,24,25,26,27,28,29,30,31]. The one-dimensional and two-dimensional hydrodynamic-ecological coupled models are not so complicated for running the models to obtain the simulation results with minimal computational cost for evaluating water quality and ecological status in lakes, but these kinds of models are not utilized to resolve the horizontal and vertical dimensions simultaneously. Therefore, three-dimensional hydrodynamic-ecological coupled models display useful tools to obtain such information for environmental assessment and management, even though they require more computational resources [31].

Due to the complex bathymetry and topography of lakes as well as the complexity of the water environment in lakes, the development of a coupled three-dimensional hydrodynamic-ecological model based on the transport and transformation processes of flow, pollutants, and ecological variables is urgently necessary. The objective of this study was to develop a three-dimensional hydrodynamic-ecological coupled model (SCHISM-Ecol) to simulate the hydrodynamics, water quality, and biomass of phytoplankton, zooplankton, and fish in Tsuei-Feng Lake (TFL), which is located in the north-central mountains of Taiwan. The coupled model was validated with observations of water depth, water temperature, water quality, and planktonic biomass. The model sensitivity analysis of parameters was conducted with the validated model to determine the important parameters affecting the biomasses of phytoplankton and zooplankton. The validated model was then applied to explore when the lake’s environment changed, which was the most crucial environmental factor affecting the planktonic biomass.

2. Materials and Methods

2.1. Description of the Study Site and Data

Tsuei-Feng Lake (TFL) is one of the largest alpine lakes in Taiwan. It is located in the mountainous area of north-central Taiwan between Taiping Mountain and Dayuan Mountain (Figure 1). TFL is surrounded by a secondary cypress forest that was replanted five decades ago after selective logging [32]. It has an altitude of 1850 m and a circumference of 600 m. It is shaped like a gourd, measuring 160 m wide at its widest point and 15 m at its narrowest point. The annual precipitation in the lake exceeds 3600 mm, and the mean annual temperature is approximately 13.4 °C. The period of full water level is from September to November, especially affected by typhoons, with an area of up to 25 ha and a depth of 7 m. During the dry season from January to April, a large area of grassland is exposed, making the lake appear as a large lake and a small lake, with a total area of 8 ha and a depth of 3 m [33]. As it is located in the Taiping Mountain National Forest Recreation Area (TMNFRA), there is a trail around TFL to provide tourists with recreation and sightseeing services. The principle of ecological engineering was employed throughout in the construction of the trail, and care was taken to avoid destroying the original surface and disturbing the underlying biological features. TFL has a more obvious inflow to the southeast of the lake and an outflow to the northwest of the lake. The deepest part of its bottom is 4.56 m, and a monitoring station buoy was set up (see Figure 1). According to Chiu et al. [32], the degree of eutrophication is between mesotrophic and eutrophic. The concentration of total phosphorus ranged from 10 to 70 µg/L, and the mean dissolved organic carbon and pH were 2.9 mg/L and 5.9, respectively.

The water depth, water temperature, and meteorological data were all obtained from the buoy station every ten minutes using autonomous sensors. The dissolved oxygen (DO) was measured in situ with a DO sensor. The lake water samples were sampled and taken back to the laboratory for measuring nutrients (ammonium nitrogen, dissolved organic nitrogen, phosphorus, and dissolved organic phosphorus), dissolved organic carbon (DOC), and ecological state variables (phytoplankton biomass and zooplankton biomass). These measured data were collected in 2017 and 2018 and were utilized for hydrodynamic-ecological model validation.

In order to collect the phytoplankton and zooplankton data, a water sample of 1 L was collected at 0–50 cm water depth from the lake surface for the determination of phytoplankton biomass. Zooplankton sample was collected in the study lake by filtering 5 L water sample obtained from 0–50 cm water depth from the lake surface with zooplankton sampling net (mesh size = 330 µm). Three replicates of phytoplankton and zooplankton samples were collected at TFL during each sampling. Both phytoplankton and zooplankton samples were fixed with formaldehyde, stored in labeled bottoms and transported back to the laboratory at 4 ℃. In the laboratory, individual water samples were passed through pre-combusted Whatmann GF/F 0.7μm filters to collect particulate organic matter (POM), which is composed of mainly (≥93%) phytoplankton in a mixture of phytoplankton and detritus. The relative importance of phytoplankton in the POM was confirmed by microscopic examination at the 500–1000x magnification of 100 mL subsample from each water sample. Each zooplankton sample was examined at 80–100x magnification using a Leica research microscope, and identified into three taxa (i.e., Copepoda, Cladocera, and Rotifera) to make a composite sample. The relative abundance of each taxon in each zooplankton sample was calculated and recorded. Prior to further analysis, the individual phytoplankton samples and zooplankton samples were dried at 60 ℃ to constant weight. After drying, samples were homogenized in an agate mortar, recorded in dry weight (to the nearest 0.1 mg), and stored in a desiccator. The biomass of the phytoplankton and zooplankton samples were determined by dividing the dry weight with the corresponding volume of the field-collected water samples (mg/L). Triplicates of phytoplankton and zooplankton samples (~2.0 mg) were put into tin capsules and weighed with an accuracy of ±1 μg. The samples were analyzed by the elemental analyzer (NA 1500, Fison, Italy) in the Plant Ecophysiological Lab of National Taiwan University. The results were expressed as % amount of carbon in each sample, and these values can be converted to carbon content (C) of the samples (mg C/L) by back-calculated with their weight.

2.2. Hydrodynamic Model

A three-dimensional semi-implicit cross-scale hydroscience integrated system model (SCHISM) was employed to simulate the hydrodynamics in the lake [34,35]. The model is an open-source supported modeling system derived from the early version of the SELFE model [36]. SCHISM adopts an efficient and accurate semi-implicit and finite-volume method with an Eulerian-Lagrangian algorithm to solve the turbulence-averaged Navier-Stokes equations under hydrostatic and Boussinesq approximations. The governing equations of SCHISM can be expressed as follows.

The continuity equation is

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0$$

(1)

The momentum equation is

$$\frac{Du}{Dt}=fv-g\frac{\partial \eta }{\partial x}-\frac{1}{{\rho }_0}\frac{\partial p_a}{\partial x}-\frac{g}{{\rho }_0}{\displaystyle {\int }_z^{\eta }\frac{\partial \rho }{\partial x}}dz+\frac{\partial }{\partial x}(K_{mx}\frac{\partial u}{\partial x})+\frac{\partial }{\partial y}(K_{my}\frac{\partial u}{\partial y})+\frac{\partial }{\partial z}(K_{mz}\frac{\partial u}{\partial z})$$

(2)

$$\frac{Dv}{Dt}=-fu-g\frac{\partial \eta }{\partial y}-\frac{1}{{\rho }_0}\frac{\partial p_a}{\partial y}-\frac{g}{{\rho }_0}{\displaystyle {\int }_z^{\eta }\frac{\partial \rho }{\partial y}}dz+\frac{\partial }{\partial x}(K_{mx}\frac{\partial v}{\partial x})+\frac{\partial }{\partial y}(K_{my}\frac{\partial v}{\partial y})+\frac{\partial }{\partial z}(K_{mz}\frac{\partial v}{\partial z})$$

(3)

The mass balance transport equation is

$$\frac{\partial C}{\partial t}+\frac{\partial (uC)}{\partial x}+\frac{\partial (vC)}{\partial y}+\frac{\partial (wC)}{\partial z}=\frac{\partial }{\partial x}(K_{cx}\frac{\partial C}{\partial x})+\frac{\partial }{\partial y}(K_{cy}\frac{\partial C}{\partial y})+\frac{\partial }{\partial z}(K_{cz}\frac{\partial C}{\partial z})$$

(4)

The equation of state is

$$\rho =\rho (T,p)$$

(5)

where x and y denote the horizontal Cartesian coordinate; z expresses the vertical coordinate with positive upward; $u,$$v,$ and $w$ represent the velocity components in the x, y, and z directions, respectively; $t$ is time; $f$ denotes Coriolis factor; $g$ expresses the acceleration of gravity; $\eta$ represents the free-surface elevation; ${\rho }_0$ is the reference water density; $\rho$ denotes the water density which considers the effects of water temperature ($T$) and hydrostatic pressure ($p$); $p_a$ is the atmospheric pressure at the free surface; $K_{mx},K_{my},$ and $K_{mz}$ represent the eddy viscosity coefficients in x, y, and z directions, respectively; $C$ denotes the water quality-ecological state variables; and $K_{cx},K_{cy},$ and $K_{cz}$ express the turbulent diffusivities in the x, y, and z directions, respectively.

The turbulence closure of SCHISM employs the generic length-scale model of Umlauf and Burchard [37] with the stability function of Kantha and Clayson [38]. Since the cross-scale, semi-implicit scheme is adopted in SCHISM, the model allows large time steps supported by lower- and higher-order numerical methods [39]. Wetting and drying processes and a quadratic bottom drag equation for bottom stress are also considered in the model. A detailed description regarding the numerical solution methods, horizontal and vertical grids, and boundary conditions can be found in Zhang et al. [34,35].

SCHISM has been widely utilized in various research fields [18,39,40,41,42,43,44,45,46,47,48], but limited cases have been utilized to simulate the hydrodynamics in lakes. Therefore, SCHISM was coupled with an ecological model to develop a hydrodynamic-ecological model (i.e., SCHISM-Ecol) for exploring the planktonic biomass in a subalpine lake.

2.3. Phytoplankton and Zooplankton Modules

The kinetic equation in the mass balance transport equation can be written as

$$\frac{\partial C}{\partial t}=KC+R$$

(6)

where $K$ denotes the kinetic rate (time⁻¹) and $R$ expresses the source/sink term (mass volume⁻¹ time⁻¹).

Phytoplankton can be divided into two state variables: periphyton and drifting phytoplankton, which include green algae, diatoms, and others. The subscript, $i$, is employed to represent these two groups. The kinetic equations for phytoplankton biomass can be expressed as

$$\frac{\partial P{B}_i}{\partial t}=(G{P}_i-BM{P}_i-PR{P}_i)\cdot P{B}_i+\frac{\partial }{\partial z}(W{S}_i\cdot P{B}_i)+\frac{WP{B}_i}{V}$$

(7)

where PB denotes the phytoplankton biomass (mass volume⁻¹); GP represents the production rate of phytoplankton (time⁻¹); BMP expresses the basal metabolism rate of phytoplankton (time⁻¹); PRP represents the predation rate on phytoplankton (time⁻¹). It denotes the rate that the phytoplankton is preyed on by zooplankton; WS is the settling velocity (Length time⁻¹); WPB denotes the external loads of phytoplankton (mass time⁻¹); and V expresses the cell volume (volume).

Since copepods frequently appear during winter and summer seasons in TFL, the zooplankton is divided into two state variables: copepods and other species. The subscript, $j$, is employed to represent these two groups. The kinetic equation for zooplankton biomass can be written as

$$\frac{\partial Z{B}_j}{\partial t}=(G{Z}_j-BM{Z}_j-PR{Z}_j)\cdot Z{B}_j+\frac{WZ{B}_j}{V}$$

(8)

where ZB denotes the zooplankton biomass (mass volume⁻¹); GZ expresses the production rate of zooplankton (time⁻¹); BMZ represents the basal metabolism of zooplankton (time⁻¹); PRZ indicates the death rate of zooplankton (time⁻¹); and WZB is the external load of zooplankton (mass time⁻¹).

Figure 2 illustrates the schematic diagram of water quality-ecological state variables. In the figure, the state variables DOC, DON, DOP, NH₄, NO₃, PO₄, DO, and COD represent dissolved organic carbon, dissolved organic nitrogen, dissolved organic phosphorus, ammonium nitrogen, nitrate nitrogen, orthophosphate, dissolved oxygen, and chemical oxygen demand, respectively. Fecal coliform and fish biomass are also included in the model simulation.

2.4. Model Setups

To build unstructured grids (triangles) for numerical simulation, bathymetric TFL data were collected from Academia Sinica in Taiwan. Surface-water modeling system (SMS) software was utilized to build unstructured grids and to greatly reduce the time for manual grid generation. In the horizontal plane, a total of 11,589 unstructured grids were generated for model simulation. The grid size is approximately 5 m. A pure S coordinate was utilized with 5 S-levels in the vertical plane (see Figure 3). A time step ($\Delta t$ = 120 s) was used in the hydrodynamic-ecological model simulations. A constant bottom roughness height, z_o = 0.001 m, was adopted in the modeling. A 60-day spin-up time was utilized to achieve an equilibrium condition.

2.5. Statistical Error for Model Performance

Three statistical indexes consisting of the mean absolute error (MAE), the Pearson correlation coefficient (CC), and the predictive skill (skill) were employed to quantify the model performance of SCHISM-Ecol. Model validation to compare the model results with observational data can achieve the minimization of MAE and the maximization of r (within positive range) and skill value. Equations (9) and (11) for MAE, CC, and skill [49] can be expressed as

$$\frac{\partial Z{B}_j}{\partial t}=(G{Z}_j-BM{Z}_j-PR{Z}_j)\cdot Z{B}_j+\frac{WZ{B}_j}{V}$$

(9)

$$CC=\frac{{\displaystyle {\sum }_{i=1}^N(C_m^i-\overline{C_m})(C_o^i-\overline{C_o})}}{\sqrt{{\displaystyle {\sum }_{i=1}^N{(C_m^i-\overline{C_m})}^2{\displaystyle {\sum }_{i=1}^N{(C_o^i-\overline{C_o})}^2}}}}$$

(10)

$$skill=\frac{{\displaystyle {\sum }_{i=1}^N{\left|C_m^i-C_o^i\right|}^2}}{{\displaystyle {\sum }_{i=1}^N{(\left|C_m^i-\overline{C_o}\right|+\left|C_o^i-\overline{C_o}\right|)}^2}}$$

(11)

where N represents the total number of modeled/observed data; $C_m^i$ denotes the modeled value at data point i; $C_o^i$ expresses the observed value at data point i; $\overline{C_m}=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{C}_m^i}$; and $\overline{C_o}=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{C}_o^i}$.

The MAE denotes the average deviation between the simulation and measurement. The Pearson correlation coefficient (CC) was employed to detect the linear relationship between the modeled results and measured data. The r value is in the range between +1 and −1 to denote a positive and negative relationship, respectively. If r = 0, no relationship exists. A predictive skill equal to 1 denotes perfect agreement of the model-measurement relationship, while a predictive skill equal to 0 expresses complete disagreement of the model measurement [45,50].

2.6. Metric for Sensitivity Analysis

Sensitivity analysis is an important procedure to understand which parameters affect the simulation results. To quantify the results of parameter sensitivity analysis, the maximum rate is adopted in this study. The equation to express the maximum rate $(MR)$ is shown as

$$MR=\frac{C_{baseline}-C_{sen}}{C_{baseline}}\times 100\%$$

(12)

where $C_b$ denotes the phytoplankton and zooplankton biomass for the baseline simulation and $C_s$ represents the phytoplankton and zooplankton biomass for the sensitivity run.

3. Model Validation

The measured data, including water depth, water temperature, and water quality-ecological state variables, were gathered from September 2017 to July 2018 to quantify the model accuracy for validating SCHISM-Ecol. Since the spin-up time of the model was set to 60 days, the model was conducted starting on 1 July 2017. The initial conditions for water depth and water temperature were specified as 1.5 m and 20 °C, respectively. The mean values of measured water quality and ecological variables were set as the initial conditions for simulating the SCHISM-Ecol model.

3.1. Water Depth

To precisely predict the water depth of the TFL, in this study, the daily discharges at the inflow and outflow locations in the lake served as an input to drive the model simulation for predicting the time-series water depth.

Figure 4a depicts the time-series inflow and outflow discharges in the TFL. Due to the multiple heavy rainfalls and Typhoon Khanun, several peak inflow discharges occurred during the periods of October 11 and 15. Based on the recorded data, the total rainfall during the periods was 1023.7 mm. Figure 4b illustrates the comparison between simulated and observed time-series water depths from 1 September 2017, to 31 July 2018, indicating that the high inflow discharges quickly responded to the high water depth in the lake. However, the simulated results reproduced the observed water depths well. It also indicated that the peak water depth reached 10.13 m at 7:00 am on 16 October 2017 (Figure 4b). It also showed that the model reproduced the observed normal and low water depths well. The values for MAE, CC, and skill were 0.03 m, 0.92, and 0.95, respectively. Based on the model-data comparison, excellent model performance was achieved.

3.2. Water Temperature

The precise prediction of water temperature in lakes is a crucial requirement for modeling water quality-ecological state variables and further water quality management. To validate the 3D coupled model, the measured water temperature data at the buoy station gathered from 1 September 2017, to 31 July 2018, were utilized for comparison with the modeling results. Meteorological data were employed to calculate heat and momentum fluxes, providing the surface boundary condition for the 3D coupled model.

Figure 5 compares the simulated water temperatures in the time series and measured results. The water temperature at 0.5 m below the water surface displayed day-night variations. The water temperature reached 22.2 °C during the summer but decreased to 4.1 °C during the winter. The water temperature during the daytime was higher than the nighttime temperature. This figure also indicated that the simulated water temperatures accurately reproduced the measured results. The values of MAE, CC, and skill were 0.29 °C, 0.96, and 0.98, respectively, between the simulated and observed water temperatures. The comparison of the observed and simulated water temperatures in vertical profiles at the buoy station is illustrated in the Supplementary Materials (Table S1). This result indicated that the simulated results satisfactorily reproduced the observed vertical water temperature on different dates, except for the measured dates on 5 September 2017. The simulated water temperatures exhibited more vertical mixing at the bottom layers than the observational data on 5 September 2017. This would be the reason that the stronger turbulent mixing in the vertical direction was calculated by the model.

Table 1. Statistical errors based on differences between simulated and measured water temperatures.

Year/Month/Date	MAE (°C)	CC	Skill
2017/09/05	0.58	0.96	0.60
2017/10/10	0.40	0.78	0.66
2017/11/21	0.18	0.91	0.56
2018/01/09	0.10	0.80	0.79
2018/03/26	0.24	0.94	0.91
2018/05/15	0.86	0.68	0.66
2018/06/12	0.90	0.78	0.86
2018/07/17	0.15	0.98	0.99

presents the statistical errors of the comparison between the simulated and observed water temperatures in the vertical profile on different measured dates. The computed water temperatures in the vertical profile satisfactorily matched the measured results during the summer and winter seasons but slightly underestimated the measured results in the bottom layer in September 2017. Based on the skill values, the model performance in the vertical profiles ranged between very good and excellent scores.

3.3. Water Quality and Ecology

To ensure that the water quality and planktonic modules in the SCHISM-Ecol model are reliable and predictive, the measured water quality and ecological variable data were employed to validate the coupled model.

Several reports have documented that the validation of water quality-ecological models is more difficult than that of hydrodynamic models since there are many parameters to be adjusted and determined in water quality-ecological models [15,51,52,53,54]. Fortunately, values for these parameters have been reported [55,56,57,58,59,60] and can be adjusted and adopted in this study. The trial-and-error method for validating the water quality-ecological model is utilized.

Figure 6 compares the measured and simulated DO, DOC, NH₄, DON, PO₄, DOP, phytoplankton biomass, and zooplankton biomass. The simulated water quality-ecological state variables quantitatively agreed with the measured data. The model results depicted that the phytoplankton and zooplankton biomass was high from May to July.

Table 2. Statistical error based on differences between simulated and measured water quality-ecological state variables.

State Variable	MAE (Unit in State Variable)	CC	Skill
DO (mg/L)	0.22	0.87	0.92
DOC (mg C/L)	0.08	0.98	0.95
NH4 (mg/L)	0.02	0.86	0.90
DON (mg/L)	0.03	0.93	0.96
PO4 (mg/L)	0.002	0.99	0.91
DOP (mg/L)	0.005	0.96	0.96
Phytoplankton biomass (mg C/L)	0.16	0.99	0.95
Zooplankton biomass (mg C/L)	0.96	0.99	0.95

shows the statistical errors between the simulated and measured results for the water quality-ecological state variables. The predictive skill and correlation coefficient are all above 0.86. Excellent model performance is achieved when the skill value exceeds 0.90. The parameters adopted in the water quality-ecological model are displayed in the Supplementary Materials (Table S1).

4. Results and Discussion

4.1. Sensitivity Analysis of Phytoplankton and Zooplankton Biomass

Furthermore, the validated model was utilized to implement sensitivity analysis to explore how the parameters adopted in the planktonic module affected phytoplankton and zooplankton biomass. The parameters used for sensitivity analysis were GP, BMP, PRP, and WS for phytoplankton and GZ, BMZ, and PRZ for zooplankton. Even though basal metabolism is an inalienable and stable parameter of an organism, the parameter is also adjusted to understand how it affects the simulation result.

The model validation from September 2017 to July 2018 served as the baseline. Two alternative scenarios were employed to investigate the effect of planktonic parameters on the biomass of phytoplankton and zooplankton. They are the planktonic parameters plus 50% and minus 50% based on the parameters used for model validation.

Figure 7 illustrates the model sensitivity results for the most important parameters, which are the predation rate on phytoplankton (PRP) and basal metabolism rate of zooplankton (BMZ), affecting the phytoplankton biomass and zooplankton biomass from 11 July to 20 July 2018. When the PRP was reduced by 50%, the phytoplankton biomass was most affected (Figure 7a), while when the BMZ was decreased by 50%, the zooplankton biomass was most affected (Figure 7b).

Table 3. Results of sensitivity analysis for phytoplankton and zooplankton parameters.

Parameter	Condition	Phytoplankton (%)	Zooplankton (%)
Growth rate of phytoplankton (GP)	+50% −50%	8.67 −21.80	−0.08 0.51
Basal metabolism rate of phytoplankton (BMP)	+50% −50%	−0.73 0.72	0.03 −0.02
Predation rate on phytoplankton (PRP)	+50% −50%	−44.36 157.12	1.96 −5.12
Settling velocity of phytoplankton (WS)	+50% −50%	−0.13 0.13	−0.01 0.01
Growth rate of zooplankton (GZ)	+50% −50%	-28.05 55.03	28.43 −28.34
Basal metabolism rate of zooplankton (BMZ)	+50% −50%	52.10 −41.98	−28.74 62.54
Mortality rate of zooplankton (PRZ)	+50% −50%	0.34 −0.34	−0.26 0.26

Note: Minus and plus denote decreasing and increasing planktonic biomass, respectively.

lists the results of sensitivity runs for different planktonic parameters. This result shows that an increase in PRP and BMZ parameters results in a decrease in phytoplankton and zooplankton biomass. The maximum rates for increasing and decreasing phytoplankton biomass referred to the PRP parameter were 157.12% and 44.36%, respectively, while the maximum rates for increasing and decreasing zooplankton biomass related to the BMZ parameter were 62.54% and 28.74%, respectively.

Many parameters existed in the water quality-ecological model to be tuned. Most studies have shown the parameters only but have not implemented sensitivity analysis for parameters [15,30,61,62]. However, there are few parameters in the planktonic model, so the procedure for performing sensitivity analysis of parameters is not too cumbersome. According to the sensitivity analysis of parameters in the planktonic model, we found that the predation rate on phytoplankton (PRP) and basal metabolism rate of zooplankton (BMZ) were the most crucial parameters to be fine-tuned in the TFL.

4.2. Effect of Water Temperature and Inflow on Phytoplankton and Zooplankton Biomass

We suppose that environmental variations, including water temperature and inflow, would be dominant factors affecting planktonic biomass. Furthermore, the validated model was employed to probe the influence of environmental variables on phytoplankton and zooplankton biomass of TFL. Three scenarios were implemented to determine the main purpose of this section.

The first scenario was enacted to alter inflow discharge to cover a strong rainfall intensity-induced high inflow discharge into the lake. In this scenario run, the normal inflow was set to 0.1 m³/s, and inflow phytoplankton biomass, zooplankton biomass, and average water temperature were specified as 9.0 mg C/L, 35.1 mg C/L, and 15.2 °C. The artificial inflow hydrograph is illustrated in Figure 8. It clearly shows that the inflow discharge reached the maximum value (=6.3 m³/s) at the 127th. The second scenario run was set to alter the inflow water temperature to the summer temperature, which was 20.4 °C, while keeping the normal inflow discharge at 0.1 m³/s. The third scenario run was specified with the inflow water temperature using winter temperature, which was 7.8 °C and kept the same inflow discharge as the second scenario run. The model was implemented in a one-month simulation to yield the planktonic biomass.

Figure 9 depicts the comparison of calculated phytoplankton biomass at the buoy station under the three scenarios. The phytoplankton biomass was significantly affected by inflow discharge and reached a maximum of 6.7 mg C/L at the 134th hour in TFL but exhibited little effect by inflow water temperature. Figure 10 illustrates the distribution of phytoplankton biomass at the surface layer at the 134th hour. Figure 10a indicates that the phytoplankton biomass near the inflow was approximately 2.5 mg C/L, and the biomass at other positions was less than 2.0 mg C/L for the baseline case. Figure 10b shows that the biomass near the inflow reached 7.5 mg C/L, and the minimum value was approximately 3.5 mg C/L at the northeast side of TFL since the inflow discharge was changed. Figure 10c displays the phytoplankton biomass difference between the changing inflow discharge scenario run and the baseline run. This reveals that the phytoplankton biomass increased in the lake, and the biomass near the outflow location increased the most, by approximately 5.0 mg C/L.

Figure 11 delineates the comparison of calculated zooplankton biomass at the buoy station for the three scenario runs. This result indicates that the inflow discharge crucially affected the zooplankton biomass, reaching a maximum of 48.7 mg C/L at the 155th hour and then recovering to biomass of 24.0 mg C/L. It can be noted that the time for phytoplankton and zooplankton biomass to reach the maximum is different due to the fact that when phytoplankton flow into the lake, they serve as food for zooplankton predation, so the time for the maximum biomass of zooplankton is slightly later than that of phytoplankton. The zooplankton biomass under the changing inflow summer water temperature and changing inflow winter water temperature scenario runs was approximately 21.0 mg C/L and 24.0 mg C/L, respectively.

Figure 12 illustrates the distribution of zooplankton biomass at the surface layer at the 155th hour. Figure 12a shows that the zooplankton biomass near the inflow reaches 25.0 mg C/L for the baseline run, and the lower biomass at the outflow side is 14.0 mg C/L. Figure 12b indicates that the zooplankton biomass was higher at the inflow location, reaching 49.0 mg C/mg, and lower at the outflow location (=40.0 mg C/L) under the changing inflow discharge scenario run. Figure 12c shows the zooplankton biomass difference between the changing inflow discharge scenario and the baseline runs. This reveals that the zooplankton biomass increased, and the biomass at the center side of the lake increased the most. Overall, the range of increasing zooplankton biomass was from 8.1 to 29.1 mg C/L.

The light regime and grazing impact are known to be important drivers of phytoplankton biomass. The same applied to fish predators in the case of zooplankton biomass. However, there are many environmental factors affecting the planktonic biomass in lakes, such as water temperature, land use, and climate change [8,28,63,64,65]. Environmental change mostly results from anthropogenic impacts [1]. Land use and climate change alter inflow discharge, nutrient loads from watersheds, and water temperature, resulting in changes in ecological status. The results of sensitivity analysis for inflow discharge and thermal alterations indicate that inflow discharge alteration plays the most important role in affecting phytoplankton and zooplankton biomass in the TFL.

4.3. Limitations

The current ecological model did not consider intraspecific and interspecific competition in plankton assemblages that might lead to inconsistent results. Moreover, seasonal difference in planktonic communities would be taken into account since tropic structure affects the total biomass of plankton. In the ecological model, the module of fish biomass has been built, but the measured of fish species and biomass were scarce. Therefore, the fish module was not driven to simulate the fish biomass. In addition, the parameters that affect phytoplankton and zooplankton biomass vary in different lakes and must be adapted to local conditions. Therefore, when implementing the ecological model, it is best to have a good understanding of the characteristics of the parameters and to prioritize parameter testing and sensitivity analysis [66,67]. Incorrect selection of parameter values may lead to incorrect simulation results.

5. Conclusions

A coupling 3D hydrodynamic-ecological model called SCHISM-Ecol was developed and implemented to simulate planktonic biomass in Tsuei-Feng Lake (TFL), Taiwan. The observational data of water depth, water temperature, water quality, and phytoplankton and zooplankton biomass measured in 2017 and 2018 were utilized to validate the coupled model. Statistical metrics were used to quantitatively evaluate the modeling results and observed data. The modeling results reasonably reproduced the observations.

The sensitivity analysis was implemented using the validated model to determine which parameter was the most influential factor in dominating planktonic biomass. The results of sensitivity analysis indicated that the predation rate on phytoplankton (PRP) and basal metabolism rate of zooplankton (BMZ) were the crucial influential parameters for phytoplankton and zooplankton biomass, respectively. Moreover, the validated model was employed to prove the impact of environmental variables on planktonic biomass. The modeling results showed that inflow discharge was the most crucial factor subjected to changes in the external environment to affect the phytoplankton and zooplankton biomass of TFL.

Future work on lake studies will focus on the cross-sector cooperation that is needed to gather scientific researchers from different fields and to work on intensive field measurements to improve the understanding of the functioning of lake ecosystems. Another issue is to utilize the validated model to explore the effects of climate change on the water quality and planktonic biomass in the lake.