Numerical Simulation of Hydrodynamics and Reaeration over a Stepped Spillway by the SPH Method

Aerated flows are characterized by complex hydrodynamics and mass-transfer processes. As a Lagrangian method, smoothed particle hydrodynamics (SPH) has a significant advantage in tracking the air-water interface in turbulent flows. This paper presents the application of an SPH method to investigate hydrodynamics and reaeration over stepped spillways. In the SPH method, the entrainment of dissolved oxygen (DO) is studied using a multiphase mass transfer SPH method for reaeration. The numerical results are compared with the hydrodynamics data from Chanson and DO data from Cheng. The simulation results show that velocity distribution and the location of free-surface aeration inception agree with the experimental results. Compared with the experimental results, the distribution of DO concentration over the stepped spillway is consistent with the measurement results. The study shows that the two-phase DO mass transfer SPH model is reliable and reasonable for simulating the hydrodynamics characteristics and reaeration process.


Introduction
Predicting the aeration process has great significance for assessing the dissolved oxygen (DO) concentration, which is critical for the self-purification of polluted water [1]. An accurate understanding of the reaeration process plays an important role in the predicting and improving the DO content in a water body. The stepped spillway, as a type of important hydraulic structure, exhibits the dual functions of energy dissipation and the reaeration process. Research on the stepped spillway has important theoretical and applied value in hydraulic engineering and environmental hydraulics.

Previous Research on the Stepped Spillway
The stepped spillway, which has dual functions of energy dissipation and reaeration, has been used for many years. The initial application of the stepped spillway was limited to overflow energy dissipation [2,3]. Later, researchers studied the hydraulic properties and air entrainment characteristics of the stepped spillway [4,5]. However, few researchers paid attention to the mass transfer process between air and water over the spillway, until such analyses were conducted by Moog [6], and experiments were conducted by Cheng [7,8]). Considering the high turbulence and strong aeration on the spillway surface, it is difficult to accurately measure the air concentration, velocity, and dissolved gas concentration by using experimental methods, especially when the flow depth is too shallow. Some numerical methods based on the Euler grid were developed to simulate the complicated hydrodynamics, air entrainments, and mass transfer characteristics over the stepped spillway [7][8][9]. For the Euler method, the volume fraction near the interface grid is not connective and has strong discontinuities, which leads to parametric oscillation. Therefore, the validity of numerical simulation based on the Euler grid method is challenged by the large deformations of free surface and violent breakup between the interface.
Due to the inherent characteristics of the Euler method, the reaeration process that accompanies the violent interactions between the air-water interface does not show good consistency with the experimental result and full traceability of the individual particles. Accordingly, new numerical methods beyond the mesh-based method should be explored to improve the simulation accuracy for violent flows.

Smooth Particle Hydrodynamics Method
As one type of Lagrangian meshless method, the smoothed particle hydrodynamics (SPH) method has various advantages when considering large deformations and breaking flows [10,11]. To date, different SPH models have been widely studied and applied in the field of fluid mechanics. Recently, SPH methods have been widely applied in single-phase flows [12][13][14][15], multiphase flows [16][17][18], and diffusion processes between different density flows [19,20]. Compared to the Euler method, the existing Lagrangian SPH method has performed well in simulating the motion characteristic of particles. Research on mass transfer between different phases (especially the gas dissolution and mass transfer process), however, is still lacking. Hence, to extend the application of the SPH method, a multiphase mass transfer SPH model for reaeration needs to be studied. For diffusion processes with the SPH method, some studies have been reported. Cleary and Monaghan [21] adopted the SPH method to simulate heat conduction. Tartakovsky et al. [19] developed an SPH model to simulate the diffusion process with a small density difference in mineral precipitation. Aristodemo et al. [20] proposed a model to simulate the convective diffusion of soluble pollutants in water with a density ratio ranging from 0.1 to 1. All of these previous studies focused on diffusion processes through molecular diffusion between liquids. No studies based on the SPH model have been developed to investigate the mass transfer processes between different phases with air and water driven by transportation, turbulent diffusion, and other reactions.
In the following sections, the discrete hydrodynamic equations and mass transfer equation for DO with SPH are illustrated. The model framework and simulation method are described. The experiments conducted by Chanson et al. [5] and Cheng [7,8] are employed to validate the model.

SPH Model for Reaeration
The main features of the SPH method were described by Monaghan [22][23][24] and Liu et al. [25]. Details of the numerical model and discretization of the equations are described below. In the SPH method, the Function A (scalar, vector or tensor) at a certain point r can be defined: where r i , r j are positions; Ω is the support domain; and h is the characteristic length that is called the smoothing length and calculated by h = σ · dp (σ is the dimensionless coefficient, dp is the particle distance). W(r i − r j , h) is the kernel function and should satisfy several conditions, such as compact support, positivity, monotonically decreasing with increasing distance, and the Dirac function property. The kernel function used in this paper is the Wendland kernel function [22]. Function A in Equation (1) can be approximated by using a certain number of particles with initial physical states at different positions [22]: where j is a neighboring particle of particle i, m j is the mass of particle j, ρ j is the density of particle j (j = 1, 2, . . . , N), and N is the number of particles in the support domain Ω. In this paper, the support domain is 2h. Considering the derivative of the kernel function, the gradient of function A can be written as: where ∇ i W r i − r j , h is the derivative of the kernel function W and the equation W(r i − r j , h) can be simplified as W ij .

Hydrodynamic Equations
The SPH form of the continuity equation is described by Monaghan [22]: where the variable u is the particle velocity at the current position. Following Colagrossi et al. [16], the momentum equations can be discretized as air and water phases separately in the two-phase flow.
The SPH form of momentum equation (water phase) is: The SPH form of momentum equation (air phase) is: The Lagrange acceleration is on the left and, on the right, the first term p is the influence of pressure for acceleration, Π is the dissipation item, g is the body force, and f c is the cohesive force, for the water phase f c = 0.
The SPH method is used to handle the problems of initially low or no dissipation. In SPH theory, Monaghan [22] proposed artificial viscosity Π ij , which has been widely used. This not only converts kinetic energy into heat energy, providing the necessary dissipation for the shockwave, but also prevents non-physical penetration of particles for different phases. In this paper, artificial viscosity Π ij is used to solve Equations (5) and (6): Variables c ij , µ ij , and ρ ij in Equation (7) can be rewritten as follows: where α is related to volume viscosity; the term related to β can prevent particle penetration under the condition of high Mach number; α ranges from 0.01 to 0.1 [20]; and u ij = u i −u j , r ij = r i −r j ; c is the sound velocity, which can be calculated by the state equation; and η = 0.1h, which can prevent the divergence of numbers when the particles are too close.
The XSPH method proposed by Monaghan [26] is used to correct the particle motion to ensure a more ordered flow and prevent physical penetration when the particles come closer together: where ε is a specific parameter ranging from 0 to 1.0; and 0.3 is adopted in this paper. Following Monaghan [23], the particle pressure is unique and determined by the density value through a modified Tait's state equation [22]: where B is the reference pressure and p b is the background pressure. Here, B a = B w = ρ 0w c 2 0w /γ w . ρ 0w and ρ 0a are the initial density for water and air, respectively. γ is the isentropic expansion coefficient (γ w = 7, γ a = 1.4). The subscript a is for the air phase; the subscript w is for the water phase. The assumed "reference pressure B" is used to solve the pressure field distribution and ensure both the compressibility of the flow field and the uniformity of particle distribution. In the initial conditions, c 0w can be calculated as follow [27]: where U wmax is the maximum field velocity of the water phase and L w is the water depth of flow field. Mach number is less than 0.1; here, it is (U wmax /c 0w ) 2 < 0.01. For Equation (10), α c is the cohesion coefficient; this represents the cohesion of the same type of particles and is set to 0 for the water phase. The term α c ρ 2 can prevent air phase diffusion and broken interfaces. As suggested by Nugent et al. [28], Equation (6) can be transformed into:

The Advection-Diffusion Equation for DO
The advection-diffusion equation for DO is introduced into the model based on two-phase hydrodynamic equations: where C is the concentration of DO, S C is the source term of DO, ν is the molecular diffusion coefficient, ν t is the turbulent diffusion coefficient, and σ t is the Schmidt number. A different σ t is adopted to conduct the simulation for the two-phase flows. In this paper, σ t is set to 1.0 [29]. Equation (13) is discretized using the SPH operator in the following section.
(1) SPH Discretization for the Diffusion Term Similar to the heat conduction equation proposed by Cleary et al. [21], the heat transfer coefficient is replaced with the diffusion coefficient to simulate the concentration of DO. To simplify the diffusion equation, the diffusion coefficient can be defined as: Following the method developed by Aristodemo et al. [20], the diffusion equation can be discretized into the following form: (2) SPH Discretization for the Convection Term The SPH operator is used to discrete the convection term: In Equation (16), the interface is not stable when the density ratio is less than 0.1. The replacement of m j /ρ i with m j /ρ j is suggested to prevent strong density gradients specifically at the fluid interface [20].
(3) SPH Discretization for the Source Term The complicated surface structure of the stepped spillway can cause collisions between the water and step, which can entrain the surrounding air in the turbulent flows to form dissolved oxygen. With the SPH method, when the water particle contacts air particles, both the contact area between the air and water and the bubble pressure will increase. As a result, the mass transfer between air and water will also increase. By referencing the classical mass transfer theory, the source term based on the SPH method can be expressed as follows: where K is the reaeration coefficient; A t is the specific surface area of bubbles in the unit volume water within a search radius h; C s is the saturated dissolved oxygen concentration under the atmospheric pressure; and C is the dissolved oxygen concentration in water. According to the mass transfer theory, when there is no air particle in the search radius, the mass transfer calculation is not conducted. Instead, the mass transfer will be conducted. For simplification, assuming that the bubble diameter is consistent in the mixing flow, the specific surface area of the bubble in unit water volume can be expressed as: where N a is the number of air particles in the search radius; A b is the surface area of a single bubble; and V w is the volume of the water. On two-dimensional numerical simulation, the variables can be rewritten as follows: According to the relationship between the dissolved gas mass transfer coefficient and the Reynolds number proposed by Banerjee et al. [30], the reaeration coefficient K can be calculated using: where θ is recommended as 0.2 and σ t is a constant value; Re is the effective Reynolds number using the formula proposed by Price [31].

Time Integration
To integrate the scheme forward in time, the Verlet explicit scheme is used [32]. Under the Courant Friedrichs Levy (CFL) condition, the unit mass force and diffusion condition derived from the viscous term are used to determine the variable time step. Detailed steps are described here.
As suggested by Monaghan [22], the time step, ∆t, of unit mass force can be described as: where f i is the acceleration of the unit mass force. Following Morris et al. [33], the time step, ∆t, caused by viscosity is: Computation of the standard time step integrates the viscosity and unit mass force effect: where C CFL is between 0.1 and 0.2. In the paper, C CFL is set equal to 0.2.

Model Framework and Simulation Method
The SPH model simulation developed in this paper contains three main steps, as shown in Figure 1. (I) Data Preparation: Before the program is executed, the data must be transformed into a format that the software can identify. The data include the particle properties, the execution parameters and the calculation domain properties; (II) Main Loop: After the preparatory works are executed, the interaction of particles is executed using hydrodynamic and mass transfer calculations. This step is the core method for determining the computing time of the SPH model; and (III) Save and Analysis: The mechanical property and mass transfer process will be saved by a time series. In the analysis software, the results will be rendered and analyzed.
where fi is the acceleration of the unit mass force. Following Morris et al. [33], the time step, Δt, caused by viscosity is: Computation of the standard time step integrates the viscosity and unit mass force effect: where CCFL is between 0.1 and 0.2. In the paper, CCFL is set equal to 0.2.

Model Framework and Simulation Method
The SPH model simulation developed in this paper contains three main steps, as shown in Figure  1. (I) Data Preparation: Before the program is executed, the data must be transformed into a format that the software can identify. The data include the particle properties, the execution parameters and the calculation domain properties; (II) Main Loop: After the preparatory works are executed, the interaction of particles is executed using hydrodynamic and mass transfer calculations. This step is the core method for determining the computing time of the SPH model; and (III) Save and Analysis: The mechanical property and mass transfer process will be saved by a time series. In the analysis software, the results will be rendered and analyzed. In each step of the SPH method, the interactions between different particles are calculated using traversal paired searches for certain particles in the computation domain, and this leads to an increased computation time. In this paper, a cell-linked list search method is used to construct the computation framework, and the computation domain is divided into "background cells" shown in Figure 2 [34]. The cell aids the numerical simulation but does not participate in it. For 2D coordinates, the cell size is the length of the support domain (2h). All of the pair-wise particle distances are less than 2h, which means all adjacent particles can only appear in the cell of the adjacent cell. For a particle that is located inside a cell, only the interactions with particles of neighboring cells must be In each step of the SPH method, the interactions between different particles are calculated using traversal paired searches for certain particles in the computation domain, and this leads to an increased computation time. In this paper, a cell-linked list search method is used to construct the computation framework, and the computation domain is divided into "background cells" shown in Figure 2 [34].
The cell aids the numerical simulation but does not participate in it. For 2D coordinates, the cell size is the length of the support domain (2h). All of the pair-wise particle distances are less than 2h, which means all adjacent particles can only appear in the cell of the adjacent cell. For a particle that is located inside a cell, only the interactions with particles of neighboring cells must be considered. Therefore, during the force and mass transfer computation, domains without adjacent particles should not be calculated, and this greatly improves the computational efficiency. considered. Therefore, during the force and mass transfer computation, domains without adjacent particles should not be calculated, and this greatly improves the computational efficiency. Judgment of air-water interfaces is the premise for studying the mass transfer of dissolved gas. In the paper, the range of the mass transfer between air and water is limited to a search radius. Before computing and analyzing the mass transfer, the particle types of different phases should be defined with different symbols. In this paper, the critical function is employed within the scope of the search radius.
The meaning of function Ce is when i and j are different phases within the search radius, the mass transfer process is computed. Otherwise, i and j are the same phases, and the mass transfer is set to 0.
In the paper, the solid boundaries are described by a set of discrete particles which are fixed in position. The momentum equations of the boundary particles are the same as water particles in the simulation. In addition, the positions of the boundary particles do not change according to the imposed forces on them, and the computer used in the simulation is an Intel ® Core™ i7-4790 CPU @ 3.60 GHz; (RAM) 16.0 GB; Windows 7. By running the case for hydrodynamics simulation for 8 s, the runtime is about 9 h.

Validation for the Hydrodynamics over the Stepped Spillway
The hydrodynamic model was validated using the experiments carried out by Chanson et al. [5]. In this experiment, the properties of air-water flow were measured with a probe. Here, the computation domain was the same as that used in the experiment.  Judgment of air-water interfaces is the premise for studying the mass transfer of dissolved gas. In the paper, the range of the mass transfer between air and water is limited to a search radius. Before computing and analyzing the mass transfer, the particle types of different phases should be defined with different symbols. In this paper, the critical function is employed within the scope of the search radius.

Description of Chanson's Experiment
The meaning of function Ce is when i and j are different phases within the search radius, the mass transfer process is computed. Otherwise, i and j are the same phases, and the mass transfer is set to 0. In the paper, the solid boundaries are described by a set of discrete particles which are fixed in position. The momentum equations of the boundary particles are the same as water particles in the simulation. In addition, the positions of the boundary particles do not change according to the imposed forces on them, and the computer used in the simulation is an Intel ® Core™ i7-4790 CPU @ 3.60 GHz; (RAM) 16.0 GB; Windows 7. By running the case for hydrodynamics simulation for 8 s, the runtime is about 9 h.

Validation for the Hydrodynamics over the Stepped Spillway
The hydrodynamic model was validated using the experiments carried out by Chanson et al. [5]. In this experiment, the properties of air-water flow were measured with a probe. Here, the computation domain was the same as that used in the experiment.

Set-Up Parameters
The parameters used for simulation of hydrodynamics are presented in Table 1.  Figure 4 shows the simulation result of the distribution of different types of particles over the stepped spillway. As shown in Figure 4, the simulation result is generally consistent with the experimental result in terms of the air water distribution and flow pattern, see Figure 3-1(B) in Chanson et al. [35]. Here, a transition flow pattern can be observed. In the third step, the air-water interface is clearly broken up and deformed. The mixed process becomes severe, and the inception point of free-surface aeration clearly appears. The numerical simulation matches the experimental. Downstream of the inception point of free-surface aeration, the kinetic energy is generated. Under this condition, water particles overcome gravity and jump away from the water surface. Thus, the mixed phenomenon can be seen clearly. Meanwhile, the nappe reattaches the water surface at the next step, and some deflected nappes fall back to the fifth and sixth step edges. In addition, when the flow passes the step edges, a certain type of fluctuation over the steps can be observed downstream of the inception point of free-surface aeration. The simulation results are in accordance with the experimental results and correctly reflect the mixing situation of the laboratory experiment. √ gq (q is the discharge per unit width, q = Q/w, where w is the width of the channel). In this paper, the discharge per unit width q is 0.058 m 3 /s; h c = 3 q 2 /g is the critical depth, where h s (m) is the step height, l s (m) is the step length, and H s (m) is the water elevation in the tank.

Set-Up Parameters
The parameters used for simulation of hydrodynamics are presented in Table 1.  Figure 4 shows the simulation result of the distribution of different types of particles over the stepped spillway. As shown in Figure 4, the simulation result is generally consistent with the experimental result in terms of the air water distribution and flow pattern, see Figure 3-1(B) in Chanson et al. [35]. Here, a transition flow pattern can be observed. In the third step, the air-water interface is clearly broken up and deformed. The mixed process becomes severe, and the inception point of free-surface aeration clearly appears. The numerical simulation matches the experimental. Downstream of the inception point of free-surface aeration, the kinetic energy is generated. Under this condition, water particles overcome gravity and jump away from the water surface. Thus, the mixed phenomenon can be seen clearly. Meanwhile, the nappe reattaches the water surface at the next step, and some deflected nappes fall back to the fifth and sixth step edges. In addition, when the flow passes the step edges, a certain type of fluctuation over the steps can be observed downstream of the inception point of free-surface aeration. The simulation results are in accordance with the experimental results and correctly reflect the mixing situation of the laboratory experiment.   (Figure 5b), the breaking degree of water increases constantly, which is typical for aeration and entrainment. Inside the steps, an expected clockwise vortex appears, and the velocity of water increases from center to the edge, which has a significant effect on energy dissipation. The main reason for velocity change inside the steps is that part of the kinetic energy of the flow will turn into heat energy, leading to velocity attenuation and energy dissipation.  Figure 6 shows the distribution of water velocity under the condition of dimensionlessness from Step 3 to Step 8. The abscissa is U/Uc, and the ordinate is y/hc, which are Uc = 0.668 m/s and hc = 0.07 m, respectively, here. As shown in Figure 6, the numerical simulation results based on the two-phase SPH model agree well with the experimental results, except the result in Step 3, and the hydrodynamic characteristics can be correctly reflected. In Step 3, the location of the inception of freesurface aeration appears. In this step, the air-water interface edge breaks up and deforms violently, which may lead to difficulties with measurement and errors in the simulation. In general, the velocity distributions vary rapidly at different step edges. The flow velocity is larger in the interval [0.3hc, 0.8hc] than it is at other locations. When y > 0.8hc, the velocity of monitoring points fluctuates significantly. Small parts of water drop much faster than do others due to the splash effect of the previous steps. Some parts of water velocity exhibit attenuation due to the collision of the step itself, which converts kinetic energy into potential energy.

Discussion about the Hydrodynamics Characteristics
Compared to the experimental data, the two-phase SPH model presented in this paper can reflect accurately and reliably the hydrodynamic characteristics.  Figure 5 shows the velocity field of the water phase on Step 2 and 5. The results obtained from analyzing the velocity field of the water on the stepped spillways are more intuitive. Upstream of the inception point of free-surface aeration (Figure 5a), an obvious cavity can be observed. Downstream of the inception point of free-surface aeration (Figure 5b), the breaking degree of water increases constantly, which is typical for aeration and entrainment. Inside the steps, an expected clockwise vortex appears, and the velocity of water increases from center to the edge, which has a significant effect on energy dissipation. The main reason for velocity change inside the steps is that part of the kinetic energy of the flow will turn into heat energy, leading to velocity attenuation and energy dissipation.  Figure 5 shows the velocity field of the water phase on Step 2 and 5. The results obtained from analyzing the velocity field of the water on the stepped spillways are more intuitive. Upstream of the inception point of free-surface aeration (Figure 5a), an obvious cavity can be observed. Downstream of the inception point of free-surface aeration (Figure 5b), the breaking degree of water increases constantly, which is typical for aeration and entrainment. Inside the steps, an expected clockwise vortex appears, and the velocity of water increases from center to the edge, which has a significant effect on energy dissipation. The main reason for velocity change inside the steps is that part of the kinetic energy of the flow will turn into heat energy, leading to velocity attenuation and energy dissipation.
(a) (b) Figure 5. The velocity field of the water phase. (a) the velocity field of the water phase on step 2, (b) the velocity field of the water phase on step 5. Figure 6 shows the distribution of water velocity under the condition of dimensionlessness from Step 3 to Step 8. The abscissa is U/Uc, and the ordinate is y/hc, which are Uc = 0.668 m/s and hc = 0.07 m, respectively, here. As shown in Figure 6, the numerical simulation results based on the two-phase SPH model agree well with the experimental results, except the result in Step 3, and the hydrodynamic characteristics can be correctly reflected. In Step 3, the location of the inception of freesurface aeration appears. In this step, the air-water interface edge breaks up and deforms violently, which may lead to difficulties with measurement and errors in the simulation. In general, the velocity distributions vary rapidly at different step edges. The flow velocity is larger in the interval [0.3hc, 0.8hc] than it is at other locations. When y > 0.8hc, the velocity of monitoring points fluctuates significantly. Small parts of water drop much faster than do others due to the splash effect of the previous steps. Some parts of water velocity exhibit attenuation due to the collision of the step itself, which converts kinetic energy into potential energy.
Compared to the experimental data, the two-phase SPH model presented in this paper can reflect accurately and reliably the hydrodynamic characteristics.  Figure 6 shows the distribution of water velocity under the condition of dimensionlessness from Step 3 to Step 8. The abscissa is U/U c , and the ordinate is y/h c , which are U c = 0.668 m/s and h c = 0.07 m, respectively, here. As shown in Figure 6, the numerical simulation results based on the two-phase SPH model agree well with the experimental results, except the result in Step 3, and the hydrodynamic characteristics can be correctly reflected. In Step 3, the location of the inception of free-surface aeration appears. In this step, the air-water interface edge breaks up and deforms violently, which may lead to difficulties with measurement and errors in the simulation. In general, the velocity distributions vary rapidly at different step edges. The flow velocity is larger in the interval [0.3h c , 0.8h c ] than it is at other locations. When y > 0.8h c , the velocity of monitoring points fluctuates significantly. Small parts of water drop much faster than do others due to the splash effect of the previous steps. Some parts of water velocity exhibit attenuation due to the collision of the step itself, which converts kinetic energy into potential energy.
Compared to the experimental data, the two-phase SPH model presented in this paper can reflect accurately and reliably the hydrodynamic characteristics.

Validation for Reaeration over the Stepped Spillway
The reaeration process over the step spillway was validated using the experiment carried out by Cheng [7,8]. The computation domain was similar to the physical model with Cheng's experiment.

Description of Cheng's Experiment
In the experiment, clear water that meets the required oxygen deficit condition is used as the experimental water. Figure 7 shows the structure of the stepped spillway. There are seven points at which the DO concentration is monitored. As shown in Figure 7, point 1 is the monitoring point for the initial DO concentration at the inlet, points 2-6 are the monitoring points on the stepped spillway surface and point 7 is the monitor point at the toe.

Validation for Reaeration over the Stepped Spillway
The reaeration process over the step spillway was validated using the experiment carried out by Cheng [7,8]. The computation domain was similar to the physical model with Cheng's experiment.

Description of Cheng's Experiment
In the experiment, clear water that meets the required oxygen deficit condition is used as the experimental water. Figure 7 shows the structure of the stepped spillway. There are seven points at which the DO concentration is monitored. As shown in Figure 7, point 1 is the monitoring point for the initial DO concentration at the inlet, points 2-6 are the monitoring points on the stepped spillway surface and point 7 is the monitor point at the toe.  Table 2 shows the structural characteristics of the stepped spillway. Where the WES curve is the surface equation of the weir crest; the weir crest meets the standard WES curve; design slope is the design grade of the spillway; T is the point of tangency of the WES curve and slope surface; steps is the number of uniform steps; hs and ls are the length and width of the step, respectively; Hs is the water elevation in the tank; q is the discharge per unit width; Cu is the concentration of DO at the weir crest; Cd is the concentration of DO at the toe; and Cs is the saturated DO concentration.

Set-up Parameters
The parameters used in the simulation are detailed in Table 3.   Table 2 shows the structural characteristics of the stepped spillway. Where the WES curve is the surface equation of the weir crest; the weir crest meets the standard WES curve; design slope is the design grade of the spillway; T is the point of tangency of the WES curve and slope surface; steps is the number of uniform steps; h s and l s are the length and width of the step, respectively; H s is the water elevation in the tank; q is the discharge per unit width; C u is the concentration of DO at the weir crest; C d is the concentration of DO at the toe; and C s is the saturated DO concentration.

Set-up Parameters
The parameters used in the simulation are detailed in Table 3.  Figure 8 shows the results obtained for the air water distribution. The free surface line gradually decreases from crest to toe, and the line is smooth. The critical water depth h c is approximately 0.14 m, which is in agreement with the experimental data. Figure 8 (upper right) shows the particle distributions for different phases on the stepped spillway. The mixed air and water flow formed on the step surfaces with violent turbulence and strong aeration. Some incompletely developed vortices that entrain air can also be seen inside the step. Figure 8 shows the results obtained for the air water distribution. The free surface line gradually decreases from crest to toe, and the line is smooth. The critical water depth hc is approximately 0.14 m, which is in agreement with the experimental data. Figure 8 (upper right) shows the particle distributions for different phases on the stepped spillway. The mixed air and water flow formed on the step surfaces with violent turbulence and strong aeration. Some incompletely developed vortices that entrain air can also be seen inside the step.    Figure 8 shows the results obtained for the air water distribution. The free surface line gradually decreases from crest to toe, and the line is smooth. The critical water depth hc is approximately 0.14 m, which is in agreement with the experimental data. Figure 8 (upper right) shows the particle distributions for different phases on the stepped spillway. The mixed air and water flow formed on the step surfaces with violent turbulence and strong aeration. Some incompletely developed vortices that entrain air can also be seen inside the step.  Data comparisons at different points are given in Table 4. As shown in Table 4, the simulation results based on the two-phase SPH mass transfer model agreed with the experiment results well.

Discussion about the Reaeration
The max relative error appeared at point 2 because the air entrainment at this point was not obvious. Another reason for the difference is that the measuring point is far from the water surface, which may lead to some differences between simulation and experiment results. The overall average relative error was 8.8%. Based on the analysis, we concluded that the two-phase SPH model established in the paper is capable of capturing the main features of the DO reaeration. The reaeration rate of DO deficit, r , is used to analyze the computational accuracy of the numerical simulation: The efficiency of the reaeration rate of DO is higher when r is greater. There is no DO reaeration rate when r = 0. Table 5 shows the DO reaeration rate, r , of numerical simulation and experimental results. The relative error is within 6%, which satisfies the accuracy requirement of the concentration field. From the analysis described, the two-phase SPH mass transfer model accurately computes the reaeration of DO on the stepped spillway. The DO concentration distribution in the numerical simulation is in accordance with the experimental data. The relative error of r is within 6%. The comparison results demonstrate the accuracy and reliability of the two-phase SPH mass transfer model.

Conclusions
To analyze the change in hydrodynamic characteristics and the dissolved oxygen concentration distribution over space and time, the advection-diffusion equation of the SPH was established in the paper. A source term was added in the two-phase SPH model. Additionally, the mass transfer coefficient calculation formula that fit the SPH method was adopted. Numerical simulations of the hydrodynamic characteristics and the mass transfer process of the stepped spillway were realized based on the two-phase mass transfer SPH model. Considering the cohesive force and source term, the two-phase mass transfer SPH model could correctly reflect the flow state, the change in velocity, the location of the inception of free-surface aeration, and the DO mass transfer process over the stepped spillway. Following Chanson's experiment, the simulation result showed that the location of the inception of free-surface aeration appeared at the third step edge. At this step, the edge of the air-water interface broke up and deformed clearly. The recirculation vortices could be observed in the downstream of the inception of free-surface aeration. The velocity was significantly greater at the edge attachments of the steps than it was inside the vortex. Following Cheng's experiment, the mass transfer mechanism was as follows: the DO concentration increased from crest to toe, and the vertical concentration gradient decreased gradually. This result implies that the DO diffusion accelerates from the turbulent aerated flows, which facilitates the uniform vertical distribution of the DO concentration.
Compared to the traditional mesh-based methods, the two-phase SPH mass transfer model established in this paper can simulate the large deformation and turbulent breakup flow properties, as well as the mass transfer process under strong turbulence conditions over the stepped spillway well.