Abstract
Dissolved oxygen (DO) reflects the self-purification ability of a water body and is also an important indicator for quantifying the water quality. The morphological changes in the cross sections of river channels will affect the hydraulic conditions, and the distribution of pollutants and DO may also be affected, possibly resulting in local oxygen deficits and pollution. To effectively predict the water quality, a coupled model is introduced in this study. The shallow water equation (SWE) is adopted to calculate the hydrodynamic processes, and the modified Streeter–Phelps model is further coupled with the SWE model to evaluate the reaeration. By applying this model, mass transportation and reaeration in rivers are analyzed. The influences of the sudden cross-sectional changes in the river channel on the DO distribution and the reaeration ability are identified. The results reveal that a certain degree of expansion in the river is conducive to reaeration and can also accelerate the consumption of pollutants through the water body’s self-purification. DO transport in two real terrains, including a mountain basin and plain river, is extensively investigated, and the results indicate that the morphological characteristics in the mountain basin will cause the concentration distribution to form inside dead zones, while in the plain, the distribution will form a fan-shaped downstream zone.
1. Introduction
With the rapid development of human society and the continuous acceleration of urbanization and industrialization, the pollution caused by waste discharged from construction, industry, municipal facilities, and residents’ domestic water is becoming increasingly serious, and the health and safety problems of residents’ daily water intake and water usage are also increasing. The dissolved oxygen concentration in water is an important indicator for quantifying water quality [1]. A sufficient DO content is also essential for maintaining the ecological balance in water and improving the decomposition function of organic matter. Dissolved oxygen allows water bodies to have stronger biochemical self-purification abilities. Additionally, the turbulent diffusion and transport effects in water streams can dilute the concentration of pollutants. The self-purification ability of rivers is largely through biomass conversion; hence, DO is indispensable for the decomposition of microorganisms. When the DO content in river channels is lower than the biochemical oxygen demand (BOD), the organic matter in the water body will be decomposed by microorganisms under anaerobic conditions, which will deteriorate the water quality. The distribution of DO and pollutants in rivers is always uneven due to the complex morphological characteristics of river channels and human influences, which may also cause oxygen deficits in local areas, ultimately resulting in local pollution accumulation.
There are currently many models based on different theories for analyzing water quality, such as STORM [2], SWWM [3], and HSPE [4]. However, many of these water quality models have difficulties in the theoretical treatment of pollutants and hydrodynamics and deficiencies in terms of accuracy and simulations. The Streeter–Phelps model of DO and BOD is relatively mature and has been widely used in pollution control, water environmental quality forecasting, project impact assessment, and water resource planning [5]. However, the early S-P model is considered to be relatively simple and mainly focuses on the study of oxygen balance, which is a one-dimensional steady-state model. The gas-transfer process is complicated, and research has not yet achieved a complete understanding [6]. With vast laboratory, field, and numerical studies, gas-transfer models have developed greatly, such as the two-film theory according to Lewis and Whitman [7] and the surface renewal theory raised by O’Connor and Dobbins [8]. Many empirical formulas have been proposed to estimate reaeration [9–13]. Nguyen et al. [14] unified the present experimental data and further raised a modified reaeration model, which was proven valid with vast parameters taken into consideration. However, most water quality models cannot achieve runoff computation independently, having to be coupled with hydrodynamic models for joint computation. Although hydrological models such as the SWAT model can calculate the runoff generation and concentration, they still have problems in calculating the scale and parameter calibration [15].
In this study, numerical simulations are undertaken to simulate the transport of DO in river channels, where the purpose of the modeling is to determine the relationship between the DO distribution pattern and the morphological changes in the river channels. SWE can be used to compute hydrodynamic problems and has a satisfactory accuracy. It has simulation advantages in the large-scale computational domain and calculations of runoff generation and is also suitable for handling practical problems. Lin applied SWE to investigate the transport of sediments [16]. Lin [17] and Zhang [18] applied SWE to analyze Malpasset dam-break flows. Ajayi et al. [19] and Lacasta et al. [20] applied the SWE model to simulate rainfall and overland flow and proved the accuracy of their models. In this paper, our model is applied to investigate the S-P curve, Gaussian hump, and dam break in order to verify its computing capability. After verification, we apply our model to study the sudden expansion/shrinkage of the river channels in a set of numerical experiments. Due to the computational efficiency of the average water depth, the model can be applied to study DO transport in large-scale terrains, including mountain streams in Palong Tsangpo and the Jinjiang River Basin, Chengdu.
2. Materials and Methods
The transportation of oxygen and pollutants is determined by hydrodynamics and the transport substances’ properties. The SWE program in this paper applies the assumption of an average water depth coupled with the modified S-P model to analyze the hydrodynamics and mass transport of the water body. The main outputs are the water depth and the concentration value in discretized grids.
This section will introduce the basic governing equations of the model, which consists of four parts: the hydrodynamic equations, mass transfer equation, reaeration equation, and discrete method. The shallow water model has been widely used and well developed by many scholars. The hydrodynamic method and the grid discretization method in this paper are being referred to in Lin et al.’s [17] and Zhang and Lin’s [18] works. The main governing equations have been introduced in their previous research, and they proved that this SWE method is indeed valid. Therefore, the hydrodynamic, mass transfer, and discrete algorithms are excluded in this article, and the reaeration equation will be mainly introduced.
2.1. Hydrodynamic Model
The hydrodynamic model was referred to in Zhang and Lin’s work [18], which was proven to be reliable and robust in his research. The eddies generated by turbulence are responsible for both momentum and mass transfer, resulting in a mixing effect that far exceeds that occurring at the molecular scale [21]. Thus, the depth-averaged turbulence model is applied to better describe the turbulent diffusion, has been widely used in previous research, and has a good performance in simulating the turbulent flow [22–24].
The shallow water equation is derived from the Navier–Stokes equation. It is assumed that the velocity change and hydrostatic pressure distribution in the vertical direction are negligible, and the water body is incompressible. By integrating the Navier–Stokes equation in the vertical direction, the basic governing equations of the SWE model can be formulated as follows:where is defined as the total water depth, is the height of the water surface line relative to the base level, is the height of the riverbed, (P = u × H) and (Q = × H) are the flux components per unit volume in the x direction and y direction, respectively, u and are the depth average velocity components in the x and y directions, respectively, S is the source/sink term of the water body, which is 0 when the inflow and outflow are not considered, is the acceleration due to gravity, and m is the Manning roughness coefficient.
The detailed derivation process is not the main point of this study, so the process will be omitted, and the final modeled equations will be given. The resulting equations and their derivation process can be found in Zhang and Lin’s [18] and Lin et al.’s [17] works, and the governing equations are summarized as follows:where i and j are the x and y directions of the Cartesian coordinate system, respectively, and is the dynamic viscosity. The superscript indicates time average value, is the turbulent eddy viscosity , and , , , and are model parameters, where , , , and .
The main function of these equations is to compute and further apply it to the following mass transfer equation.
2.2. Mass Transfer Equation
The three-dimensional form of the advection-diffusion equation deduced from Fick’s law iswhere is the concentration of pollutants or oxygen and , , and are the flow diffusion coefficients. On this basis, it is necessary to derive the mass transfer equation in a two-dimensional form.
For the two-dimensional model, the boundary condition can be formulated as follows (the letters with hats indicate the average depth value):
If the equation is fully expanded after being integrated along the water depth, the form will be very complicated. In this model, some of the expanded subitems are included in the source/sink items for convenience and unity. Different components are distinguished as different forms in the source/sink items. Omitting its derivation process, it is
The advection-diffusion equation for molecular diffusion and the average water depth can be expressed as
Fluids in natural river channels are mostly turbulent, where the turbulent diffusion coefficient is much larger than the molecular diffusion coefficient. In this study, the turbulence coefficient is used instead of the molecular diffusion coefficient. Therefore, the molecular diffusion equation for turbulent diffusion can be transformed aswhere is the turbulent eddy viscosity. Molecular diffusion is neglected, and only turbulent diffusion is considered. On this basis, we assumed that turbulent eddies responsible for momentum transfer are the same as those controlling the mass transfer, which means that the turbulent Schmidt number can be set as [25].
2.3. Reaeration Model
The gas-transfer process depends on the physicochemical characteristics of the substance being transferred and on the interaction between the turbulence in the atmosphere and the water body [6, 26]. In this study, the shallow water model is combined with two-film theories to simplify the reaeration, the Lamont and Scott model is also applied to better describe the turbulence reaeration [27], and Nguyen’s reaeration model based on experiments is applied to modify and supplement its molecular diffusion part [14]. For a specific grid, the change in the oxygen mass in the grid per unit time can be expressed in the following form:where is the source of oxygen, indicating the oxygen flux from air into the water body per unit time, and O is the oxygen concentration.
Currently, much research on reaeration is based on the two-film theory. According to Lewis and Whitman [7], it is assumed that there is a static thin film between the gas and liquid, and oxygen transfers through the molecular diffusion of this film. Nguyen et al. [14] further revised this model, and its form can be concluded as a unified expression:where Cs and C0 characterize the saturation concentration of oxygen in the water body and the actual concentration at the time, respectively, reflecting the degree of oxygen deficiency, and are the reaeration coefficient and the surface mass transfer coefficient, respectively, indicating the reaeration ability, and h is the water depth in this grid.
It is impossible to cover everything in actual modeling, so this research aimed to relatively simply apply a gas-transfer model as the source term. In most cases, influencing factors, such as atmospheric conditions and hydraulic slopes, can be ignored [28–32]. However, some factors truly should not be oversimplified or ignored, such as the temperature, Schmidt number (), and turbulence conditions included in this research.
Therefore, the reaeration coefficient in this paper can be concluded as the following function:where T is the temperature.
Based on the two-film theory, O’Connor and Dobbins [8] assumed that the membrane is constantly renewed and that there is a gradient of oxygen concentrations in the vertical direction. The liquid layer saturated with oxygen on the surface will be replaced by the unsaturated liquid due to turbulence, forming a new liquid film, which is the surface renewal theory. To make the formula practical, they analyzed the experimental data measured by the river and found the typical form:where is the reaeration coefficient in the reaeration model.
The formula is further simplified and imitated, and many semiempirical formulas are obtained, such as the Churchill formula, Thackston-Krenkel formula, Benett and Rathbun formula, and Owens formula [8–12].
This model refers to the work of Nguyen et al. [14] who conducted experiments to study the reaeration process. The derivation process is omitted in this paper, and the treatment of the molecular diffusion surface mass transfer coefficient also refers to Nguyen. At a specific temperature, the influencing factors of the molecular diffusion reaeration coefficient mainly include the Schmidt number and water depth [17], among which can be expressed as follows:where is the reaeration coefficient of the average water depth. By applying Nguyen’s experimental data [14], and can be represented as
This model adopts the small-eddy models of Scott et al. and Lamont and Scott [27, 33] for treatment of the turbulent part, which can be obtained as
In summary, with the two formulas added together, the reaeration coefficient of the shallow water model can be obtained as
This formula not only follows the small-eddy model but also introduces the latest stationary water reaeration model, which solves the previous dilemma where it was not applicable to stationary waters. Moreover, additional molecular diffusion also makes the theory more in line with the actual situation. Dm is determined by the Stokes-Einstein equation in the following formula:where M = 2.26, x = 25.6 cm2/mol in the water solvent, and V is the molar volume of oxygen.
2.4. Discrete Method
This model applies the finite difference method to discretize the equations, the staggering leapfrog grid to divide the computational domain, and the two-step Lax–Wendroff method to discretize the advection term of the advection-diffusion equation [34]. These methods in the model can obtain good calculation results for most cases. If numerical oscillations occur in advection-dominant situations, errors can be avoided by applying the WENO format [35].
3. Verification
The purposes of this study are to investigate the ability of coupled SWE to simulate DO transport in a shallow flow and to discuss the effect of grids on discretized computation. Thus, available analytical, numerical, and experimental results will be adopted to verify the present model. Their different results will be analyzed to check the accuracy of this model, and their consistency with theoretical or reliable numerical results in previous research will be considered as criteria of successful performance. A water quality test of the Streeter–Phelps model is set in 3.1, a 2D Gaussian hump test is set in 3.2, and a dam-break experiment containing mass transfer is set in 3.3. The hydrodynamic module of the SWE is quite mature for computing both experimental and engineering problems. Lin et al. [17] and Zhang and Lin [18] conducted many verification cases, and this model proved accurate. In addition, Zhang and Lin [18] solved numerical oscillations in the dam break by adding artificial viscosity items. Therefore, this article will omit the hydrodynamic verification part, mainly focusing on verification of the S-P model and mass transfer.
3.1. Streeter–Phelps Model Test
The Streeter–Phelps model was proposed in 1925, and so far it has still been used as a water quality modeling tool for water pollution research [5]. This model attributes the reduction of DO in the water body to the degradation of BOD. Since then, this model has been continually improved and expanded. This model mainly refers to a one-dimensional steady-state case, and the equation of the Street–Phelps model satisfies the following form:where represents the function of the BOD change, represents the function of the DO change, is the oxygen consumption coefficient, is the reaeration coefficient, and u is the flow velocity.
The analytical solution can be obtained aswhere is the initial BOD value and is the initial DO value.
In the literature on the S-P model, the coefficients are determined based on different assumptions. Our model is also a kind of S-P model but takes different coefficients. The purpose of this verification case is to compare the numerical results with those of previous scholars and to investigate the discretized computation ability of the S-P model. The same initial settings and coefficients are required to be modified temporarily and in correspondence with the verification example. In this case, the initial settings of the model proposed by Bencala and Walters [36] in the study of Uvas Creek are applied: , , and , and the coefficients adopt the settings proposed by Chapra [37]: and . The initial condition setting, and , and a computational domain of x = 2000 in the model are set up to study the changes of BOD and DO. Since the S-P model is one-dimensional, the grid sizes are and , and the grid numbers are 1000 in x and 1 in y. The parameters are listed in Table 1.
The DO and BOD changes along the way are plotted in Figure 1, and the curves clearly fit the Streeter–Phelps model’s classic form with a sag curve. This model also fits the improved Streeter–Phelps model in Chapra’s paper [37] in the discretized computation domain, showing that our model has a successful simulation of the BOD and DO concentration distribution with the theoretical results, which is also capable of simulating the BOD and DO distribution in the flume along the river channel.
3.2. Two-Dimensional Gaussian Hump Advection and Diffusion Verification
The purpose of this case is to compare the numerical results with the theoretical results and to discuss whether they are in agreement. Furthermore, the Gaussian hump is able to test the accuracy of the rectangular grid system and check the numerical stability, which is worth paying attention to in the CFD area. The case is referred to in the work of Noye and Tan who solved the theoretical solution form of the Gaussian hump in two dimensions [38], where the mass iswhere W represents the concentration of mass, t represents time, u and are the velocities in the x and y directions, respectively, D is a diffusion coefficient, and and are initial coordinates of the hump. The parameters are listed in Table 2.
The numerical computational domain of 4 m × 4 m is discretized into 320 × 320 uniform grids, with the average discretized grid sizes being and , respectively. The initial concentration is set according to the theoretical solution (20). A fixed time step is adopted, with parameters , , and D = 0.01. The initial velocities of the fluids are and . The concentration contours at t = 0 s, t = 1 s, and t = 2 s are plotted in Figure 2.
It can be clearly seen that the hump is generated at (x, y) = (0.5 m, 0.5 m), while t = 0. Due to the equal velocities and in the x and y directions, the hump moves along the line y = x and reaches (1.5 m, 1.5 m) at t = 1 s and reaches (2.5 m, 2.5 m) at t = 2 s. The transfer results are consistent with the theoretical results.
A cross-sectional diagram along the line y = x is drawn, and a comparison with the theoretical solution is shown in Figure 3, which shows that the results of the hump diffusion match the theoretical value well. This calculation example shows that the joint model consisting of SWE and the mass transfer model on two-dimensional problems in discrete grids still has accurate results. This model solves the mass transfer equation accurately and conforms to theoretical results, indicating that the model has a satisfying and reliable computing ability to solve advection and diffusion problems, such as diffusion of DO and pollutants. It thus provides a strong basis for further calculation and analysis of practical problems.
(a)
(b)
3.3. Dispersion of Transport after Dam Break
The numerical results of the abovementioned two cases are in agreement with the theoretical results, so this section will focus on the accuracy of the coupled computation. This case is a dam break containing pollutants, and it aims to verify the model’s ability to predict the transport and diffusion of oxygen or pollutants with obstacles and rapidly changing flow conditions. The dam-break case has been investigated by many authors [39–42]. Our settings are in agreement with these, and our results are compared. The model is set to a flat bottom in a 1400 m × 1400 m square domain with no-slip boundaries, without a water source/sink or inflow/outflow. The water depth on the lower side is 0.5 m, while the water depth on the higher side is 1.0 m. A dam break occurs instantaneously. With the aeration effect of the atmosphere when the water body suddenly changes ignored and the convective diffusion and turbulent effect of the water considered, the initial DO concentration is distributed in a hump around the gate and can be defined aswhere W represents the DO concentration (mg/L). The 3D view’s initial bottom and oxygen concentration setting are demonstrated in Figure 4. The water depth and concentration on either side are different, and a high concentration is set near the gate at 650 m and 600 m.
(a)
(b)
The 1400 m × 1400 m numerical computational domain is discretized into 500 × 500 uniform grids, with the average discretized grid sizes being and , respectively. The fixed time step is adopted. 3D contours at t = 200 s are plotted in Figure 5. The water depth is distributed in a symmetrical form in Figure 5(a) with two vortices behind the gates, while the concentration in Figure 5(b) is distributed asymmetrically. As the dam broke, 2 eddies with different concentrations were formed behind the gate, which can be considered to be dead zones. Dead zones are commonly due to geometrical irregularities in the flow boundary conditions. The eddies in this case have the same size and shape due to the symmetry of the dam. The dead zones cause the concentration on both sides to increase at different degrees because the initial concentration is asymmetrical. The lower side with a higher initial concentration will accumulate more DO due to the recirculation effect.
(a)
(b)
In Figure 6, the contours of the water depth and concentration at t = 200 s are compared with Li’s [43]. Li’s result has a consistent characteristic with previous scholars who have investigated this case. By comparing our result with Li’s result [43], the water depth and the concentration distribution are mainly inconsistent with his research. Since our model assumes that the turbulent eddy viscosity equals the diffusion coefficient to describe the turbulent mixing of solute, the concentration distribution will be different in some turbulent flumes. While this assumption can be considered in approximations and experiments [25], the result is still of a reference value. The form of dead zones is clear in our research. Due to momentum exchange across the interface with the main flume from the breach, flow patterns in dead zones are characterized by recirculating flows, and the flow velocity in the main stream is slowed [44]. The recirculation flows also cause the mass to accumulate to different degrees due to the asymmetric initial concentration. In this case, our model proves that it is accurate in simulating the mass transfer in turbulent flumes.
(a)
(b)
4. Results and Discussion
In river channels, cross section changes will affect flow patterns and possibly form dead zones. Additionally, the mass transport will be significantly affected by the exchange processes between the main flow in the channel and the dead zones [45]. To investigate the influence of sudden changes in cross sections on the DO distribution and to discuss how dead zones affect DO accumulation, several model tests are set up. The different degrees of shrinkage/expansion in different positions are compared. Although the complete S-P model allows us to discuss the influence of more factors or to study more indicators, such as NBOD, CBOD, and salinity, to establish a more complete water quality evaluation system, and to ensure the simplicity of model research, we only investigate DO as a preliminary indicator of water quality, which also has guiding significance for further establishment of a well-rounded water quality model.
Sudden changes in the cross section often appear in natural rivers due to morphological irregularities, such as cavities in sand or gravel beds, side arms, embayment, or obstacles [46]. Some artificial constructions, such as hydropower stations, can also transform the turbulent state of the river, form dead zones, and eventually affect the oxygen concentration distribution. The uneven distribution of DO may cause a local oxygen deficit, weaken the self-purification ability, and cause pollution. Therefore, we can rely on this well-performing model to predict the DO transmission in rivers and discuss the relevant influence of buildings or river changes.
In this section, to discuss the effects of the abrupt cross section change in detail, several numerical simulation cases of an open channel flow with transferring DO are undertaken, the subbed and wall are set as a smooth interface, the boundary condition of the sidewalls is set to no slip, and the inflow and outflow boundaries are set as free. The friction of the sidewalls is ignored, but the recirculation and backflow effect are studied. A 900 m × 60 m numerical computational domain is set, which is discretized into 180 × 30 uniform grids. At the inlet, an inflow type is set with a uniform velocity profile, and the initial hydraulic conditions are flow velocity u = 0.06 m/s, water depth h = 0.15 m, and coefficient . The average discretized grid sizes are and . In addition, the fixed time step is .
The initial concentration C of DO is 3 mg/l, the oxygen saturation concentration Cs is set as 9.09 mg/l, and the reaeration coefficient is calculated according to formula (19). The water inflow is set from the left and outflow from the right. 10 tests are set, with different degrees and locations of shrink/expansion for a comparison. The setting of the tests is summarized in Table 3, and the parameters are listed in Table 4.
The parameters are chosen to better visualize the computational results, and a lower velocity is better for observing the reaeration effect in a limited length. The initial concentration should not be too close to saturation to more clearly observe the reaeration effects and how the concentration is distributed.
Since the sudden expansion of the cross section has little effect on the upstream distribution, there is no need to change the lateral position of the expansion. The DO distribution contour (left) and the reaeration coefficient distribution contour (right) in the computational domain of the control group (without section modification) are shown in Figure 7. In the original open channel flow (test #, which is also the reference case in Figure 7), because the model also considers the influence of turbulence, the relatively strong kinetic energy near the central axis promotes the transport and redistribution of DO, while more DO accumulates near the sidewall.
(a)
(b)
4.1. Analysis of Sudden Shrink
The distribution contours of the sudden shrinkage are plotted in Figure 8. Compared to Figure 7, it can be observed that the sudden change in the cross section undoubtedly changes the original flow conditions, resulting in the redistribution of DO. The shrunken position can hinder the flume and solute from transferring along the main stream. The flume will be hindered by this obstruction, causing fluctuations in the water level and solute concentration. Before the shrunken position, the flow velocity rapidly decreases, and the water depth increases, causing the DO value to fall sharply. When the cross section shrinks by a greater amount, the influence on the distribution of DO and reaeration coefficient is relatively great.
(a)
(b)
To further discuss the influence of the shrink location, tests 4, 5, and 6, which shrunk at x = 400 m, are plotted in Figure 9. The mutation position has no significant influence on the entire distribution, but the shrunk positions are moved forward, so the entire volume of the channel is reduced, and the entire concentration rises less dramatically. However, tests 4, 5, and 6 allow us to discuss the distribution downstream, and it can be seen that the channel narrowed by 3/4 will have a more well-proportioned distribution downstream.
(a)
(b)
The DO values of the cross sections at x = 225 m, 450 m, 675 m, and 900 m are output in Figure 10. Viewing the distribution at the downstream outlet when the cross section is shrunk by 3/4, the rapid shrinkage has a redistribution impact on the flow, but the turbulence is not particularly strong. Because the reaeration at this time is relatively low, the increase in the local DO value at the sudden change is dominated by advection and diffusion rather than reaeration. According to the distribution contour of the reaeration coefficient in Figure 9(b), the reaeration coefficient distribution at the rear end of the section that is shrunk by 1/2 is better, and the reaeration capacity is strong and relatively uniform. However, due to the sharp shrinkage of the cross section, the flow velocity increases and the reaeration time for the flume transferring along the same distance is shorter. Thus, the reaeration of the section shrunk by 1/4 is better, the reaeration effect is greater, and the distribution is more uniform and less likely to form low reaeration zones.
(a)
(b)
(c)
(d)
4.2. Analysis of Sudden Expansion
The DO and the reaeration coefficient distribution is displayed in Figure 11, and the same computational domain and initial hydraulic conditions are set as in Section 4.1. Since the sudden expansion has little effect on the upstream flow, the downstream changes are mainly discussed in this example setting. Figure 11 shows that the sudden expansion of the cross section will cause a flow separation. Due to the sudden change in the geometric boundary, the flumes are separated from the main stream. Their directions are also changed, eventually forming a recirculating pattern. This phenomenon is recognized as the backward-facing step flow, which produces dead zones [47]. In this case, the position and size of the dead zones are related to the degree of expansion. A larger dead zone area will usually be generated after the larger sudden expansion of the cross section. The substances transported upstream will have a certain degree of concentration accumulation in the dead zone. With a larger dead zone, more concentration accumulates, and the flume will contain a greater reaeration coefficient.
(a)
(b)
The DO distributions at the x = 225 m, 450 m, 675 m, and 900 m cross sections are also output for comparison in Figure 12. The 3/4 expanded section is selectable for its largest dead zone, and the accumulation capacity in its dead zone is stronger than those in other cases. On the other hand, the DO concentration downstream has decreased in these 3 cases, and the downstream has a higher reaeration coefficient than when being expanded before. Thus, expansion is conducive for the reaeration in flumes, but the dead zone will also accumulate contaminants, and, downstream of test 8, which has a more uniform reaeration coefficient, the 1/2 expanded channel may have better self-purification ability.
(a)
(b)
(c)
(d)
In summary, the sudden shrinkage of the cross section will hinder the transport of pollutants and DO along the main stream but will not cause a remarkable accumulation in shrunken areas. The sudden expansion of the cross section will cause flow separation, resulting in a backward-facing step flow, and generate dead zones due to the recirculation effect. Stable mass accumulation in dead zones may deteriorate the environment, but at the same time the dead zones will also accumulate DO and drive the water body to have a strong reaeration capacity.
4.3. Analysis of Natural Terrains
In naturally formed mountain rivers or urban rivers, due to the irregularity and complexity in the morphology and hydraulic conditions at both riverbanks and riverbeds, the transport and distribution of dissolved oxygen are also quite complicated [44]. This complex flow pattern will force more backward-facing step flows and dead zones with different scales and shapes. Dead zones accumulate and store solutes transported by the main stream, some of which are later released into the main flow. These dramatic changes are more frequent in natural rivers, making the investigation and prediction of DO transport in large-scale terrain more necessary. Because the coupling model is based on the average water depth, it has a satisfactory efficiency and accuracy in modeling the DO transport in a large-scale topography. At the same time, the model can also consider the influence of natural geometric characteristics and turbulence in calculating the DO distribution and reaeration in water. Therefore, two real terrain cases are considered. In Section 4.3.1, a river in Palong Tsangpo, Xizang, China, which is a mountainous area with steep slopes, is studied. In Section 4.3.2, an additional river basin is also considered, which is the Jinjiang River Basin in Chengdu, China.
4.3.1. Case Study 1: Palong Tsangpo Mountain River
Figure 13(a) shows the terrain of the Palong Tsangpo mountainous area in China. Figures 13(b) and 13(c) are zoomed-in zones, which illustrate the topographic contour and the water depth of the initial dam in Palong Tsangpo, respectively. The computational domain of x = 15540 m and y = 14600 m is taken, which is discretized into 518 × 730 uniform grids with average grids and and time step (the grids have been selected under the principle of convergence and the CFL condition). The boundaries are set as free outflow. A 600 m × 600 m water-bearing dam is set with a height of 200 m and a DO concentration of 5 mg/L. After an instantaneous break, an oxygen-free water source is installed upstream of the dam to an input flow of 400 m³/s to accelerate the transport of oxygen-containing water. In the following description, dz represents the water depth (m), and represents the DO concentration (mg/L).
(a)
(b)
(c)
The parameters in Table 5 are selected mainly according to the convergence and performance, grid size, and time step needed to meet the CFL conditions. Thus, setting smaller grid sizes is more conducive to acquiring accurate computing results in the CFD. However, the computing results will converge with a certain grid size, and the computing result will not be significantly improved while unnecessarily wasting computing resources. The water volume and DO content in the dam will also affect the computing result. The larger water volume will cause the water depth to rise higher in the early stages, and the DO content in the dam can affect the concentration in dead zones. A higher Q will accelerate the flume and mass transfer, raise the water depth, and dilute the DO. The criterion for setting the parameters is to obtain a clear observation.
The entire dam-break process is demonstrated in Figure 14. The process can be divided into 4 stages. In the first stage (t = 2 min and t = 4 min), as the dam breaks in an instant, a rapid change is produced and fluid with high velocity spreads rapidly. During the second stage (t = 10 min and t = 20 min), the water propagates down the ravine between the mountains and reaches the bifurcation position, while the upstream depth remains unchanged. In the third stage (t = 40 min and t = 80 min), the river flow is divided into two streams, and the flow rate begins to decrease. The upper part continues to flow, but because the lower part of the riverbed has an inverted slope, the height gradually increases, and thus the water depth decreases. In the ending stage (t = 120 min and t = 160 min), although there is still a flow input from the upstream, the water depth in this part no longer increases significantly, and the water depth develops into a steady state.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Except for the water depth characteristics, the DO transfer process is further investigated, as demonstrated in Figure 15. In the first stage (t = 2 min and t = 4 min), as the dam breaks the strong turbulence mixes with the upstream nonoxygen stream, resulting in a DO redistribution. The upstream water flow produces a huge amount of kinetic energy in an instant, causing part of the water body and DO to climb along the mountain. During the second stage (t = 10 min and t = 20 min), DO gathers and spreads at the front end of the water flow. In the third stage (t = 40 min and t = 80 min), the water flow spreads to the bifurcation point, and DO propagates in different directions, forming separated flows. The lower part has a higher oxygen concentration due to the lower water depth, and its velocity is relatively low. The flow rate of the upward part is slightly higher, and a dead zone begins to appear at the end of this stage. In the last stage (t = 120 min and t = 160 min), with the stumbling effect of the mountainous area at the bifurcation, a more obvious dead zone appears and continues to develop, and more dead zones with different scales are developed due to the river channel’s complex morphological characteristics. The DO gradually reaches a steady distribution state in this period. In general, the central stream contains a low DO concentration, it is not hindered by the river boundary, and the flow velocity accelerates the transport. However, the velocity is low in the reverse slope area with a low flow velocity, as are the backside of the mountain and expansion positions in the channel. There will be a concentration accumulation in the dead zones, and the result is consistent with the test results in the open channel. In this large-scale computation, uncertainty may be related to these complex geometric conditions and turbulence effects.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
4.3.2. Case Study 2: Jinjiang River Basin
The downstream impact of the DO transport in rivers in densely populated plains is analyzed in this section. Figure 16(a) shows the terrain of the Jinjiang River Basin in the Chengdu Plain, China. Figures 16(b) and 16(c) are the corresponding zoomed-in views, demonstrating the contour and the water depth of the initial dam in the Jinjiang River Basin, respectively. The computational domain is selected as x = 3675 m and y = 5862.5 m. The entire domain is discretized into 294×469 uniform grids, with average grid sizes of and and time steps (the grids have been selected under the principle of convergence and CFL conditions). The boundaries are set as a free outflow. A 15 m high water-bearing dam with a DO concentration of 5 mg/L is set at the upper channel. An oxygen-free water source is set upstream of the dam to input a flow rate of 250 m³/s to accelerate the water shifts.
(a)
(b)
(c)
The parameters in Table 6 are selected mainly according to the convergence and performance, grid size, and time step needed to meet the CFL conditions.
Figure 17 illustrates the process of water evolution after the dam breaks. The entire process can be divided into 4 stages. Since the Jinjiang River Basin in the Chengdu Plain is relatively flat, the height difference is not significant, and the water volume contained by the dam is also low, so there are no drastic downstream changes like in mountain rivers. In the first stage (t = 4 min), the dam breaks instantaneously, and the potential energy is transformed into kinetic energy, which brings about rapid changes and the diffusion of water bodies along the river channel in a short amount of time. In the second stage (t = 20 min), the water flow develops along the curve-shaped channel, while more flumes accumulate in gullies. In the third stage (t = 40 min), the curved pattern continues to grow along the front of the stream. Since the river channel downstream of the dam has been bent, the corresponding change in the curved flow pattern is also reflected in the water depth contour. In the last stage (t = 60 min), because the water body is relatively flat downstream and the islands and beaches are more intricately distributed, the water body begins to spread to lower elevations in a scattered form.
(a)
(b)
(c)
(d)
The water depth changing process as well as the DO transfer process is analyzed, as shown in Figure 18. In the initial stage (t = 4 min), the strong kinetic effect and momentum exchange force the DO to spread to adjacent gullies and shoals to force dead zones. The DO might accumulate in these shoals in the early stages, but they have little change with time. Most of the DO converges in front of the main stream and is transported along the main riverbed. In the second stage (t = 20 min), as the dam is not far from the corner of the river course, the diffusion of DO in water also spreads downward with a tortuous pattern. The riverbed has little height difference, and the terrain is relatively fragmented, so the DO will also spread into the shallows along with the water body. In the third stage (t = 40 min), the tortuous terrain on both sides significantly promotes the redistribution of DO. It can be observed that the DO in separated shallows or gullies remains unchanged. In the final stage (t = 60 min), the terrain downstream is relatively fragmented, and the water depth decreases, so the DO distribution is quite scattered, and as a whole it is distributed in a fan-shaped downstream open channel.
(a)
(b)
(c)
(d)
5. Conclusions
In this paper, coupled SWE and Streeter–Phelps models were used to investigate reaeration. The accuracy of the proposed model in predicting the distribution of pollutants and DO was verified by comparing the analytical results to numerical ones. On this basis, the sudden shrinkage and expansion, the geometric changes in the river channel, and their effects on DO distribution were studied. This study determined that a sudden change in the cross sections will lead to the redistribution and dispersion of the DO due to turbulent mixing. The initial concentration, water volume, and velocity can also affect the final concentration and formation of dead zones. The large change in the reaeration of the sudden shrinkage position is mainly due to the advection and diffusion effects and geometric conditions. The expansion has little effect on the upstream distribution. Although the downstream geometry forms a backward-facing step flow and further produces a dead zone, the total dispersion of DO has not been significantly hindered. Thus, river expansion will strengthen the river’s self-purification ability.
In the computational results of the large-scale terrains, the DO distribution also corresponded with the expansion/shrinkage tests. The complex geometric characteristics of mountains will make the entire DO distribution change significantly, and dead zones are often located at the junctions of rivers or after expanding areas. The DO distribution in plain areas is more scattered and fan-shaped due to the flat and continuous terrains. This model is proven to have calculation and prediction functions for the temporal and spatial changes in DO distribution and reaeration in large-scale real terrains with steady statistical results.
Our model currently not only considers a relatively simple water quality but also builds a solid foundation for establishing a more comprehensive and multi-index water quality model with a reliable hydrodynamic simulation performance in the future.
This well-performed model enlightens us about the feasibility of numerical investigations for predicting the mass transport before construction, ecological restoration, and river improvement to provide a reference for early phase planning. This may reveal the possible impact on the self-purification capacity in the design stage of hydraulic structures. Furthermore, a numerical investigation can be a promising method for assessing the environmentally friendly performance and risks of artificial constructions, such as hydropower stations, groins, artificial islands, and dams. It also enables us to analyze the DO and pollutant distribution, locate regions with oxygen deficits, further formulate remediation plans, and construct improvements to raise the self-purification capacity of these regions and reduce pollution risks. Designers may adopt more reasonable construction plans to strengthen the self-purification capacity of water bodies and choose a structure that is conducive to reaeration and pollutant dissipation.
Data Availability
The water surface and concentration data used to support the findings of this study are available from the corresponding author upon request.
Disclosure
This paper was further developed on the basis of the master dissertations of Jiaxin Wu (first author) and Xiaoxiang Yu (second author) at Sichuan University.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
Dr. Xin Jin and Dr. Jinbo Tang offered great help with the data sources, ideas, and revision of this research.