Development of an Ecosystem Model Considering Sediment Redox Processes in Enclosed Water Bodies

Abstract

This study investigates enclosed water bodies to understand the intricate interactions among physical, chemical, and biological processes. A one-dimensional model, integrating sediment and ecosystem components, was developed to analyze oxygen depletion, nutrient dynamics, and the influence of sediment organic matter on oxygen concentrations due to stratification. Oxygen and nutrient concentrations were derived by combining results from a one-dimensional vertical diffusion model with the ecosystem model. The dissolution, adsorption, and desorption of inorganic substances in bottom mud react differently under anaerobic and aerobic conditions, necessitating separate formulations for reduction and oxidation reactions in the bottom mud model. In Lake Biwa, Japan’s largest lake, dissolved oxygen levels have decreased despite regulatory efforts since the 1970s, prompting an investigation into the causes of oxygen-depleted water masses. The model was thus applied to examine oxygen consumption in Lake Biwa’s deep sediments, shedding light on anoxia in enclosed water bodies. This study emphasizes the significance of bottom sediments in water quality models and their implications for ecosystems. Furthermore, the adaptive nature of this model allows for its application in various water bodies, including freshwater reservoirs, lakes, and saltwater bays.

1. Introduction

The characteristic feature of enclosed water bodies is the intertwining of physical, chemical, and biological processes (phenomena related to flow, turbulence, temperature stratification, oxygen, organic matter, and organisms). In enclosed water bodies, freshwater from rivers generally enters from the surface, while seawater from the open ocean enters from the lower layers. Consequently, in areas where there is a significant load of organic matter from land and marine aquaculture, the lake or sea waters stagnate due to weather and topographical conditions. This often leads to stratification in the marine environment, especially during the high-temperature period in summer, resulting in the formation of oxygen-depleted water masses in the bottom layer. If there are harbor trenches or dredged areas in these regions, the formation rate of oxygen-depleted water masses significantly increases. In lakes such as Lake Biwa, Ikeda Lake, and marine areas like Ise Bay, Mikawa Bay, and the eastern part of Osaka Bay, as well as the Seto Inland Sea, red tide occurs in the surface layer, and oxygen-depleted water masses occur in the bottom layer during the summer, causing significant damage to ecosystems [1,2,3]. Causes of oxygen-depleted water masses in enclosed water bodies include nutrient loading from land, organic pollution associated with aquaculture in stagnant water areas, and endogenous/exogenous organic pollution derived from vegetation [4,5,6]. Weather changes due to global warming and alterations in the use of river water also lead to a decrease in the supply of low-temperature, high-oxygen-concentration river water and a decrease in the utilization of river water [7].

Prolonged oxygen deficiency in the bottom layer environment results in an increase in anaerobic organisms, an increase in hydrogen sulfide and sulfide production by sulfate-reducing bacteria, and the release of nutrients from the sediment. The released nutrients are transported to the euphotic layer, utilized for plankton growth, organic matter is produced, and this organic matter eventually settles into the substrate, generating the sediment. Such interactions between water and substrate have a significant impact on water quality in enclosed water bodies where outflow loads are expected to be high [8,9,10].

In the past years, for addressing such challenges, mathematical approaches combining ecosystem models such as water quality prediction models with detailed water environment monitoring have proven effective [11,12,13,14,15]. Among them, ecosystem models consider biochemical and ecological processes in material cycling in water areas, modeling the interactions between organisms and environmental factors, enabling us to propose better water management plans and restoration programs through better predictive capabilities.

Accurate long-term prediction of water quality changes requires models that simultaneously solve for water and sediment nutrient conditions [16,17,18,19,20,21]. Some models of sediments and various water quality models that include them have already been developed [22,23,24,25]. Chapelle et al. [26] and Tamsalu et al. [27] have modeled nitrogen and phosphorus in sediments and combined them with water quality models, focusing on Vilainc Bay and Finland Bay, respectively. While these substrate models are based on the material cycling of nitrogen and phosphorus in sediments, Baretta et al. [28] and Ruardij and Raaphorst [29] targeted the North Sea, and Horie [30] targeted Kure Bay. They have proposed models that directly express the dynamics of benthic organisms, delving into the ecosystem of intertidal zones and the impact on benthic organisms in eutrophication measures. Considering the objectives, these models have been extensively studied, such as the analysis of the ecosystem of tidal flats and the response of water quality and nutrient efflux fluxes associated with load reduction in the inner bay areas.

In Lake Biwa, which has the largest lake surface area in Japan, oxygen-depleted water masses in the deep lake bottom have been noticeable in recent years. Until the 1970s, the annual minimum dissolved oxygen concentration at the bottom of Lake Biwa decreased due to eutrophication associated with increased organic material and nutrient loading from the land. Even though organic material and nutrient loading from the land were regulated by the early 1980s, dissolved oxygen concentrations continued to decline, sometimes falling below the general criterion for hypoxia, 2 mg/L. One of the causes is attributed to recent climate change [1], but the mechanisms by which interactions in physical, chemical, and biological processes induce hypoxia are not clearly understood. Ohte et al. [31] analyzed the dissolved oxygen dynamics in the deep layer of the northern basin of Lake Biwa during stratification periods based on oxygen isotope ratios, estimating that oxygen consumption by lake sediment accounted for approximately 60% of total oxygen consumption, and that by the water column accounted for about 40%. This result indicates that oxygen consumption by lake sediment is the dominant factor in dissolved oxygen consumption in the deep layer of Lake Biwa’s northern basin. Therefore, understanding the oxygen consumption dynamics in lake sediment in this area is essential for elucidating the formation process of oxygen-depleted water masses in the deep lake bottom of the northern basin.

In this study, a bottom sediment model is constructed to reproduce nutrient sediment remobilization and bottom sediment and oxygen consumption in the bottom sediment in order to elucidate anoxic phenomena in an enclosed water body. The model is applied to Lake Biwa, which is an anoxic lake, and its characteristics are discussed. The developed model is applicable to all enclosed water bodies of reservoirs and lakes, which are freshwater cases, and can also be applied to enclosed bays, which are saltwater cases.

2. Methodology

2.1. Simulation Model

In this study, the new vertical one-dimensional model, based on the three-dimensional hydrodynamic model [32] and ecosystem model built by Koue et al. [33], was developed. The new ecosystem model considered 11 state variables, including “phytoplankton (Chlorophyll a, b, c)”, “zooplankton”, “suspended organic matter”, “dissolved organic matter”, “ammonium nitrogen (NH₄-N)”, “nitrite nitrogen (NO₂-N)”, “nitrate nitrogen (NO₃-N)”, “phosphate phosphorus (PO₄-P)”, “dissolved oxygen (DO)”, “total nitrogen (T-N)”, and “total phosphorus (T-P)” and calculated concentration changes. Additionally, to replicate oxygen consumption in sediment, a sediment model was added as a sub model in this study.

The constructed model was applied to the northern observation point of Lake Biwa, the largest freshwater lake in Japan, located at Imazu-oki (35°23′41″ N, 136°07′57″ E), and its reproducibility was evaluated. Lake Biwa has a surface area of 674 km², an average depth of 41 m, and a maximum depth of 103 m near Imazu-oki. The total inflow of river water, groundwater contributions, and direct surface precipitation into Lake Biwa is estimated to be approximately 5.89 billion m³ per year. Given that the storage capacity of Lake Biwa is 27.5 billion m³, the calculated retention time stands at approximately 4.7 years.

2.1.1. Computational Domain

Figure 1 shows the simulation domain of this model. At the Imazu-oki monitoring point, an observation station in the northern part of Lake Biwa, a vertical unequal mesh was installed from the surface to a depth of about 90 m in order to evaluate the water temperature stratification and bottom sediment environment in detail. From the surface to a depth of 20 m, the grid size was set at 0.5 m. From a depth of 20 m to 55 m, the grid size was increased from 0.5 m to a maximum of 2.5 m in 0.1 m increments. From 55 m to 87 m depth, the grid size was decreased from 2.5 m by 0.1 m to a minimum of 0.1 m. From a depth of 87 m and below to the bottom of the lake at a depth of 90 m, the grid size was set at 0.1 m. The vertical axis (z-axis) was set upwards.

2.1.2. Computational Period

The calculation period was set from 1 April 2006 to 31 March 2014. The first year was designated as a warm-up period, and the subsequent 7 years (from 1 April 2007 to 31 March 2014) were considered the evaluation period. The time step was set to 1 s.

2.2. Hydrodynamics Model

The one-dimensional model was constructed based on the three-dimensional hydrodynamic model of Lake Biwa developed by Koue et al. [32]. The details are as follows.

2.2.1. Basic Equations

The one-dimensional diffusion model used in this study consisted of a vertical one-dimensional diffusion equation, a hydrostatic pressure approximation equation, and a diffusion equation for water temperature. The basic equations are shown below.

Vertical one-dimensional diffusion equation is described as follows:

$$\frac{\partial u}{\partial t}={\nu }_z\frac{{\partial }^{2}u}{\partial z^2}$$

(1)

$$\frac{\partial v}{\partial t}={\nu }_z\frac{{\partial }^{2}v}{\partial z^2}$$

(2)

Hydrostatic pressure approximation is given as follows:

$$0=-\frac{1}{{\rho }_0}\frac{\partial p}{\partial z}-\frac{\rho }{{\rho }_0}g$$

(3)

Heat diffusion equation is described as follows:

$$\frac{\partial T}{\partial t}={\kappa }_z\frac{{\partial }^{2}T}{\partial z^2}$$

(4)

Density is calculated as follows:

$$\rho =1000.07-0.00469{\left(T-273.15\right)}^2-0.0035\left(T-273.15\right)$$

(5)

Vertical eddy viscosity and eddy diffusion coefficient are calculated as follows:

$$Ri=-\frac{g}{{\rho }_0}·\frac{\frac{\partial \rho }{\partial z}}{{\left(\frac{\partial W}{\partial z}\right)}^2}$$

(6)

$${\nu }_z=\frac{0.0001}{1.0+5.2Ri}$$

(7)

$${\kappa }_z=\frac{0.0001}{{\left(1.0+\frac{10}{3}\times Ri\right)}^{\frac{3}{2}}}$$

(8)

u, v: horizontal velocity (m/s), W: Magnitude of horizontal flow velocity (m/s), T: water temperature (K), p: pressure (N/m²), $\rho$: Density of lake water (kg/m³), ${\rho }_0$= 103: Density of reference lake water (kg/m³), g = 9.8: Gravitational acceleration (m/s²), νz: Vertical eddy viscosity (m²/s), κ_z: Vertical eddy diffusivity (m²/s).

In summer, the thermocline is formed in Lake Biwa at a depth of 10 to 20 m. Therefore, ${\mathsf{\nu }}_{\mathrm{z}}$ and ${\mathsf{\kappa }}_{\mathrm{z}}$ are obtained using Richardson number Ri, a dimensionless number that represents the instability of the flow where the density is changing in a layered manner.

2.2.2. Meteorological Data

The meso-scale model grid point value (MSM-GPV) data, produced by Japan Meteorological Agency, were used for heat balance on the lake surface and wind-induced surface shear stress. The data items used include “temperature (K)”, “pressure (Pa)”, “wind direction and speed (m/s)”, and “relative humidity (%)”. For solar radiation, observations from the Hikone Regional Meteorological Observatory in Figure 1 were used, provided at hourly intervals.

Figure 2 shows the position of GPV data within the calculation domain of the Lake Biwa watershed model. The grid spacing in GPV data is 5 km, and in this study, the meteorological data at Imazu-oki used a weighted average inverse proportional to the square of the distance to obtain interpolated values. These meteorological data were updated at hourly intervals and used in the calculation of the hydrodynamic model. Hourly updates enable a more detailed reproduction of heat balance and flow on the lake surface.

2.2.3. Initial Conditions

The initial condition was set to a flow velocity of 0 m/s. The initial water temperature on 1 April 2006 was set by linear interpolation of data from 20 March 2006 and 10 April 2006, which were observed twice a month by the Lake Biwa Environmental Research Institute at Imazu-oki at a depth of 0.5 m, 5 m, 10 m, 15 m, 20 m, 30 m, 40 m, 60 m, 80 m, and approximately 90 m, out of data observed at 10 sampling point locations.

2.3. Water Quality Ecosystem Model

The one-dimensional model for long-term analysis was created based on the three-dimensional water quality model of Lake Biwa developed by Koue et al. [33]. Furthermore, in order to analyze inorganic nutrient dynamics and dissolved oxygen consumption in detail, inorganic nitrogen was divided into the state variables of “ammonia-form nitrogen”, “nitrite-form nitrogen”, and “nitrate-form nitrogen” and inorganic phosphorus into “phosphate-form phosphorus.” The vertical flow field was solved by a one-dimensional diffusion equation, and using the results of the vertical one-dimensional diffusion model, a new model was developed by adding a bottom sediment model that considers anaerobic and aerobic processes in bottom mud sediments to a suspended ecosystem model that considers three types of phytoplankton.

2.3.1. Basic Equations

A vertical one-dimensional model is used in the water quality calculations to determine the average concentration of each layer. Each constituent is affected by changes due to chemical and biological processes and by diffusion. The time variation is represented by the following one-dimensional diffusion equation:

$$\frac{\partial C_i}{\partial t}={\kappa }_z\frac{{\partial }^2{C}_i}{\partial z^2}+Q_{C_i}$$

(9)

The left side represents the time-varying term and the right side represents the diffusion term and changes due to chemical and biological processes ($Q_{C_i}$). The state variables are the 11 components shown in

Table 1. State Variables in Ecosystem Model.

Element	Symbol	Unit
Phytoplankton	$C_{PX}$	μgC/L
Zooplankton	$C_Z$	μgC/L
Suspended organic matter	$C_{POC}$	mg/m3
Dissolved organic matter	$C_{DOC}$	mg/m3
Ammonia nitrogen	$C_{N{H}_{4}N}$	mg/m3
Nitrite nitrogen	$C_{N{O}_{2}N}$	mg/m3
Nitrate nitrogen	$C_{N{O}_{3}N}$	mg/m3
Phosphate-form phosphorus	$C_{P{O}_{4}P}$	mg/m3
Dissolved oxygen	$C_{DO}$	mg/l
Total nitrogen	$C_{TN}$	mg/m3
Total phosphorus	$C_{TP}$	mg/m3
Total carbon	$C_{TC}$	mg/m3

. A summary of the chemical and biological processes for each state variable is shown in

Table 2. Concentration balance for each state variable in the ecosystem model.

$C_{PX}$	Growth by photosynthesis—Extracellular secretion—Predation by zooplankton—Phytoplankton respiration—Dying—Sedimentation
$C_Z$	Growth by phytoplankton feeding—Zooplankton feces—Zooplankton excretion by respiration—Zooplankton death
$C_{POC}$	Phytoplankton die-off + zooplankton feces + natural mortality of zooplankton—mineralization of suspended-form organic matter—generation of decomposition surplus + sedimentation of suspended-form organic matter
$C_{DOC}$	Extracellular secretion + decomposition surplus generation—mineralization of dissolved organic matter by aerobic bacteria
$C_{N{H}_{4}N}$	-Phytoplankton ingestion by photosynthesis + phytoplankton respiration + zooplankton excretion by respiration + mineralization of suspended organic matter + mineralization of dissolved organic matter by aerobic bacteria + leaching of nitrogen from bottom sediment—nitrification of ammonia-form nitrogen
$C_{N{O}_{2}N}$	Nitrification of ammonia nitrogen—Nitrite nitrogen
$C_{N{O}_{3}N}$	-Phytoplankton uptake by photosynthesis + nitrification and denitrification of nitrite nitrogen
$C_{P{O}_{4}P}$	-Phytoplankton uptake by photosynthesis + phytoplankton respiration + zooplankton excretion by respiration + mineralization of suspended organic matter + mineralization of dissolved organic matter by aerobic bacteria + leaching of phosphorus from bottom sediment
$C_{DO}$	Supply by photosynthesis—Consumption by phytoplankton respiration—Consumption by zooplankton respiration—Oxygen consumption by suspended organic matter—Oxygen consumption by dissolved organic matter—Consumption in bottom mud + aeration—Oxygen consumption associated with nitrification of ammonia-form nitrogen—Oxygen consumption associated with nitrite-form nitrogen
$C_{TN}$	N/C ratio in phytoplankton x time variation of phytoplankton + N/C ratio in zooplankton × time variation of zooplankton + N/C ratio in suspended organic matter × time variation of suspended organic matter + N/C ratio in dissolved organic matter × time variation of dissolved organic matter + time variation of ammonia-form nitrogen + time variation of nitrite-form nitrogen + nitrate-form nitrogen
$C_{TP}$	P/C ratio in phytoplankton × time variation of phytoplankton + P/C ratio in zooplankton × time variation of zooplankton + P/C ratio in suspended organic matter × time variation of suspended organic matter + P/C ratio in dissolved organic matter × time variation of dissolved organic matter + time variation of phosphate form phosphorus
$C_{TC}$	Time variation of phytoplankton + time variation of zooplankton + time variation of suspended organic matter + time variation of dissolved organic matter

.

2.3.2. Initial Conditions

The data for water temperature, flow direction, and vertical flow velocity used in the ecosystem model were updated every 3 h from the calculation results of the vertical one-dimensional diffusion model. Initial conditions for nutrients such as nitrogen and phosphorus were set based on the data observed twice a month at Imazu-oki by the Lake Biwa Environmental Research Institute.

The amount of change in chemical and biological processes ($Q_{C_i}$) for each state variable was calculated as follows. The variables in each equation are shown in

Table 3. List of symbols of ecosystem model.

${\alpha }_{1-13}$	Velocity coefficient
${\mathsf{\beta }}_{1-13}$	Temperature coefficient
$K_{N_X}$	Half-saturation constant of DO for nitrogen
$K_{P_X}$	Half-saturation constant of DO for phosphorus
$\epsilon$	Light Attenuation Rate
${I_X}_{opt}$	Optimum light intensity
${T_X}_{opt}$	Optimum temperature
${\left[Chl.a:C\right]}_{pp}$	Chl.a/C ratio in phytoplankton
${\left[Chl.b:C\right]}_{pp}$	Chl.b/C ratio in phytoplankton
${\left[Chl.c:C\right]}_{pp}$	Chl.c/C ratio in phytoplankton
$\mathsf{\Omega }$	Threshold concentration of food
$\lambda$	Ivlev constant
$V{M}_{px}$	Sedimentation rate of plankton
$e$	Assimilation efficiency
$g$r	Total growth rate
$D{O}_{1,4,5}$	Half-saturation value of DO
$D{O}_{2,3}$	Half saturation value of oxygen limit
$D{O}_6$	DO threshold for denitrification
$\zeta$	Percentage of decomposition surplus generated
$V{M}_{POC}$	Sedimentation rate of POC
${\left[P:C\right]}_{PP}$	P/C ratio in phytoplankton
${\left[P:C\right]}_{ZP}$	P/C ratio in zooplankton
${\left[P:C\right]}_{POC}$	P/C ratio in suspended organic matter
${\left[P:C\right]}_{DOC}$	P/C ratio in dissolved organic matter
${\gamma }_P$	Inhibition of phosphorus leaching by DO
${\gamma }_1$	NO3N/NH4N distribution ratio in inorganic nitrogen absorption
${\left[N:C\right]}_{PP}$	N/C ratio in phytoplankton
${\left[N:C\right]}_{ZP}$	N/C ratio in zooplankton
${\left[N:C\right]}_{POC}$	N/C ratio in suspended organic matter
${\left[N:C\right]}_{DOC}$	N/C ratio in dissolved organic matter
${\gamma }_N$	Inhibition of nitrogen leaching by DO
${\left[TOD:C\right]}_{PP}$	TOD/C ratio in phytoplankton
${\left[TOD:C\right]}_{ZP}$	TOD/C ratio in zooplankton
${\left[TOD:C\right]}_{POC}$	TOD/C ratio in suspended organic matter
${\left[TOD:C\right]}_{DOC}$	TOD/C ratio in dissolved organic matter
$k_a$	Reaeration coefficient

.

• Phytoplankton;

In this model, phytoplankton is divided into chlorophyll a, b, c to simplify handling. The time variation of phytoplankton concentration is expressed by the following equation:

$$\frac{d{C}_{px}}{dt}=B1-B2-B3-B4-B5-B6$$

(10)

where B1 is growth by photosynthesis, B2 is extracellular secretion, B3 is feeding by zooplankton, B4 is respiration, B5 is mortality, and B6 is the amount of change due to sedimentation.

The amount of increase B1 is given by the growth rate of phytoplankton as a function of water temperature T (°C), and by considering the control function min(f(N), f(P)) for nutrients, the control function f(I) for light intensity I (W/m²) in the water, and the control function f(T) for temperature, the following formula is used.

$$B1={\alpha }_1·\exp \left({\beta }_1·T\right)·\left\{min\left(f\left(N\right),f\left(P\right)\right)·f\left(I\right)·f\left(T\right)\right\}·C_{px}$$

(11)

The control function for nutrients is expressed as follows:

$$min\left(f\left(N\right),f\left(P\right)\right)=\min \left(\frac{C_{N{H}_{4}N}+C_{N{O}_{2}N}+C_{N{O}_{3}N}}{K_{N_X}+C_{N{H}_{4}N}+C_{N{O}_{2}N}+C_{N{O}_{3}N}},\frac{C_{P{O}_{4}P}}{K_{P_X}+C_{P{O}_{4}P}}\right)$$

(12)

The control function f(I) for light intensity in water is expressed as follows:

$$f\left(I\right)=\frac{light}{{I_X}_{opt}}·\exp \left(1.0-\frac{light}{{I_X}_{opt}}\right)$$

(13)

where light is the attenuation of light at a depth of zh(m) and is expressed by the following equation:

$$\mathit{light}=I·\mathit{exp}(-\epsilon ·\mathit{zh})$$

(14)

The control function of temperature, f(T), is expressed as follows:

$$f\left(T\right)=\frac{T}{{T_X}_{opt}}·\exp \left(1.0-\frac{T}{{T_X}_{opt}}\right)$$

(15)

Extracellular secretion B2, a phenomenon in which a part of organic substances fixed by photosynthesis is discharged out of the cell in the form of dissolved organic substances, is expressed by the following equation [34,35]:

$$B2=0.135·\exp \left(-0.00201·{\left[Chl.a:C\right]}_{pp}·C_{PX}\right)·B1$$

(16)

The zooplankton diet B3 is expressed using the Ivlev equation [36], which takes into account the existence of a threshold [37] for the concentration of filter-feeding (phytoplankton) of filter-feeding zooplankton, as follows:

$$B3={\alpha }_2·\exp \left({\beta }_2·T\right)[1.0-\exp \left\{\lambda \left(\mathsf{\Omega }-C_{PX}\right)\right\}·C_Z$$

(17)

where $\mathsf{\Omega }$ (mg/m³) is the threshold value for food concentration, and B3 = 0 when $\mathsf{\Omega }\ge C_{PX}$.

Respiration rate B4 and mortality rate B5 are expressed by the following equation considering temperature dependence.

$$B4={\alpha }_3·exp\left({\beta }_3·T\right)·C_{PX}$$

(18)

$$B5={\alpha }_4·exp\left({\beta }_4·T\right)·C_{PX}$$

(19)

The phytoplankton sedimentation rate B6 is given by the following equation using its sedimentation rate $V{M}_{PX}$ (m/d), phytoplankton concentration $Px\left(k\right)$ (mg/m³) at the position of vertical coordinate k, and distance $dZ\left(k\right)$(m) from the lake bottom.

$$B6=V{M}_{PX}·\frac{C_{PX}\left(k+1\right)}{dZ\left(k+1\right)}-V{M}_{PX}·\frac{C_{PX}\left(k\right)}{dZ\left(k\right)}$$

(20)

2.

Zooplankton;

The time variation of phytoplankton concentration is expressed by the following equation:

$$\frac{d{C}_z}{dt}=B3-B7-B8-B9$$

(21)

where B7 is the amount of feces, B8 is respiration, and B9 is the amount of variation due to natural mortality. The amount of feces excreted B7 and respiration B8 are expressed by the following equation using the zooplankton feeding rate B3 and the indices of assimilation efficiency e and total growth efficiency gr

$$B7=\left(1-e\right)·B3$$

(22)

$$B8=\left(e-gr\right)·B3$$

(23)

Natural mortality of zooplankton B9 is assumed to depend primarily on water temperature and dissolved oxygen and is given by:

$$B9={\alpha }_5·\exp \left({\beta }_{5}T\right)·\frac{D{O}_1+C_{DO}}{C_{DO}}·C_Z$$

(24)

3.

Suspended organic matter;

Suspended organic matter is granular organic matter (detritus) in the non-living state, and its concentration over time is expressed by the following equation:

$$\frac{d{C}_{POC}}{dt}=B5+B7+B9-B10-B11-B12$$

(25)

where B10 represents the mineralization of suspended organic matter, B11 represents the generation of decomposition surplus, and B12 represents the amount of variation due to sedimentation of suspended organic matter.

The mineralization of organic matter corresponds to respiration by aerobic bacteria, and is primarily a function of water temperature and dissolved oxygen concentration, and is given by the following equation:

$$B10={\alpha }_6\exp \left({\beta }_{6}T\right)·\frac{C_{DO}}{D{O}_2+C_{DO}}·C_{POC}$$

(26)

The formation of decomposition surplus is a phenomenon in which persistent suspended organic matter, such as humic substances, is not mineralized but is converted to dissolved organic matter. The amount of variation in this process, B11, is expressed by the following equation using the ratio of surplus production to mineralization of suspended organic matter, ζ:

$$B11=\zeta ·B10$$

(27)

The sedimentation rate B12 is given by the following equation using its sedimentation rate VM_POC (m/d) and the concentration of suspended organic matter C_POC (k) (mg/m³) at the vertical coordinate k

$$B12=V{M}_{POC}·\frac{C_{POC}\left(k+1\right)}{dZ\left(k+1\right)}-V{M}_{POC}·\frac{C_{POC}\left(k\right)}{dZ\left(k\right)}$$

(28)

4.

Dissolved organic matter;

The time variation of dissolved organic matter is expressed by the following equation:

$$\frac{d{C}_{DOC}}{dt}=B2+B11-B13$$

(29)

where B13 represents the amount of consumption due to mineralization of dissolved organic matter by aerobic bacteria, and is expressed as the following equation based on the dependence on water temperature and dissolved oxygen concentration [38].

$$B13={\alpha }_7\exp \left({\beta }_{7}T\right)·\frac{C_{DO}}{D{O}_3+C_{DO}}·C_{DOC}$$

(30)

5.

Phosphate phosphorus;

The time variation of phosphate-form phosphorus is expressed by the following equation:

$$\begin{gathered} \frac{\partial C_{P{O}_{4}P}}{\partial t}=-{\left[P:C\right]}_{PX}B1+{\left[P:C\right]}_{PX}B4+{\left[P:C\right]}_{ZP}B8+{\left[P:C\right]}_{POC}B10 \\ +{\left[P:C\right]}_{DOC}B13+B14 \end{gathered}$$

(31)

where B14 represents the amount of variation due to leaching of phosphorus from the bottom sediment.

$$B14=\frac{{\alpha }_8\exp \left({\beta }_{8}T-{\gamma }_P·C_{DO}\right)}{dZ\left(k\right)}$$

(32)

6.

Ammonia nitrogen;

The time variation of ammonia nitrogen is expressed by the following equation:

$$\begin{gathered} \frac{\partial C_{NH4N}}{\partial t}=-\left(1-{\gamma }_1\right){\left[N:C\right]}_{PX}B1+{\left[N:C\right]}_{PX}B4+{\left[N:C\right]}_{ZP}B8+{\left[N:C\right]}_{POC}B10 \\ +{\left[N:C\right]}_{DOC}B13+B15-B16 \end{gathered}$$

(33)

where B15 represents the leaching of nitrogen (ammonia) from the bottom sediment and B16 represents the variation of ammonia nitrogen by nitrification.

$$B15=\frac{{\alpha }_9\exp \left({\beta }_{9}T-{\gamma }_N·C_{DO}\right)}{dZ\left(k\right)}$$

(34)

The process of ammonia digestion is related to the activity of nitrifying bacteria and is expressed by the following equation considering the dependence on water temperature and dissolved oxygen concentration:

$$B16={\alpha }_{10}\exp \left({\beta }_{10}T\right)\frac{C_{DO}}{D{O}_4+C_{DO}}{C}_{N{H}_{4}N}$$

(35)

7.

Nitrite nitrogen;

The time variation of nitrite nitrogen is expressed by the following equation:

$$\frac{\partial C_{N{O}_{2}N}}{\partial t}=B16-B17$$

(36)

where B17 represents the amount of variation of nitrite nitrogen due to nitrification.

$$B17={\alpha }_{11}\exp \left({\beta }_{11}T\right)\frac{C_{DO}}{D{O}_5+C_{DO}}{C}_{N{O}_{2}N}$$

(37)

8.

Nitrate nitrogen;

The time variation of nitrate nitrogen is expressed by the following equation:

$$\frac{\partial C_{N{O}_{3}N}}{\partial t}=-{\gamma }_1{\left[N:C\right]}_{PP}B1+B17-B18$$

(38)

where B18 represents the amount of variation due to denitrification. Denitrification is a biochemical process in which nitrate nitrogen is reduced to nitrogen gas under anaerobic conditions, and by setting $D{O}_6$ (mg/L) as the limit value of dissolved oxygen concentration at which denitrification occurs, the denitrification amount, B18, can be expressed by the following formula.

$$\begin{gathered} B18={\alpha }_{12}\exp \left({\beta }_{12}T\right)C_{N{O}_{3}N} \left(C_{DO}\le D{O}_6\right) \\ 0 \left(C_{DO}>D{O}_6\right) \end{gathered}$$

(39)

9.

Dissolved oxygen concentration;

The time variation of dissolved oxygen concentration is expressed by converting the concentration variation due to all biochemical reaction terms related to the production and consumption of oxygen into the amount of oxygen.

$$\begin{gathered} \frac{\partial C_{DO}}{\partial t}={\left[TOD:C\right]}_{PX}B1-{\left[TOD:C\right]}_{PX}B4-{\left[TOD:C\right]}_{ZP}B8-{\left[TOD:C\right]}_{POC}B10 \\ -{\left[TOD:C\right]}_{DOC}B13-B19+B20-B21-B22 \end{gathered}$$

(40)

where B19 is the oxygen consumption by the bottom sediment, B20 is the supply by reaeration, and B21 and B22 are the oxygen consumption associated with nitrification of ammonia nitrogen and nitrite nitrogen, respectively.

The oxygen consumption of bottom sediment, B19, as well as the leaching flux of nutrients (nitrogen and phosphorus), is expressed as follows, taking into account the temperature dependence.

$$B19=\frac{{\alpha }_{13}\exp \left({\beta }_{13}\left(T-T_B\right)\right)}{dZ\left(k\right)}\times {10}^{-3}$$

(41)

The oxygen supply by reaeration, B20, is expressed by the following equation using the saturated dissolved oxygen $D{O}_s$ (mg/L) and the reaeration coefficient $k_a\left(d^{-1}\right)$.

$$B20=k_a\left(D{O}_s-C_{DO}\right)$$

(42)

$$D{O}_s=\frac{32DOS{S}_S}{22.4\left(1+T/273\right)}$$

(43)

$$DOS{S}_s=10.291-0.2809T+0.006009T^2-0.000063T^3$$

(44)

The oxygen consumption B21 and B22 due to nitrification of ammonia-form nitrogen and nitrite-form nitrogen are expressed by the following equation:

$$B21=-3.429\times {10}^{-3}\times B16$$

(45)

$$B22=-1.1429\times {10}^{-3}\times B17$$

(46)

10.

Total nitrogen, total phosphorus, and total organic carbon;

The time variation of total nitrogen, total phosphorus, and total organic carbon content is expressed as the sum of the time variation of the state variables related to these variables.

$$\begin{gathered} \frac{\partial C_{TN}}{\partial t}={\left[N:C\right]}_{PX}·\left(\frac{d{C}_{PX}}{dt}\right)+{\left[N:C\right]}_{ZP}·\left(\frac{d{C}_{ZP}}{dt}\right)+{\left[N:C\right]}_{POC}·\left(\frac{d{C}_{POC}}{dt}\right) \\ +{\left[N:C\right]}_{DOC}·\left(\frac{d{C}_{DOC}}{dt}\right)+\frac{d{C}_{N{H}_{4}N}}{dt}+\frac{d{C}_{N{O}_{2}N}}{dt}+\frac{d{C}_{N{O}_{3}N}}{dt} \end{gathered}$$

(47)

$$\begin{gathered} \frac{\partial C_{TP}}{\partial t}={\left[P:C\right]}_{PX}·\left(\frac{d{C}_{PX}}{dt}\right)+{\left[P:C\right]}_{ZP}·\left(\frac{d{C}_{ZP}}{dt}\right)+{\left[P:C\right]}_{POC}·\left(\frac{d{C}_{POC}}{dt}\right) \\ +{\left[P:C\right]}_{DOC}·\left(\frac{d{C}_{DOC}}{dt}\right)+\frac{d{C}_{P{O}_{4}P}}{dt} \end{gathered}$$

(48)

$$\frac{\partial C_{TC}}{\partial t}=\frac{d{C}_{PX}}{dt}+\frac{d{C}_{ZP}}{dt}+\frac{d{C}_{POC}}{dt}+\frac{d{C}_{DOC}}{dt}$$

(49)

2.4. Sediment Model

A bottom sediment model was constructed to model the behavior of nitrogen, phosphorus, and organic matter in bottom sediments. The state variables used in the bottom sediment model are listed in

Table 4. Sediment Model State Variables.

Element		Symbol	Unit
Sediment	Inorganic nitrogen	$C_{MIN}$	mg/m3 dry mud
	Inorganic phosphorus	$C_{MIP}$	mg/m3 dry mud
	Organic nitrogen	$C_{MON}$	mg/m3 dry mud
	Organic phosphorus	$C_{MOP}$	mg/m3 dry mud
	Organic matter	$C_{MOC}$	mg/m3 dry mud
Pore water	Ammonia nitrogen	$C_{WN{H}_{4}N}$	mg/m3
	Nitrite nitrogen	$C_{WN{O}_{2}N}$	mg/m3
	Nitrate-nitrogen	$C_{WN{O}_{3}N}$	mg/m3
	Inorganic phosphorus	$C_{WIP}$	mg/m3
	Dissolved organic matter	$C_{WOC}$	mg/m3
	Dissolved oxygen	$C_{WDO}$	mg/L

. The chemical and biological processes of each state variable are summarized in

Table 5. Concentration balance for each state variable in the sediment model.

$C_{MIN}$	Sedimentation + decomposition of organic nitrogen in bottom mud (aerobic/anaerobic) + adsorption/desorption onto mud particles
$C_{MIP}$	Sedimentation + decomposition of organic phosphorus in bottom mud + adsorption/desorption on mud particles—adsorption/desorption on pore water (aerobic/anaerobic)
$C_{MON}$	Sedimentation—Decomposition of organic nitrogen in bottom mud (aerobic/anaerobic)
$C_{MOP}$	Sedimentation—Decomposition of organic phosphorus in bottom mud
$C_{MOC}$	Sedimentation—Dissolution of organic matter in bottom sediment
$C_{WN{H}_{4}N}$	Decomposition of organic nitrogen in bottom mud—Adsorption/desorption on mud particles + Diffusion—Nitrification of ammonia nitrogen
$C_{WN{O}_{2}N}$	Diffusion + Nitrification of ammonia nitrogen—Nitrification of nitrite nitrogen
$C_{WN{O}_{3}N}$	Diffusion + Nitrification of nitrite nitrogen—Denitrification of nitrate nitrogen
$C_{WIP}$	Decomposition of organic phosphorus in bottom mud—Adsorption/desorption on mud particles + diffusion + desorption from bottom mud (aerobic/anaerobic)
$C_{WOC}$	Dissolution of organic matter in sludge + diffusion—Decomposition of organic matter in pore water
$C_{WDO}$	Diffusion—Oxygen consumption by nitrification of ammonia nitrogen—Oxygen consumption by nitrite nitrogen—Oxygen consumption by decomposition of organic matter in pore water

.

The amount of change in chemical and biological processes ($Q_{C_i}$) for each state variable was calculated as follows. The variables in each equation are shown in

Table 6. List of symbols of the sediment model.

$S_{IN}$	Deposition rate of inorganic nitrogen
$D_{SN}$	Decomposition rate of organic nitrogen
${\theta }_{N,P,DEO,DEA,C}^{\left(T-20\right)}$	Temperature coefficient
$D{O}_{SON}$	DO half-saturation constant in oxidative decomposition of nitrogen
$\varphi$	porosity
$C_{EXN}$	EXN (exchangeable nitrogen)/N ratio in bottom sediment
$T_{AD}$	Time to adsorption equilibrium
$S_{IP}$	Deposition rate of inorganic phosphorus
$D_{SP}$	Degradation rate of organic phosphorus
$C_{EXP}$	EXP (exchangeable phosphorus)/P ratio in bottom sediment
$D{O}_{CR}$	DO threshold for aerobic and anaerobic (phosphorus adsorption)
$D{E}_O$	Desorption rate constant (under aerobic conditions)
$D{E}_A$	Desorption rate constant (under anaerobic conditions)
$S_{ON}$	Organic nitrogen deposition rate
$S_{OP}$	Organic phosphorus deposition rate
$S_{OC}$	Organic matter deposition rate
$D_{SC}$	Organic dissolution rate
DC	Diffusion coefficient in pore water
NR	Nitrification rate constant for NH4N
$D{O}_{NR,NR2,DNR}$	Half-saturation constant for DO concentration
NR2	Nitrification rate constant for NO2N
DNR	Denitrification rate constant for NO3N
$N{H}_4{N}_{DO}$	Oxygen consumption for nitrification of NH4N
$N{O}_2{N}_{DO}$	Oxygen consumption during nitrification of NO2N
$C_{DO}$	Oxygen consumption by decomposition of organic matter in pore water

.

• Inorganic nitrogen in bottom sediment;

The time variation of inorganic nitrogen in bottom sediment is expressed as follows:

$$\frac{d{C}_{MIN}}{dt}=BB1+BB2+BB3$$

(50)

BB1 represents sedimentation, BB2 represents decomposition of organic nitrogen in bottom sediment, and BB3 represents adsorption/desorption onto mud particles.

BB1 is expressed using the deposition rate of inorganic nitrogen, S_IN, as follows:

$$BB1=S_{IN}\frac{\partial C_{MIN}}{\partial z}$$

(51)

To reproduce the difference in degradation rates between anaerobic and aerobic conditions, BB2 is expressed using the degradation rate of organic nitrogen in the bottom mud under aerobic conditions, $D_{SN1}$, and the degradation rate of organic nitrogen in the bottom mud under anaerobic conditions, $D_{SN2}$, as follows,

$$\begin{gathered} BB2=\frac{1}{1+\lambda \alpha }{D}_{SN1}·{\theta }_N^{\left(T-20\right)}·C_{MON}·\frac{C_{DO}}{D{O}_{SON}+C_{DO}}\mathrm{aerobic}\mathrm{conditions} \\ \frac{1}{1+\lambda \alpha }{D}_{SN2}·{\theta }_N^{\left(T-20\right)}·C_{MON}·\frac{D{O}_{SON}}{D{O}_{SON}+C_{DO}}\mathrm{anaerobic}\mathrm{condition} \end{gathered}$$

(52)

where λ is the pore water volume—dry mud ratio, expressed using the porosity ϕ as follows:

$$\lambda =\frac{\varphi }{2.65\left(1-\varphi \right)}·{10}^{-3}$$

(53)

In addition, α is the distribution ratio and is expressed as follows:

$$\mathsf{\alpha }=12.4·{0.745}^{DO}·{\theta }_N^{\left(T-20\right)}$$

(54)

Part of the inorganic nitrogen, phosphorus, and organic matter produced by the decomposition of organic matter in the bottom mud is dissolved in the pore water, and the rest is adsorbed on the mud particles. BB3 adsorption on mud particles is expressed by the following equation.

$$BB3=\frac{\frac{C_{MIN}}{\alpha }-C_{EXN}}{1+\frac{1}{\lambda \alpha }}·\frac{1}{T_{AD}}$$

(55)

where $T_{AD}$ is the time required to reach adsorption equilibrium.

2.

Inorganic phosphorus in bottom sediment;

Inorganic phosphorus can be roughly divided into insoluble (Fe-P, Al-P, Ca-P) and dissolved (${\mathrm{PO}}_{4}P$). Insoluble inorganic phosphorus is mainly found in sediments and suspended solids, while soluble reactive phosphorus is found in water. In the sediment model, inorganic phosphorus in the sediment ($C_{MIP}$) refers to insoluble inorganic phosphorus, while inorganic phosphorus in the pore water ($C_{WIP}$) refers to soluble reactive phosphorus.

The time variation of inorganic phosphorus in bottom sediment is expressed as follows:

$$\frac{d{C}_{MIP}}{dt}=BB1+BB4+BB5-BB6$$

(56)

BB1 represents sedimentation, BB4 represents decomposition of organic phosphorus in bottom mud, BB5 represents adsorption and desorption to mud particles, and BB6 represents desorption to pore water.

Deposition BB1 is expressed using the deposition rate of inorganic phosphorus, $S_{IP}$, as follows:

$$BB1=S_{IP}\frac{\partial C_{MIP}}{\partial z}$$

(57)

Degradation of organic phosphorus in the bottom mud BB4 is expressed using the degradation rate of organic phosphorus in the bottom mud, $D_{SP}$, as follows:

$$BB4=\frac{1}{1+\lambda \alpha }{D}_{SP}·{\theta }_P^{\left(T-20\right)}·C_{MOP}$$

(58)

Adsorption BB5 on mud particles is expressed as follows:

$$BB5=\frac{\frac{C_{MIP}}{\alpha }-C_{EXP}}{1+\frac{1}{\lambda \alpha }}·\frac{1}{T_{AD}}$$

(59)

Insoluble inorganic phosphorus associated with metal salts (mainly iron) is stable under aerobic conditions, but under anaerobic conditions, iron is quickly reduced and dissolved into the water as soluble reactive phosphorus. Therefore, BB6 is adsorbed into pore water under aerobic conditions when the dissolved oxygen concentration is above $D{O}_{CR}$, and under anaerobic conditions when it is below $D{O}_{CR}$, and under anaerobic conditions when the dissolved oxygen concentration is below $D{O}_{CR}$. When the dissolved oxygen concentration is below $D{O}_{CR}$, and under anaerobic condition, it can be expressed as the following equation:

$$BB6=D{E}_O·C_{EXP}·{\theta }_{DEO}^{\left(T-20\right)} \left(C_{DO}\ge D{O}_{CR} \right)D{E}_A·C_{EXP}·{\theta }_{DEA}^{\left(T-20\right)} (C_{DO}<D{O}_{CR} )$$

(60)

where $D{E}_O$ is the rate constant of desorption under aerobic conditions and $D{E}_A$ is the rate constant of desorption under anaerobic conditions.

3.

Organic nitrogen in bottom sediment;

The time variation of organic nitrogen in bottom sediment is expressed as follows:

$$\frac{d{C}_{MON}}{dt}=BB1-BB2$$

(61)

BB1 is sedimentation and BB2 is decomposition of organic nitrogen in the bottom sediment.

The deposition of BB1 is expressed by the following equation using the deposition rate of organic-form nitrogen, $S_{ON}$.

$$BB1=S_{ON}\frac{\partial C_{MON}}{\partial z}$$

(62)

The decomposition of organic-form nitrogen in the bottom mud of BB2 is expressed by the following equation:

$$BB2=D_{SN1}·{\theta }_N^{\left(T-20\right)}·C_{MON}·\frac{C_{DO}}{D{O}_{SON}+C_{DO}}\mathrm{aerobic}\mathrm{conditions}{D}_{SN2}·{\theta }_N^{\left(T-20\right)}·C_{MON}·\frac{D{O}_{SON}}{D{O}_{SON}+C_{DO}}\mathrm{anaerobic}\mathrm{condition}$$

(63)

4.

Organic phosphorus in bottom sediment;

The time variation of organic phosphorus in bottom sediment is expressed as follows:

$$\frac{d{C}_{MOP}}{dt}=BB1-BB4$$

(64)

BB1 is sedimentation and BB4 is decomposition of organic phosphorus in the bottom sediment.

The sedimentation of BB1 is expressed by the following equation using the sedimentation rate of organic-form phosphorus, $S_{OP}$:

$$BB1=S_{OP}\frac{\partial C_{MOP}}{\partial z}$$

(65)

The degradation of organic phosphorus in BB4 bottom mud is expressed by the following equation:

$$BB4=D_{SP}·{\theta }_P^{\left(T-20\right)}·C_{MOP}$$

(66)

5.

Organic matter in bottom sediment;

The time variation of organic matter in the bottom sediment is expressed as follows:

$$\frac{d{C}_{MOC}}{dt}=BB1-BB7$$

(67)

BB1 represents sedimentation and BB7 represents dissolution of organic matter in the bottom sediment.

The sedimentation of BB1 is represented by the following equation using the sedimentation rate of organic matter, $S_{OC}$:

$$BB1=S_{OC}\frac{\partial C_{MOC}}{\partial z}$$

(68)

Dissolution BB7 of organic matter in the bottom mud is expressed using the dissolution rate of organic matter D_SC by the following equation:

$$BB7=D_{SC}·{\theta }_C^{\left(T-20\right)}·C_{MOC}$$

(69)

6.

Ammonia nitrogen in pore water;

The time variation of ammonia nitrogen in pore water is expressed as follows:

$$\frac{d{C}_{WN{H}_{4}N}}{dt}=BB2-BB3+BB8-BB9$$

(70)

BB2 represents decomposition of organic nitrogen in bottom mud, BB3 represents adsorption/desorption on mud particles, BB8 represents diffusion, and BB9 represents nitrification of ammonia nitrogen.

The decomposition of organic nitrogen in BB2 bottom mud is expressed by the following equation:

$$\begin{gathered} BB2=\frac{\alpha }{1+\lambda \alpha }{D}_{SN}·{\theta }_P^{\left(T-20\right)}·C_{MON}·\frac{C_{DO}}{D{O}_{SON}+C_{DO}}\mathrm{aerobic}\mathrm{conditions} \\ \frac{\alpha }{1+\lambda \alpha }{D}_{SN}·{\theta }_P^{\left(T-20\right)}·C_{MON}·\frac{D{O}_{SON}}{D{O}_{SON}+C_{DO}}\mathrm{anaerobic}\mathrm{conditions} \end{gathered}$$

(71)

The adsorption/desorption on BB3 mud particles is expressed by the following equation:

$$BB3=\frac{1}{\lambda }·\frac{\frac{C_{MIN}}{\alpha }-C_{EXN}}{1+\frac{1}{\lambda \alpha }}·\frac{1}{T_{AD}}$$

(72)

BB8 diffusion is expressed by the following equation using the diffusion coefficient DC in pore water,

$$BB8=\frac{\partial }{\partial z}\left(\lambda ·DC·\frac{\partial C_{WN{H}_{4}N}}{\partial z}\right)$$

(73)

BB9 nitrification of ammonia nitrogen is expressed using the nitrification rate constant for ammonia nitrogen, NR, and the half-saturation constant for DO concentration, $K_{NR}$, as follows:

$$\mathrm{BB}9=\mathrm{NR}·\frac{C_{WDO}}{\left(D{O}_{NR}+C_{WDO}\right)}·C_{WN{H}_{4}N}$$

(74)

7.

Nitrite nitrogen in pore water;

The time variation of nitrite nitrogen in pore water is expressed as follows:

$$\frac{d{C}_{WN{O}_{2}N}}{dt}=BB8+BB9-BB10$$

(75)

BB8 represents diffusion, BB9 represents nitrification of ammonia nitrogen, and BB10 represents nitrite nitrogen.

BB8 diffusion is represented by the following equation using the diffusion coefficient DC in pore water.

$$BB8=DC·\lambda ·\frac{\partial }{\partial z}\left(\frac{\partial C_{WN{O}_{2}N}}{\partial z}\right)$$

(76)

The nitrification of BB9 ammonia nitrogen is expressed by the following equation:

$$\mathrm{BB}9=\mathrm{NR}·C_{WNH4N}$$

(77)

BB10 nitrification of nitrite nitrogen is expressed by the following equation using the nitrification rate constant NR2 for nitrite nitrogen and the half-saturation constant $K_{NR2}$ for DO concentration,

$$\mathrm{BB}10=\mathrm{NR}2·\frac{C_{WDO}}{\left(D{O}_{NR2}+C_{WDO}\right)}·C_{WN{O}_{2}N}$$

(78)

8.

Nitrate nitrogen in pore water;

The time variation of nitrate nitrogen in pore water is expressed as follows:

$$\frac{d{C}_{WN{O}_{3}N}}{dt}=BB8+BB10-BB11$$

(79)

BB8 represents diffusion, BB10 nitrite nitrogen nitrification, and BB11 nitrate nitrogen denitrification.

BB8 diffusion is represented by the following equation using the diffusion coefficient DC in pore water,

$$BB8=DC·\lambda ·\frac{\partial }{\partial z}\left(\frac{\partial C_{WN{O}_{3}N}}{\partial z}\right)$$

(80)

The nitrification of BB10 nitrite nitrogen is expressed as follows:

$$\mathrm{BB}10=\mathrm{NR}2·C_{WN{O}_{2}N}$$

(81)

BB11 denitrification of nitrate nitrogen is expressed by using the denitrification rate constant DNR of nitrate nitrogen and the half-saturation constant $K_{DNR}$ against DO concentration as follows:

$$\mathrm{BB}11=\mathrm{DNR}·\frac{C_{WDO}}{\left(D{O}_{DNR}+C_{WDO}\right)}·C_{WN{O}_{3}N}$$

(82)

9.

Inorganic phosphorus in pore water;

The time variation of inorganic phosphorus in pore water is expressed as follows:

$$\frac{d{C}_{WIP}}{dt}=BB4-BB5+BB6-BB8$$

(83)

BB4 represents the decomposition of organic phosphorus in bottom mud, BB5 represents adsorption/desorption on mud particles, BB6 represents desorption from bottom mud, and BB8 represents diffusion.

The decomposition of organic phosphorus in BB4 bottom mud is expressed by the following equation:

$$BB4=\frac{\alpha }{1+\lambda \alpha }{D}_{SP}·{\theta }_P^{\left(T-20\right)}·C_{MOP}$$

(84)

BB5 adsorption/desorption on mud particles is expressed as follows:

$$BB5=\frac{1}{\lambda }·\frac{\frac{C_{MIP}}{\alpha }-C_{EXP}}{1+\frac{1}{\lambda \alpha }}·\frac{1}{T_{AD}}$$

(85)

Desorption from BB6 bottom mud is expressed by the following equation:

$$\begin{gathered} BB6=D{E}_O·C_{EXP}·{\theta }_{DEO}^{\left(T-20\right)} \left(C_{DO}\ge D{O}_{CR} \right) \\ D{E}_A·C_{EXP}·{\theta }_{DEA}^{\left(T-20\right)} (C_{DO}<D{O}_{CR} ) \end{gathered}$$

(86)

BB8 diffusion is expressed as follows:

$$BB8=DC·\lambda ·\frac{\partial }{\partial z}\left(\frac{\partial C_{WIP}}{\partial z}\right)$$

(87)

10.

Dissolved organic matter in pore water;

The time variation of dissolved organic matter in pore water is expressed as follows:

$$\frac{d{C}_{WOC}}{dt}=BB7+BB8-BB12$$

(88)

BB7 represents dissolution of organic matter in bottom sediment, BB8 represents diffusion, and BB12 represents dissolution of organic matter in pore water.

Dissolution of organic matter in BB7 bottom mud is expressed by the following equation:

$$BB7=\frac{1}{1+\lambda }{D}_{SC}·{\theta }_C^{\left(T-20\right)}·C_{MOC}$$

(89)

BB8 diffusion is expressed as follows:

$$BB8=DC·\lambda ·\frac{\partial }{\partial z}\left(\frac{\partial C_{WOC}}{\partial z}\right)$$

(90)

The dissolution of organic matter in BB12 pore water is expressed by the following equation using the dissolution rate of organic matter in pore water $D{E}_{SC}$:

$$BB12=D{E}_{SC}·{\theta }_C^{\left(T-20\right)}·C_{WOC}$$

(91)

11.

Dissolved oxygen in pore water;

The time variation of dissolved oxygen in pore water is expressed as follows:

$$\frac{d{C}_{WDO}}{dt}=BB8-BB13-BB14-BB15$$

(92)

BB8 represents diffusion, BB13 represents oxygen consumption due to nitrification of ammonia nitrogen, BB14 represents oxygen consumption due to nitrification of nitrite nitrogen, and BB15 represents oxygen consumption due to decomposition of organic matter in pore water.

BB8 diffusion is expressed by the following equation:

$$BB8=DC·\lambda ·\frac{\partial }{\partial z}\left(\frac{\partial C_{WDO}}{\partial z}\right)$$

(93)

BB13, BB14, and BB15 are expressed as follows using the oxygen consumption $N{H}_4{N}_{DO}$ from nitrification of ammonia nitrogen, $N{O}_2{N}_{DO}$ from nitrite nitrogen, and $C_{DO}$ from decomposition of organic matter in pore water.

$$BB13=\mathrm{NR}·\frac{C_{WDO}}{\left(D{O}_{NR}+C_{WDO}\right)}·C_{WN{H}_{4}N}·N{H}_4{N}_{DO}$$

(94)

$$BB14=\mathrm{NR}2·\frac{C_{WDO}}{\left(D{O}_{NR2}+C_{WDO}\right)}{C}_{WN{O}_{2}N}·N{O}_2{N}_{DO}$$

(95)

$$BB15=-D{E}_{SC}·{\theta }_N^{\left(T-20\right)}·C_{WOC}·C_{DO}$$

(96)

3. Results and Discussion

Figure 3a,b show seasonal and inter-annual change in vertical distribution of observed and simulated water temperature at Imazu-oki from April 2007 to March 2014. Regarding seasonal changes in the vertical distribution of water temperature, the observed results and calculated values generally agreed. As for the inter-annual change in water temperature, although there was a slight difference between the observed and calculated values for water temperature in the thermocline, the observed and calculated values for the surface and bottom layers generally agreed. In 2008, 2010, and 2013, it was observed that air temperatures tended to be higher than usual in the summer, and that water temperature stratification became stronger. Although it was a cool summer in 2009, a similar trend was observed when comparing calculated values and observed values. Figure 4 shows a comparison of seasonal change in vertical distribution of observed and simulated (a) NH₄N, (b) NO₂N, (c) NO₃N, (d) PO₄P, (e) T-N, (f) T-P, and (g) DO at Imazu-oki from April 2007 March in 2008. Figure 5 shows a comparison of inter-annual change in vertical distribution of observed and simulated phytoplankton (a) Chlorophyll a, (b) b, (c) c, (d) NH₄N, (e) NO₂N, (f) NO₃N, (g) PO₄P, (h) T-N, (i) T-P, and (j) DO in the surface layer (0 m depth), thermocline (20 m depth), and deep layer (90 m depth) at Imazu-oki from April 2007 to March 2014.

In order to verify the accuracy of this ecosystem model, the numerical calculation results were compared with the measured data on the present amount of nutrients and DO that seem to have the greatest impact on the water quality environment in the inland lake. The numerical results generally represent seasonal increases and decreases, and seem to be in agreement with the quantitative results.

3.1. Seasonal and Inter-Annual Change in Phytoplankton

The main factors for phytoplankton proliferation are the nutrients nitrogen and phosphorus in the water [39]. The community structure of phytoplankton in the central part of Imazu-oki shows an increasing trend of cyanobacteria and flagellates (Figure 5b,c), and the proportion of small phytoplankton increases. Cyanobacteria have chlorophyll a and b and chlorophyll b in abundance. Diatoms have chlorophyll a, c, and especially chlorophyll c in abundance. Green algae have chlorophyll a and b, with especially high chlorophyll a. Based on the analysis of lake sediments and the monitoring data conducted by Shiga Prefecture, the phytoplankton biomass in Lake Biwa has tended to increase year by year since the 1960s, but began to decrease in the 1990s. This decrease was mainly due to a decrease in green algae, while cyanobacteria biomass increased or remained high after 1990. In addition to these changes in biomass, community structure has also been shown to have changed over the past several decades [40]. The wind speed and air temperatures in winter, in addition to nutrient loads, have affected phytoplankton abundance and species composition [41]. The pattern of seasonal change is characterized by an increase in yellow flagellates (chlorophyll c) in the spring and a large increase in green algae (chlorophyll a) in the fall.

3.2. Seasonal and Inter-Annual Change in Inorganic Nitrogen

Inorganic nitrogen exists in the form of ammonia nitrogen (NH₄⁺), nitrite nitrogen (NO₂⁻), and nitrate nitrogen (NO₃⁻). Nitrate ions are abundant in rain and river inflows, and after rainfall, nitrate ion concentrations in lake water can be high [42]. This inorganic nitrogen is taken up by phytoplankton as nutrients and assimilated into their somatic cells. Organic nitrogen, on the other hand, is present in the lake as amino acids and urea, which are decomposed by bacteria and animals through predation on organic matter contained in the lake plankton, the carcasses of other organisms, and fallen leaves that come in from outside the lake [43]. In many lakes, an increase in ammonium ions is often observed in the bottom layer during summer due to the decomposition of organic matter in the bottom sediment [44]. In addition to this, ammonia and urea are emitted directly from organisms. These are then taken up again as a source of nitrogen for phytoplankton. The increase in ammonium ions in the bottom sediment during the summer is reproduced in the simulation results (Figure 4a). Among phytoplankton, cyanobacteria directly utilize molecular nitrogen gas dissolved in water to fix nitrogen when nutrient nitrogen compounds are insufficient [45]. This is consistent with the simulation results (Figure 5a), in which chlorophyll growth is limited in the summer months when inorganic nitrogen species tend to decrease in the surface layer.

In addition, the form of nitrogen in lakes and marshes changes as a result of changes in the redox state of the water. Under oxidative conditions, ammonium ions are oxidized by bacteria, and their forms change to nitrite ions and then to nitrate ions, as shown in Figure 4a–c. This is referred to as nitrification by both ammonia-oxidizing bacteria and nitrite-oxidizing bacteria active under aerobic conditions. On the other hand, under reducing conditions, when the dissolved oxygen in the water almost is absent, nitrate ions are reduced by bacteria to nitrogen gas via nitrite ions and nitrous oxide (N₂O). This process is referred to as denitrification. Denitrification begins when the water just above the lake bottom is deoxidized. In general, the reaction mechanism of nitrification and denitrification is used as a powerful treatment method to decompose and remove nitrogen from sewage, such as domestic wastewater. Thus, nitrogen circulates in lakes in various chemical forms [46].

In the bottom layer, nitrogen nutrients were higher in the summer months of June and September, with a similar trend for phosphate phosphorus (Figure 4a–c and Figure 5d–g). Regarding the observation results of inorganic nitrogen, nitrate nitrogen is detected without ammonia, despite the low oxygen concentration in the bottom layer. This may indicate a problem in the sample handling or analysis method, since ammonia should have been detected [47,48]. The concentration of substances in the bottom sediment before the dissolution was high, and the difference in concentration between the bottom sediment and the water directly above it was large, suggesting that the dissolution rate increased. This might be because the decomposition of organic nitrogen to produce ammonia nitrogen and the leaching of ammonia nitrogen from the bottom sediment into the pore water proceeded smoothly [44].

For nitrate nitrogen (NO₂N), the trend of low in the warm season and high in the cold season was observed from the surface layer to the thermocline, and the trend of high in the warm season and low in the cold season was observed in the deep layer, and this trend was also similar to the calculated values. The following are possible reasons for this. Ammonia- and nitrate-form nitrogen are incorporated into phytoplankton cells during their proliferation. After the death of phytoplankton, nitrogen is released into the water in the ammonia form through decomposition and predation, and then reverts to the nitrate form through nitrification. During the active growth period of phytoplankton, the uptake rate of nitrate ions exceeds the regression rate, and the concentration of nitrate nitrogen in the surface water layer during the stratification period decreases as phytoplankton that have taken up nitrate ions settle [49].

3.3. Seasonal and Inter-Annual Change in Inorganic Phosphorus

Soluble reactive phosphorus that flows into lakes becomes suspended phosphorus. The suspended phosphorus is also relatively light in specific gravity, but it settles and accumulates in the bottom sediment through processes such as decomposition in the water. Phosphorus from the bottom sediment is supplied to the water by the current. The particles that are swept up adsorb and desorb phosphate-form phosphorus in the water, and solid–liquid phosphorus exchange occurs in the water as well [50].

Phosphate-form phosphorus in the solid phase adsorbed on suspended solids and bottom sediments is thought to be released to the liquid phase (water and pore water) as environmental conditions such as pH and oxidation-reduction potential change, and the phosphorus in the suspended solids and bottom sediments is a source for phosphate-form phosphorus in the water. Phosphate phosphorus is the most readily available form of phosphorus for living organisms, and an internal phosphorus cycle is established in lakes in which released phosphate phosphorus is reused by phytoplankton and converted to bioavailable phosphorus. These processes are thought to involve a complex combination of physical processes, such as the supply of suspended sediment derived from rivers and bottom sediments, chemical processes associated with changes in pH, oxidation-reduction potential, etc., and biological processes, such as the synthesis of compounds [51]. One key to these biochemical reactions lies in the various compounds contained in bottom sediments and suspended solids. For example, a long-standing theory on inorganic phosphorus states that when reduced near the bottom of a lake or in bottom sediment, Fe³⁺ is reduced to Fe²⁺, and at the same time phosphorus adsorbed on Fe³⁺ is released into the pore water and water [52]. On the other hand, the reaction of phosphate-form phosphorus in water, especially in the water column of shallow lakes, is generally under aerobic conditions, so the release of phosphorus into the water requires aerobic phosphorus leaching conditions. Even in such cases, sediment sorption complex, pH, and other conditions present in the suspension have a significant influence on the phosphorus exchange between water and suspension [52]. On the other hand, inorganic phosphorus in the adsorbed form on Ca²⁺ is a stable form that is unlikely to be a source of phosphate-form phosphorus in water.

The reason for the higher phosphate-form phosphorus concentrations in the surface layer in 2009, 2010, and 2013 is that, as shown in Figure 5g, in addition to the leaching of phosphate-form phosphorus from the bottom mud in previous years, the cold surface water of the winter season mixed vertically to the bottom layer in the overturning, carrying the phosphate-form phosphorus layer to the surface layer. Phosphate-form phosphorus concentrations in the bottom layer in the fall were also higher in 2009 and 2013 than in the previous year because water stagnated in the bottom layer. Although the increases in T-N and T-P concentrations in the bottom layer over the seasons were not at problematic levels, phosphate-form phosphorus concentrations increased in the bottom layer in the fall of 2009, and when vertical mixing after overturning was promoted, the phosphate-form phosphorus concentration in the surface layer of North Lake was equal to the maximum value ever recorded (Figure 5g).

Although both the measured and calculated values varied widely, the seasonal changes in phytoplankton were generally higher in the warm season and lower in the cold season, and the seasonal changes in phosphate phosphorus were lower in the warm season and higher in the cold season, which were almost similar trends. The annual changes in T-P were also similar for the measured and calculated values.

3.4. Seasonal and Inter-Annual Change in Total Nitrogen and Phosphorus

T-N showed the same seasonal variation as inorganic nitrogen (Figure 4e and Figure 5h), but the variation was less than that of inorganic nitrogen, which decreased during the summer and increased from winter to spring. The concentration of T-N was higher in the bottom layer than in the surface layer during the summer. The reason why T-N concentrations were a little low from January to February in the bottom layer (Figure 4e and Figure 5h) is that in winter, vertical mixing supplies oxygen to the bottom layer, forming an oxidized layer that prevents leaching [53]. On the other hand, T-P increased in summer and decreased from winter to spring, opposite to the seasonal change in T-N as shown in Figure 5i. Comparing the concentrations of T-P in the surface and bottom layers, the concentration of T-P in the bottom layer was higher than in the surface layer in the observation unlike the simulation (Figure 4f), which was also due to the sedimentation of nutrients from the surface layer and the leaching of nutrients from the bottom mud.

3.5. Seasonal and Inter-Annual Change in Dissolved Oxygen

As the air temperature rose, the stratification was formed (Figure 3a,b), and the dissolved oxygen concentration, which had been supplied by the atmosphere and diffused from the surface to the bottom layer, began to decrease. When the thermocline with a large water temperature gradient was formed, the bottom layer was no longer supplied with oxygen, but only decomposed organic matter by bacteria, resulting in oxygen consumption (Figure 4g). In the surface layer, oxygen was supplied through photosynthesis by plant cells within the range of light penetration, and oxygen supply from the atmosphere was the main factor, so oxygen depletion was less likely to occur than in the bottom layer. Oxygen depletion continued until around September, when a strong water temperature stratification was formed. From fall through winter, the water temperature stratification was destroyed as the air temperature dropped, oxygen was supplied from the surface layer, and the oxygen concentration recovered, and during the full layer circulation period in winter, the oxygen concentration in the bottom layer approached that of the surface layer. This is consistent with the results of Ito et al. [54].

Figure 5j shows the change in dissolved oxygen concentration over time from April 2007 to March 2014. The saturated dissolved oxygen concentration is dependent on the water temperature because the surface water temperature changed with the air temperature in each year [55], and the model results are consistent with this (Figure 3b). The stirring of the surface layer by the wind affected the DO concentration fluctuations in the thermocline and the surface layer. The oxygen concentration in the bottom layer decreased in 2008, 2010, and 2013 (Figure 5j). This is because when the water temperature in the surface layer rose above normal due to the change in air temperature in each year (Figure 3b), the stratification became stronger, and it was difficult for oxygen to be supplied from the surface layer. Oxygen consumption was also high in years with high sedimentation of organic matter. The comparison of observed and calculated values shows that these trends are generally consistent.

4. Conclusions

The vertical flow field was resolved using a vertical one-dimensional diffusion model, and a novel water quality model was developed by incorporating a bottom sediment model that accounts for anaerobic and aerobic processes within bottom mud sediments into a suspended ecosystem model that considers three types of phytoplankton. This comprehensive model includes parameters such as organic and inorganic nitrogen and phosphorus in bottom sediments, organic and inorganic nitrogen and phosphorus in pore water, adsorbed nitrogen and phosphorus, and dissolved oxygen. To predict these spatio-temporal changes, processes such as sedimentation, diffusion, decomposition, adsorption/desorption, nitrification, and denitrification were considered.

The behavior of inorganic nutrients—specifically ammonia nitrogen, nitrate nitrogen, nitrite nitrogen, and phosphate phosphorus—along with seasonal variations in dissolved oxygen concentrations, was compared with observational data. Notably, the model reproduced the seasonal and inter-annual variations in dissolved oxygen (DO) concentrations in the surface layer, the thermocline, and the bottom layer. The seasonal variation of NH₄-N in the bottom pore water in the study area corresponded to fluctuations in bottom sediment total organic carbon (TOC), indicating that the supply of organic matter to the lake bottom and water temperature significantly influenced these variations. NO₃-N in the bottom pore water was depleted during the summer/fall (July–October) and winter (January) periods, primarily due to the depletion of NH₄-N, which serves as the substrate for denitrification and nitrification, respectively. Seasonal variations in PO₄-P in bottom mud pore water were closely related to changes in the redox environment of the bottom mud.

The comparison between the model outcomes and observed data demonstrated that the dissolved oxygen concentration in the bottom layer of Lake Biwa was well reproduced. Furthermore, the model effectively accounted for the phenomena resulting from the decomposition of organic-form nitrogen in the sediment over time under aerobic conditions and the effects of changes in inorganic-form nitrogen and phosphate-form phosphorus in the pore water on the bottom layer. Although some aspects could not be adequately reproduced, necessitating further studies, the developed model accurately represented the behavior of inorganic nutrients in the bottom layer under both aerobic and anaerobic conditions and reproduced oxygen variability. The model provided significant insights into the interactions of nutrients with water, pore water, and bottom sediments. Given that in many enclosed water bodies, newly settled sediments form new bottom layers, promoting anoxia, this model will be instrumental in analyzing the effects of anoxia as the vertical flow field near the lake bottom changes with climate change.

By integrating this model with a water quality model based on a one-dimensional diffusion framework and a suspended ecosystem model that incorporates the bottom sediment system, this model can be applied not only to Lake Biwa but also to other lakes and marshes. Additionally, by considering saltwater dynamics, it can be used for comparative studies on water quality and sediment depletion issues in other enclosed water bodies such as oceans, lakes, and reservoirs.