Numerical Investigation of Dissolved Oxygen Transportation through a Coupled SWE and Streeter–Phelps Model

Dissolved oxygen (DO) reflects the self-purification ability of a water body and is also an important indicator for quantifying the water quality. ,e 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. ,e 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. ,e influences of the sudden crosssectional changes in the river channel on the DO distribution and the reaeration ability are identified. ,e 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.


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. e 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 selfpurification abilities. Additionally, the turbulent diffusion and transport effects in water streams can dilute the concentration of pollutants. e 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. e 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.
ere 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. e 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. e gas-transfer process is complicated, and research has not yet achieved a complete understanding [6]. With vast laboratory, field, and numerical studies, gastransfer models have developed greatly, such as the twofilm 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][10][11][12][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.

Materials and Methods
e transportation of oxygen and pollutants is determined by hydrodynamics and the transport substances' properties.
e 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. e main outputs are the water depth and the concentration value in discretized grids.
is 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. e shallow water model has been widely used and well developed by many scholars. e 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. e main governing equations have been introduced in their previous research, and they proved that this SWE method is indeed valid. erefore, the hydrodynamic, mass transfer, and discrete algorithms are excluded in this article, and the reaeration equation will be mainly introduced.

Hydrodynamic Model.
e hydrodynamic model was referred to in Zhang and Lin's work [18], which was proven to be reliable and robust in his research. e 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]. us, the depth-averaged k − ε 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][23][24]. e 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 H(� η − z b ) is defined as the total water depth, η is the height of the water surface line relative to the base level, z b is the height of the riverbed, (P � u × H) and (Q � v × H) are the flux components per unit volume in the x direction and y direction, respectively, u and v 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, g is the acceleration due to gravity, and m is the Manning roughness coefficient. e 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. e resulting k − ε 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. e superscript indicates time average value, ] t is the turbulent eddy viscosity ] t � C μ (k 2 /ε), and σ k , σ ε , C ε1 , and C ε2 are model parameters, where σ k � 1.0, σ ε � 1.3, C ε1 � 1.44, and C ε2 � 1.92. e main function of these equations is to compute ] t and further apply it to the following mass transfer equation.

Mass Transfer Equation.
e three-dimensional form of the advection-diffusion equation deduced from Fick's law is where C is the concentration of pollutants or oxygen and D x , D y , and D z are the flow diffusion coefficients. On this basis, it is necessary to derive the mass transfer equation in a twodimensional 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 e 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. erefore, the molecular diffusion equation for turbulent diffusion can be transformed as where ] t 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 Sc t � 1 [25].

Reaeration
Model. e 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 S O 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 C s and C 0 characterize the saturation concentration of oxygen in the water body and the actual concentration at the time, respectively, reflecting the degree of oxygen deficiency, k 2 and k L 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-Mathematical Problems in Engineering transfer model as the source term. In most cases, influencing factors, such as atmospheric conditions and hydraulic slopes, can be ignored [28][29][30][31][32]. However, some factors truly should not be oversimplified or ignored, such as the temperature, Schmidt number (Sc � (]/D)), and turbulence conditions included in this research. erefore, 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. e 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 k 2 is the reaeration coefficient in the reaeration model. e formula is further simplified and imitated, and many semiempirical formulas are obtained, such as the Churchill formula, ackston-Krenkel formula, Benett and Rathbun formula, and Owens formula [8][9][10][11][12].
is model refers to the work of Nguyen et al. [14] who conducted experiments to study the reaeration process. e 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 k 2 can be expressed as follows: where k 2 is the reaeration coefficient of the average water depth. By applying Nguyen's experimental data [14], λ 1 ′ and λ 2 can be represented as is 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 is 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. D m is determined by the Stokes-Einstein equation in the following formula: where M � 2.26, x � 25.6 cm 2 /mol in the water solvent, and V is the molar volume of oxygen.

Discrete Method.
is 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]. ese methods in the model can obtain good calculation results for most cases. If numerical oscillations occur in advectiondominant situations, errors can be avoided by applying the WENO format [35].

Verification
e 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. us, available analytical, numerical, and experimental results will be adopted to verify the present model. eir 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. e 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. erefore, this article will omit the hydrodynamic verification part, mainly focusing on verification of the S-P model and mass transfer.

Streeter-Phelps Model Test.
e 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].
is 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. is model mainly refers to a one-dimensional steady-state case, and the equation of the Street-Phelps model satisfies the following form: where L(x) represents the function of the BOD change, D(x) represents the function of the DO change, k 1 is the oxygen consumption coefficient, k 2 is the reaeration coefficient, and u is the flow velocity. e analytical solution can be obtained as where L 0 is the initial BOD value and D 0 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. e 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. e 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: Q � 0.013 (m 3 /s), A � 0.4 m 2 , and u � 2808 m/day, and the coefficients adopt the settings proposed by Chapra [37]: k 1 � 3 day − 1 and k 2 � 10 day − 1 . e initial condition setting, L 0 � 10 mg/L and D 0 � 8 mg/L, 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 Δx � 2 m and Δy � 1 m, and the grid numbers are 1000 in x and 1 in y. e parameters are listed in Table 1.
e 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. is 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.

Two-Dimensional Gaussian Hump Advection and Diffusion Verification.
e 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. e 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 is , (20) where W represents the concentration of mass, t represents time, u and v are the velocities in the x and y directions, respectively, D is a diffusion coefficient, and x 0 and y 0 are initial coordinates of the hump. e parameters are listed in Table 2. e numerical computational domain of 4 m × 4 m is discretized into 320 × 320 uniform grids, with the average discretized grid sizes being Δx � 0.025 m and Δy � 0.025 m, respectively. e initial concentration is set according to the theoretical solution (20). A fixed time step Δt � 0.002 s is adopted, with parameters x 0 � 0.5 m, y 0 � 0.5 m, and D � 0.01. e initial velocities of the fluids are u � 1.0 m/s and v � 1.0 m/s. e 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 u � 1.0 m/s and v � 1.0 m/s 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. e 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. is 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. is 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.

Dispersion of Transport after Dam Break.
e numerical results of the abovementioned two cases are in agreement with the theoretical results, so this section will focus on the Mathematical Problems in Engineering 5 accuracy of the coupled computation. is 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. e dam-break case has been investigated by many authors [39][40][41][42]. Our settings are in agreement with these, and our results are compared. e 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. e 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 as where W represents the DO concentration (mg/L). e 3D view's initial bottom and oxygen concentration setting are demonstrated in Figure 4. e water depth and concentration on either side are different, and a high concentration is set near the gate at 650 m and 600 m. e 1400 m × 1400 m numerical computational domain is discretized into 500 × 500 uniform grids, with the average discretized grid sizes being Δx � 2.8 m and Δy � 2.8 m, respectively. e fixed time step Δt � 0.2 s is adopted. 3D contours at t � 200 s are plotted in Figure 5. e 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. e eddies in this case have the same size and shape due to the symmetry of the dam. e dead zones cause the concentration on both sides to increase at different degrees because the initial concentration is    Mathematical Problems in Engineering asymmetrical. e lower side with a higher initial concentration will accumulate more DO due to the recirculation effect. 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. e 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]. e 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.

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. e 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. e uneven distribution of DO may cause a local oxygen deficit, weaken the self-purification ability, and cause pollution. erefore, 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. e 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 D m � 2.1824 × 10 − 9 m 2 /s. e average discretized grid sizes are Δx � 5 m and Δy � 2 m. In addition, the fixed time step is Δt � 0.5 s.
e initial concentration C of DO is 3 mg/l, the oxygen saturation concentration C s is set as 9.09 mg/l, and the reaeration coefficient is calculated according to formula (19). e 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. e setting of the tests is summarized in Table 3, and the parameters are listed in Table 4. e parameters are chosen to better visualize the computational results, and a lower velocity is better for observing the reaeration effect in a limited length. e 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. e 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.

Analysis of Sudden Shrink.
e 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. e shrunken position can hinder the flume and solute from transferring along the main stream. e 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.
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. e 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.
e 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.
us, 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.

Analysis of Sudden Expansion.
e 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. eir directions are also changed, eventually forming a recirculating pattern. is   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. e 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. e DO distributions at the x � 225 m, 450 m, 675 m, and 900 m cross sections are also output for comparison in Figure 12. e 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. us, 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.
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. e 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.

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]. is 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. ese 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. erefore, 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.  . e boundaries are set as free outflow. A 600 m × 600 m waterbearing 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 w represents the DO concentration (mg/L). e parameters in Table 5 are selected mainly according to the convergence and performance, grid size, and time step needed to meet the CFL conditions. us, 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. e water volume and DO content in the dam will also affect the computing result. e 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. e criterion for setting the parameters is to obtain a clear observation.  Figure 14. e 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. e 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.

Case Study 1: Palong Tsangpo Mountain River.
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. e 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. e lower part has a higher oxygen concentration due to the lower water depth, and its velocity is relatively low. e 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. e 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  ere 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.

Case Study 2: Jinjiang River Basin.
e 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.
e computational domain is selected as x � 3675 m and y � 5862.5 m. e entire domain is discretized into 294×469 uniform grids, with average grid sizes of Δx � 12.5 m and Δy � 12.5 m and time stepsΔt � 0.02 s (the grids have been selected under the principle of convergence and CFL conditions). e 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. e 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. e 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.
e water depth changing process as well as the DO transfer process is analyzed, as shown in Figure 18. In the   diffusion of DO in water also spreads downward with a tortuous pattern. e 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.

Conclusions
In this paper, coupled SWE and Streeter-Phelps models were used to investigate reaeration. e 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. is study determined that a sudden change in the cross sections will lead to the redistribution and dispersion of the DO due to turbulent mixing. e initial concentration, water volume, and velocity can also affect the final concentration and formation of dead zones. e large change in the reaeration of the sudden shrinkage position is mainly due to the advection and diffusion effects and geometric conditions. e expansion has little effect on the upstream distribution. Although the downstream geometry forms a backwardfacing step flow and further produces a dead zone, the total dispersion of DO has not been significantly hindered. us, 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. e 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. e DO distribution in plain areas is more scattered and fan-shaped due to the flat and continuous terrains. is 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.
is 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. is may reveal the possible impact on the selfpurification 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 e water surface and concentration data used to support the findings of this study are available from the corresponding author upon request. Disclosure is 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
e authors declare that there are no conflicts of interest regarding the publication of this paper.