Numerical Simulation of Phosphorus Release with Sediment Suspension under Hydrodynamic Condition in Mochou Lake , China

Phosphorus is a major cause of lake eutrophication. Understanding the characteristics regarding the release of phosphorus from sediments under hydrodynamic conditions is critical for the regulation of lake water quality. In this work, the effects of sediment suspension on the release characteristics of phosphorus from sediment were investigated under different hydrodynamic conditions. The experimental results showed that in the experimental process, the phosphorus was at first released quickly into the overlying water but then slowed down. Furthermore, the process of dissolved phosphorus (DP) release under hydrodynamic conditions with and without sediment suspension was simulated using a lattice Boltzmann method. The simulation showed satisfying results.


Introduction
Sewage, industrial discharge, and agricultural runoff contribute to most sources of phosphorus released into shallow lakes, and most of the phosphorus is accumulated in lake sediments, which come back to overlying water under certain conditions [1,2].Investigations on large and shallow lakes [3,4] indicated that wind waves are critical to sediment resuspension processes.Dynamic behaviors of suspended sediments and wind-wave effects have been studied in past decades and has shown that the suspended sediment has a strong effect on the release of phosphorus [5].Phosphorus is a major factor that leads to the eutrophication of lakes as it can be absorbed by the vegetation of a lake, and studies have shown that strong wave disturbance can double the phosphorus concentration in the overlying water of lakes [6,7].Therefore, understanding the release characteristics of the phosphorus from sediments under hydrodynamic conditions is critical for the regulation of lake water quality.
Two methods were proposed to investigate the influence of hydrodynamic forces on the release of phosphorus from sediment.One is to measure the concentration of phosphorus in natural lakes under natural hydrodynamic conditions, and then to analyze the effect of the hydrodynamic force on the phosphorus release from the sediment.The other is to take sediment samples from lakes and then simulate the hydrodynamic conditions as natural conditions using laboratory instruments [8].In terms of the research on natural lakes, Nilsson et al. measured the phosphorus concentration in the overlying water of sediment under different hydrodynamic conditions in recirculating channel and found that the change of hydrodynamic force could significantly affect the phosphorus concentration in overlying water [9].Previous research has shown that sediment concentration has a significant correlation with nutrients [10].Wang et al. [11] found a positive correlation between total phosphorus (TP) and total suspended sediment (TSS) concentrations.Moreover, numerical methods have been developed to simulate the phosphorus concentrations in lakes and rivers [12].Some of these methods consider the processes of adsorption-desorption and the release from bed sediment [13].
In recent years, the lattice Boltzmann method (LBM) has been widely applied to the study of diffusion.Wolf-Gladrow [14] gave the basic equations and the parameter selection method of LBM for solving diffusion problems.Jiaung et al. [15] used the LBM to simulate the process of heat conduction, where the thermal diffusivity was controlled by changing the relaxation time.Chen showed that the LBM is also applicable to simulate the diffusion problem in porous media with uniform and heterogeneous porosities, such as lake sediment [16].
In this study, laboratory experiments were conducted to investigate the release characteristics of phosphorus from sediment in shallow lakes.Experiments with and without sediment suspension were designed.The main objectives were as follows: (1) to investigate the release characteristics of the TP and dissolved phosphorus (DP) from the sediment with and without sediment suspension, and (2) to simulate the DP release with and without sediment suspension using LBM.

Laboratory Experiments
The experiment was divided into two parts: a phosphorus release experiment with sediment suspension and a phosphorus release experiment without sediment suspension.The difference between them was the diameter of sediment, which, depending on its size, would either suspend or not under a hydrodynamic condition.

Experimental Device
The structure of the experimental device is shown in Figure 1a,b.A Plexiglas container with a diameter of 30 cm and a height of 50 cm was mounted with a variable speed motor that ran the propeller to adjust the hydrodynamic conditions in the device.

Experimental Sediment
The sediment was collected from Site 1 and Site 2 (Figure 2), respectively.Due to the previous lake slope protection project, a large amount of large diameter sediment was filled in the shore of Mochou Lake in 2015, the sediment diameter of Site 1 along the shore of Mochou Lake is larger than that of Site 2. The diameter of the sediment at Site 2 (D2) ranged from 0.15 cm to 0.23 cm, with a median diameter of 0.18 cm.The sediment was too coarse to suspend, even at the highest propeller speed (300 rad/min).The sediment diameter (D1) of Site 1 is given in Table 1.The diameter of sediment in Site 1 is small enough to have it suspend in a hydrodynamic condition.Mochou Lake has a history of 1500 years, with an average water depth of 2.5 m, a lake surface area of 0.3 km 2 , and a maximum water depth of 4 m.On 21 July 2018, we used a Petersen grab to dig a surface sediment at a depth of 0 cm to 10 cm in Mochou Lake.This lake is frequently influenced by winds, with wind-induced currents transporting dissolved matters (e.g., nutrients).There are two main parts of the sampling device: a sampling grab and a pulling rope.The sampling grab is made of a high-quality alloy material, as shown in Figure 1c.The one-time sampling sediment volume was 1-5 L, and the diameter of the sampling grab was 18 cm.

The Experimental Method
Before the experiment, the sediment samples were fully stirred, and the thickness was controlled at 10 cm.A total of 20 L of deionized water was slowly added to the Plexiglas container, and the container was placed for 24 h prior to the experiment.The speed of the propeller for the various steps was set to 100 rad/min, 200 rad/min, and 300 rad/min, respectively.The experimental parameters are shown in Table 2.
Water 2019, 11, 370 4 of 15 Three 100 g sediment samples were taken from the sediment bed before and after each experiment run.The TP concentration of the sediment was determined by Chinese environmental standard HJ632-2011: Soil Determination of Total Phosphorus by alkali fusion-Mo-Sb anti-spectrophotometry method.
For each experiment run, three 30 mL water samples were extracted from each experimental run for each hour of total experimental time.The DP and TP in the water was determined using the Mo-Sb anti-spectrophotometry method of Yuan [17], the samples for TP determination needed to be digested in advance, and the samples for DP determination needed to be filtered without being digested.In addition, the particulate phosphorus (PP) was estimated as the difference between the TP and the DP (i.e., PP = TP − DP).The TSS concentration was calculated using a drying method.The beaker with 20 mL sampling water was placed in an oven at 115 • C till all the water evaporated.The calculated equation was as follows: where W 1 is the weight after drying, W 2 is the weight of the beaker, and V is the volume of sampling water.At Site 1, the flow velocity (u z ), positioned in the water 0.5 m above the bed, was measured using a LS1206B intelligent portable velocimeter which was made by Nanjing Shunlaida Measurement and Control Equipment Co., Ltd.(Nanjing, China).The velocity measurement error was less than 1.5%, and each sampling time was 5 min.The wind speed (u w ), positioned 0.2 m above the water surface, was measured using an AR866 handheld thermosensitive anemometer made by Suzhou R.B.T Measurement and Control Technology Co., Ltd.(Suzhou, China).The wind speed measurement error was less than 1%, the measurement range was from 0.3 m/s to 30 m/s, and each sampling time was 5 min.The data was measured at 10:00-12:00 am, 2:00-4:00 pm on 22, 24 and 27 July 2018, respectively.Three groups of corresponding velocity and wind speed per hour were measured at Site 1.A statistical analysis of the field data (Figure 3) indicates that the occurrences of the u z were less than 0.25 m/s and that the u w was less than 8 m/s.
Based on the similarity principle, the rotational speeds are selected so that the Froude number in the laboratory condition (Fr m ) is equal to that in field conditions (Fr p ) [18].
Under laboratory conditions, where u m is the tangential speed in the laboratory experiment, h m is the height above the sediment surface in the laboratory experiment (herein, h m = 0.25 m), and g is the gravitational acceleration.Under the field condition Water 2019, 11, 370 where h z is the height above the sediment surface in the field (herein, h z = 0.5 m) under laboratory conditions.
From the general relation between the tangential and angular speeds, the laboratory experimental rotational speed is where P s is the propeller speed, and r is the radius of the stirring rod (herein, r = 0.045 m).This study examined the relationship between u m and u w that corresponds to the dynamic condition caused by the propeller speed to the wind speed.When measuring, the maximum wind speed can be less than 8 m/s, but sometimes the maximum wind speed in Nanjing can be larger than 8 m/s.The field wind speed varied from 4.65 m/s to 12.74 m/s, and P s varied from 0.075 to 0.225 (Table 3).As a result, the propeller speeds were set to 100, 200, and 300 rad/min in the laboratory experiment to simulate the prevailing wind speed of 4.65 m/s to 12.74 m/s in the lake.Huang et al. [18] examined the shear stress caused by lake currents, because they are considered to have similar dynamic effects to the stirring rod on the sediment resuspension, and because the dynamic simulator can produce shear stresses of water currents in the laboratory.The propeller speeds were selected to simulate the prevailing bottom flow velocities of 0 m/s to 0.08 m/s (calculated by Equation ( 6)), as in Taihu Lake.As a result, the blade stirrer operated at propeller speeds of 0, 100 rad/min, 200 rad/min, 300 rad/min, and 400 rad/min in the laboratory experiment.We used the same method to examine which propeller speeds to select to simulate the prevailing bottom flow velocities of 0 m/s to 0.04 m/s in the lake.As a result, the blade stirrer operated at the propeller speeds of 0, 100 rad/min, 200 rad/min, and 300 rad/min in this laboratory experiment.
The bottom boundary velocity is computed as where κ is the von Kármán's constant, and z 0 is the bottom boundary roughness, (we assume this based on existing literature [1,8], z 0 = 0.02 m).
where  is the tangential speed in the laboratory experiment, ℎ is the height above the sediment surface in the laboratory experiment (herein, ℎ = 0.25 m), and  is the gravitational acceleration.
Under the field condition where ℎ is the height above the sediment surface in the field (herein, ℎ = 0.5 m) under laboratory conditions.
From the general relation between the tangential and angular speeds, the laboratory experimental rotational speed is where  is the propeller speed, and  is the radius of the stirring rod (herein,  = 0.045 m).This study examined the relationship between  and  that corresponds to the dynamic condition caused by the propeller speed to the wind speed.When measuring, the maximum wind speed can be less than 8 m/s, but sometimes the maximum wind speed in Nanjing can be larger than 8 m/s.The field wind speed varied from 4.65 m/s to 12.74 m/s, and  varied from 0.075 to 0.225 (Table 3).As a result, the propeller speeds were set to 100, 200, and 300 rad/min in the laboratory experiment to simulate the prevailing wind speed of 4.65 m/s to 12.74 m/s in the lake.Huang et al. [18] examined the shear stress caused by lake currents, because they are considered to have similar dynamic effects to the stirring rod on the sediment resuspension, and because the dynamic simulator can produce shear stresses of water currents in the laboratory.The propeller speeds were selected to simulate the prevailing bottom flow velocities of 0 m/s to 0.08 m/s (calculated by Equation ( 6)), as in Taihu Lake.As a result, the blade stirrer operated at propeller speeds of 0, 100 rad/min, 200 rad/min, 300 rad/min, and 400 rad/min in the laboratory experiment.We used the same method to examine which propeller speeds to select to simulate the prevailing bottom flow velocities of 0 m/s to 0.04 m/s in the lake.As a result, the blade stirrer operated at the propeller speeds of 0, 100 rad/min, 200 rad/min, and 300 rad/min in this laboratory experiment.
The bottom boundary velocity is computed as where  is the von Kármán's constant, and  is the bottom boundary roughness, (we assume this based on existing literature [1,8],  = 0.02 m).

Simulation Method
Generally, the nutrient release in lakes undergoes the following process: the nutrient is released from sediment, the resuspended sediment is desorbed, and then it is diffused in overlying water.For sediment release, it follows the process of diffusion, adsorption, and desorption.

LBM
LBM includes two steps: migration and collision.In the migration step, the particles move in a certain direction and at a certain speed to the adjacent nodes of the grid.The migration process is expressed as: where f k is the particle distribution function in terms of a discrete particle along the direction k; f k is the value of a particle before a migration along the direction k.
Theoretically, the particle collision process is very complex and difficult to solve.Bhatnagar et al. proposed a simple BGK collision operator for a discretized LBM equation, which can be expressed as: where f eq k is the distribution of the equilibrium function along the direction k, and ω is the relaxation frequency.
Bringing Equation ( 6) into Equation ( 5), we have: In LBM, the solution region is divided into many lattices.In this paper, two dimensions with a nine-direction (D2Q4) lattice configuration were used (Figure 4).f eq k (x, z, t) = w k T(x, z, t) (10) where f k is the particle distribution function in terms of a discrete particle along the direction k, f eq k is the equilibrium distribution function along the direction k; w 1 = w 2 = w 3 = w 4 = 0.25, and T is the nutrients concentration.The LBM function for sediment in overlying water can be described as follows,

Boundary Conditions
By applying LBM, the free surface of overlying water can be defined as the thermal insulating boundary.The nutrient concentration gradient of the free surface is 0, and the container wall is defined as rebound boundary.Please refer to Mohamad [26] for an in-depth explanation of the boundary method.

Model Parameters
The parameters in this model are given in Table 5.

Phosphorus Release from Sediment
The phosphorus release in lakes includes the phosphorus release from sediment and the phosphorus diffusion in overlying water.For sediment release, it follows the process of diffusion, adsorption, and desorption.The diffusion rate can be measured by the nutrient concentration gradient in the pore water of sediment with Fick's first law, and the adsorption and desorption can be defined by a source item [19].
The phosphorus release from sediment can be expressed as: where c d is the DP concentration in water (mg/L), t is time (s), ϕ is the porosity of sediment, z is vertical axis originated (at the sediment-water interface, z = 0), D zs is the diffusion coefficient of phosphorus in sediment (m 2 /s), ρ b ∂c s ∂t is a source term for phosphorus adsorption and desorption by the sediment bed [20], c s is the quantity of phosphorus adsorption in the sediment bed (mg/kg) and ρ b is the density of sediment (kg/m 3 ).
The Lagergren first-order (LFO) equation is commonly used for describing the adsorption and desorption and for their kinetics research, which is expressed as [21]: where b 1 is the first-order rate constant: absorption efficiency of sediment bed (g/(Ls)), and c se is the sediment contamination level (mg/kg).Yuan et al. [22] assumed that the desorption amount of the sediment bed are equal to the amount added to the solution.Then, they modified the LFO equation as: where c e is the equilibrium concentration of TP in water (mg/L) [23].
The modified LFO equation only considers the constant hydrodynamic condition where c e remains unchanged in a constant condition.Therefore, in an airtight container without phosphorus input, the concentration of phosphorus in overlying water and sediment would be constant values.However, under the action of shear velocity, c e will vary with the change of the hydrodynamic condition (or the adsorption rate decreases and c e increases with the increase of shear velocity [24]).
Here the second term on the right of Equation ( 1) is modified as: D zs can be expressed as: where D zm is the molecular diffusion coefficient in water (m 2 /s) and it varies with the targeting solution; m = 3 is a constant [25].
The source term can be added in the LBM function [26], and the LBM function can be described as: where ∆t is the time step and ω p is the relaxation frequency in pore water; The relationship between the diffusion coefficient and ω p can be obtained by Chapman-Enskog expansion [26].
where p presents the dimension of the model.(p = 1 presents one dimension, p = 2 presents two dimensions, and p = 3 presents three dimensions).
In this paper, we adopted a two-dimension model, and ω can be obtained as follows:

Nutrient Diffusion in Overlying Water
In overlying water, the formulations can be simply described as a diffusion process, as the phosphorus release from TSS and the biochemical reactions are negligible.The governing equations can be expressed as [27]: where D zt is the turbulent diffusion coefficient (m 2 /s), A is the area of water-sediment interface (m 2 ), v is the kinematic viscosity of water (m 2 /s), and n = 3 is a constant [28]; S ∂c p ∂t is a source term for phosphorus adsorption and desorption by TSS, S is the TSS concentration in the overlying water, and c p is the quantity of phosphorus adsorption by TSS in the overlying water.
We assumed that the desorption amount of TSS is equal to the amount added to the solution.
Then, S ∂c p ∂t can be modified as: where b 2 is the first-order rate constant, that is, the absorption efficiency of TSS (g/(Ls)).
The LBM function and ω w can be obtained with the same method as with the sediment: ) where u * is the shear velocity, and the u * generated by the propeller speed is given by Chandler [29]:

.4. Nutrient Desorption from Resuspended Sediment
In this research, the sediment concentration in the overlying water was assumed to be uniform.Here, the dynamic release amount caused by sediment suspension is much larger than the constant release amount caused by diffusion [30], while the vertical velocity component is relatively smaller than the convection velocity component [31].Therefore, only the effects of sediment re-suspension and sedimentation are considered, and Equation ( 4) could be simplified into: The settling velocity of sediment particles was obtained by Song et al. [32]: Water 2019, 11, 370 9 of 15 where D is the particle diameter and D * is the dimensionless particle diameter.∆ = ρ s ρ − 1, where ρ s and ρ represent the density of particles and the density of the fluid, respectively.K = 9 × 10 −3 (kg/m 3 ); S s is the deposition flux; S r is the re-suspension flux; u c is the critical velocity of sediment, which can be expressed as: Van Rijn [21] related the dimensionless particle parameter to a critical sediment mobility parameter (d cr ), which is elaborated in Table 4.The LBM function for sediment in overlying water can be described as follows,

.5. Boundary Conditions
By applying LBM, the free surface of overlying water can be defined as the thermal insulating boundary.The nutrient concentration gradient of the free surface is 0, and the container wall is defined as rebound boundary.Please refer to Mohamad [26] for an in-depth explanation of the boundary method.

Model Parameters
The parameters in this model are given in Table 5.To conclude the above contents, the sediment and phosphorus transport can be calculated using LBM, as shown in Figure 5.

Experimental Results
For Run 4-6, there is very little sediment in the samples, and all of the measured TSS concentration are less than 0.02 g/L, which means that the phosphorus absorbed by TSS in the overlying water can be ignored.Therefore, only the TSS concentration of Run 1-3 are listed in Table 6.The TSS concentration is higher with the higher propeller speed.
The TP concentrations are shown in Table 7; the TP concentrations of Run 1 to Run 6 ranged from 0.059 mg/L to 0.119 mg/L, 0.063 mg/L to 0.191 mg/L, 0.065 mg/L to 0.24 mg/L, 0.031 mg/L to 0.128 mg/L, 0.04 mg/L to 0.143 mg/L, and 0.031 mg/L to 0.151 mg/L, respectively, for the 8 experimental hours.The TP concentration in the water increased quickly in the first four hours, then slowed down in the last four hours.The DP concentrations are shown in Table 8; the DP concentrations of Run 1 to Run 6 ranged from 0.048 mg/L to 0.104 mg/L, 0.034 mg/L to 0.163 mg/L, 0.034 mg/L to 0.207 mg/L, 0.037 mg/L to 0.123 mg/L, 0.039 mg/L to 0.136 mg/L, and 0.03 mg/L to 0.146 mg/L, respectively, for the 8 experimental hours.The DP concentration in the water increased quickly in the first four hours, then slowed down in the last four hours.
The variations of TP concentration in the sediment before and after experiment are shown in Figure 6.We can see that the initial TP concentration of Run 1-3 is higher than Run 4-6.After the experiment TP, decreased by 18.5%, 29.7%, 38.8%, 36.7,46.2%, and 48.9%, from Run 1 to Run 6, respectively.

Experimental Results
For Run 4-6, there is very little sediment in the samples, and all of the measured TSS concentration are less than 0.02 g/L, which means that the phosphorus absorbed by TSS in the overlying water can be ignored.Therefore, only the TSS concentration of Run 1-3 are listed in Table 6.The TSS concentration is higher with the higher propeller speed.
The TP concentrations are shown in Table 7; the TP concentrations of Run 1 to Run 6 ranged from 0.059 mg/L to 0.119 mg/L, 0.063 mg/L to 0.191 mg/L, 0.065 mg/L to 0.24 mg/L, 0.031 mg/L to 0.128 mg/L, 0.04 mg/L to 0.143 mg/L, and 0.031 mg/L to 0.151 mg/L, respectively, for the 8 experimental hours.The TP concentration in the water increased quickly in the first four hours, then slowed down in the last four hours.The DP concentrations are shown in Table 8; the DP concentrations of Run 1 to Run 6 ranged from 0.048 mg/L to 0.104 mg/L, 0.034 mg/L to 0.163 mg/L, 0.034 mg/L to 0.207 mg/L, 0.037 mg/L to 0.123 mg/L, 0.039 mg/L to 0.136 mg/L, and 0.03 mg/L to 0.146 mg/L, respectively, for the 8 experimental hours.The DP concentration in the water increased quickly in the first four hours, then slowed down in the last four hours.Froude number in the field condition f k particle distribution function in terms of a discrete particle along the direction k f k the value of a particle before migration along the direction k. f

Figure 1 .
Figure 1.The structure of experimental device and sampling device.(a) structure of the experimental device; (b) Sketch of the experimental device; (c) sampling grab.

Figure 2 .
Figure 2. Sample location and structure of the experimental device.

Figure 1 .
Figure 1.The structure of experimental device and sampling device.(a) structure of the experimental device; (b) Sketch of the experimental device; (c) sampling grab.

Figure 1 .
Figure 1.The structure of experimental device and sampling device.(a) structure of the experimental device; (b) Sketch of the experimental device; (c) sampling grab.

Figure 2 .
Figure 2. Sample location and structure of the experimental device.

Figure 2 .
Figure 2. Sample location and structure of the experimental device.

Figure 3 .
Figure 3.The relation between field flow velocity and wind speed.Figure 3. The relation between field flow velocity and wind speed.

Figure 3 .
Figure 3.The relation between field flow velocity and wind speed.Figure 3. The relation between field flow velocity and wind speed.

Figure 4 .
Figure 4. Nutrient migration process in the experimental device.

Figure 4 .
Figure 4. Nutrient migration process in the experimental device.

s density of particles ω p relaxation frequency in pore water ∆t time step v kinematic viscosity of water ρ density of fluid ϕ porosity of sediment ω relaxation frequency
eq k equilibrium distribution function along the direction k g gravitational acceleration P s propeller speed r radius of the stirring rod S TSS concentration in the overlying water S r re-suspension flux S s deposition flux u m tangential speed in the laboratory experiment u w wind speed positioned at 0.2 m above the water surface u z flow velocity positioned at 0.5 m above the bed under field conditions ρ

Table 1 .
Grain size distribution of the sediment.

Table 1 .
Grain size distribution of the sediment.

Table 1 .
Grain size distribution of the sediment.
2.1.4.Selection of Propeller Speeds Based on Field Data

Table 3 .
The contrast table of wind velocity and rotational speeds computed.
Note: u m was calculated by Equation (4); u w was calculated by the Equation in Figure3; u z was calculated by Equations (1)-(3).

Table 4 .
Variation of  with  *

Table 5 .
Parameters in the model.

Table 5 .
Parameters in the model.

Table 7 .
The TP concentration under different hydrodynamic conditions (mg/L).

Table 8 .
The DP concentration under different hydrodynamic conditions (mg/L).