Optimal In-Stream Structure Design through Considering Nitrogen Removal in Hyporheic Zone

: The hyporheic zone (HZ), the region beneath or alongside a streambed, can play a vital role in a stream ecosystem. Previous studies have examined the impacts of in-stream structures on the HZ and river restoration; however, studies on optimizing the design of in-stream structures are still lacking. Therefore, this study aims to propose a method for optimizing the design of in-stream structures (e.g., weirs) through comprehensively considering both nitrogen removal amount ( NRA ) and nitrogen removal ratio ( NRR ) in the HZ based on numerical modelling. The Hydrologic Engineering Center’s River Analysis System (HEC-RAS) and COMSOL Multiphysics are employed for surface water and hyporheic ﬂow simulations, respectively, and these two models are coupled by the hydraulic head along the surface of the streambed. The NRA and NRR are both closely related with residence time ( RT ), while the NRA is also inﬂuenced by hyporheic ﬂux. Using the model outputs under di ﬀ erent scenarios, regression equations for estimating the relevant variables (e.g., the maximum upstream distance in the subsurface ﬂow inﬂuenced by the weir, the RT , and the hyporheic ﬂux) are proposed. Then, the cumulative NRA ( CNRA ) and NRR can be calculated, and an objective function is formulated as the product of the normalized CNRA and NRR . The results show that the optimal height of the weir can be obtained based on the proposed method, and the validation shows the good general performance of this method. Sensitivity analysis indicates that the optimal height generally can be sensitive to the river discharge, i.e., the optimal height increases when the river discharge increases and vice versa. In addition, it is observed that, in the case of the optimal height, hyporheic ﬂux increases when the slope increases while the inﬂuence of depth to bedrock on hyporheic ﬂux is not signiﬁcant. This study enhances our understanding of the optimal in-stream structure design, and potentially beneﬁts river restoration in the face of continual degradation caused by human activities.


Introduction
The hyporheic zone (HZ), which is the region beneath or alongside a streambed and characterized by the active mixing of shallow groundwater and surface water, acts as an active/connecting ecotone between groundwater and surface water [1][2][3][4][5][6][7]. The HZ forms a unique environment for macroinvertebrates and biogeochemical reactions, and is paramount to the biological and ecological functions of streams [2,8]. For example, the biogeochemical reactions in the HZ produce unique water Summarizing the studies mentioned above, hyporheic flux is always an important indicator for the evaluation of the in-stream structure in hyporheic restoration. Some in-stream structures have been proposed to improve the hyporheic flux. However, most focused on the HZ changes induced by the river restoration structures, rather than optimizing the design of the structures. Besides designing for new structures, retrofitting existing in-stream structures to improve the hyporheic exchange has also been studied. Moreover, the function of the HZ on pollutant removal has been considered in some recent hyporheic restoration projects. However, there is still a lack of simplified methods which can be widely used for preliminarily evaluating the suitability of a site and optimizing the design of structures for hyporheic restoration.
To fill this gap, this study proposes a method to optimize the design of in-stream structures by considering the nitrogen removal in the HZ based on numerical modelling. The weir is selected as the representative of in-stream structures since it is regarded as the most effective in improving the hyporheic flux [15]. Among the different factors that can influence the performance of the weir in hyporheic restoration, the height of the weir is one of the most important. Hence, in this study, the proposed method aims to evaluate the HZ performance influenced by different heights of the weir using several indicators (e.g., the maximum upstream distance in the subsurface flow influenced by the weir, the RT, and the hyporheic flux). It also specifically considers the denitrification function of the HZ. Overall, the proposed methods and the results obtained will enhance our understanding and facilitate future studies regarding the evaluation of the nitrogen removal ratio in the HZ, which can be important to the stream ecosystem and benefit the restoration of the stream and the HZ [4,14,17,30,31].

Materials and Methods
The first step in this study is to build a numerical model that couples surface water and hyporheic flow using the hydraulic head along the surface of the streambed. Second, regression equations to estimate the relevant indicators for the optimal design of the weir are proposed. Third, the height of the weir is optimized, and the results obtained are validated.

Numerical Model
In order to simulate the flow regime which is influenced by the weir, HEC-RAS (Hydrologic Engineering Center's River Analysis System) is employed to simulate surface water flow in this study [32]. Using HEC-RAS, the hydraulic head along the surface of the streambed can be calculated based on the energy equation, and is used as the top boundary condition for subsurface water flow simulation. At the same time, COMSOL Multiphysics is employed to simulate the subsurface flow based on the Brinkman-Darcy equation. Both the surface and subsurface flow simulations are conducted in the steady state.

Hydrologic Engineering Center's River Analysis System (HEC-RAS) for Surface Water Simulation
In the surface water simulation, this study basically follows the study of Hester and Doyle [15]. The length of the domain is 20 m, and the weir is located in the middle of the domain (see Figure 1); the width of the stream is 3 m, and the maximum depth of the stream is 1 m, with the rectangular cross-section as illustrated in Figure 1; the Manning coefficient of the streambed is 0.03. The boundary condition upstream of the weir is set as the calculated critical depth of the flow. In order to detect the exact upstream location at which the depth is critical, the standard step method, which is a computational technique to estimate one-dimensional surface water profiles in open channels with gradually varied flow under steady state conditions, is used in HEC-RAS. In this method, by applying the energy equation from the location of the weir where the water height (above the weir) correlated with the designed stream discharge, the exact upstream location of critical depth can be obtained through an iterative process. In addition, the downstream boundary condition, that is, the right end of the domain, is set as the normal depth (see Figure 1). Based on such information, the hydraulic Water 2020, 12, 1399 4 of 16 head along the surface of the streambed can be calculated when the river discharge is designated. The energy equation for the simulation of surface water is as follows [33]: where Z 1 and Z 2 are elevations of the main channel inverts (m), Y 1 and Y 2 are depths of water at cross sections (m), V 1 and V 2 are average velocities at cross sections (m/s), c 1 and c 2 are velocity weighting coefficients, g is gravitational acceleration (m/s 2 ), and h e is energy head loss (m).
Water 2020, 12, x FOR PEER REVIEW 4 of 16 where Z1 and Z2 are elevations of the main channel inverts (m), Y1 and Y2 are depths of water at cross sections (m), V1 and V2 are average velocities at cross sections (m/s), c1 and c2 are velocity weighting coefficients, g is gravitational acceleration (m/s 2 ), and he is energy head loss (m).

COMSOL Multiphysics for Hyporheic Flow Simulation
Imposing the hydraulic heads along the surface of the streambed as the upper boundary condition, the hyporheic flow is simulated using COMSOL Multiphysics. In this study, the weir and the surface flow conditions are assumed to govern the hydraulic head distribution along the surface of the streambed; the feedback from the subsurface is, therefore, ignored. The left and right boundary conditions are assumed to have hydraulic heads corresponding to the total hydraulic head of the surface water at those locations, while there is no flow at the bottom boundary of the domain ( Figure  1). The Brinkman-Darcy equation for the simulation of subsurface flow is as follows [34]: where ρ is the density of the fluid (kg/m 3 ), u is the flow velocity (m/s), t is the time (s), p is the pressure (Pa), μ is the dynamic viscosity (Pa•s), μe is the effective viscosity (Pa•s), k is the intrinsic permeability (m 2 ), ε is the Brinkman porosity, and μu/k is the Darcy term (Pa/m).

Nitrogen Transport/Removal Calculation
For each subsurface flow path, the concentration of the nitrogen can be assumed to be related to the RT using the first order reaction expression [35]:

COMSOL Multiphysics for Hyporheic Flow Simulation
Imposing the hydraulic heads along the surface of the streambed as the upper boundary condition, the hyporheic flow is simulated using COMSOL Multiphysics. In this study, the weir and the surface flow conditions are assumed to govern the hydraulic head distribution along the surface of the streambed; the feedback from the subsurface is, therefore, ignored. The left and right boundary conditions are assumed to have hydraulic heads corresponding to the total hydraulic head of the surface water at those locations, while there is no flow at the bottom boundary of the domain (Figure 1). The Brinkman-Darcy equation for the simulation of subsurface flow is as follows [34]: where ρ is the density of the fluid (kg/m 3 ), u is the flow velocity (m/s), t is the time (s), p is the pressure (Pa), µ is the dynamic viscosity (Pa·s), µ e is the effective viscosity (Pa·s), k is the intrinsic permeability (m 2 ), ε is the Brinkman porosity, and µu/k is the Darcy term (Pa/m). For each subsurface flow path, the concentration of the nitrogen can be assumed to be related to the RT using the first order reaction expression [35]: where C 0 is the initial concentration of the nitrogen (i.e., 0.307 mol/m 3 in this study), C RT is the concentration of the nitrogen (mol/m 3 ) with the reaction time of RT, C is the change in the nitrogen concentration in the time duration of RT, and r is the reaction rate constant (day −1 ). It is worth noting that the unit of the RT in Equation (3) should be day. Thus, the nitrogen removal amount (NRA, mol) for each flow path can be expressed as follows: Then, the cumulative NRA (CNRA, mol) can be calculated as the sum of the NRA for all the flow paths within the maximum upstream distance in the subsurface flow influenced by the weir (L max ): (5) where N denotes the total number of the flow paths, and i denotes the i-th path line. Let NRR be the nitrogen removal ratio, which is regarded as another performance indicator. Here, the NRR of the HZ is defined as the difference between 1 and the ratio of the total nitrogen (TN, mol/m 3 ) in the HZ at the RTmean (i.e., the mean RT) and the total initial nitrogen in the HZ. Thus, the NRR can be computed as follows: Here, RTmean is used as the reaction time because it is considered to be the most single representative value given the variation of RT under different HZ scenarios [6].

Modeling Scenarios
To examine the impact of weir height on the HZ, various modeling scenarios are adopted, as listed in Table 1. In this study, the height of the weir varies from 0.3 m to 0.7 m with equal intervals of 0.1 m, and the variables to be held constant are defined as slope of 0.01, river discharge of 1 m 3 /s, and depth to bedrock of 5 m. It is worth noting that, in this study, the h values used for determining the optimal h lie within the range of the modeling scenarios (i.e., [0.3, 0.7] m), but with equal intervals of 0.01 m. For designing the optimal weir in hyporheic restoration, hyporheic flux and biochemical reactions (e.g., nitrification and denitrification) are two important factors. Previous studies have found that Water 2020, 12, 1399 6 of 16 higher in-stream structures can induce more hyporheic flux, and denitrification is closely related to the RT. Therefore, this study examines the optimal structure design to enhance denitrification in the HZ by increasing hyporheic flux and RT. Specifically, both the hyporheic flux and the RT impacted by the height of the weir are investigated.

Relevant Indicators
In order to make the results of this study more widely applicable and simplify the evaluation procedures of other studies, this study derives empirical equations based on the simulation results to describe the changes in the RT and hyporheic flux at the streambed with the weir in the stream.
This study assumes that the maximum upstream distance in the subsurface flow influenced by the weir (L max ) can be related to the height of the weir (h). Therefore, the L max can be expressed as follows: Based on the numerical modeling results, this study traces the subsurface flow paths of the particles that fall within the L max . It means that particle tracing is first used in the numerical simulations and then the RT values for the subsurface flow paths are used to develop the coefficients of the regression equations. This study assumes that (1) the RT and the hyporheic flux in the upstream (Q u ) can be expressed as functions of x using the regression method, where x is the distance measured from the left boundary of the domain, and (2) the coefficients of the regression equations can be related to h. Previous studies have shown that the regression equations in the form of power law can be the best to describe the RT and Q u [36,37] and, therefore, the RT and Q u can be expressed as follows: where a 1 (h) and b 1 (h) are the coefficients of the regression equation for the RT, and a 2 (h) and b 2 (h) are the coefficients of the regression equation for the Q u . Based on the simulation results under different scenarios with different h values, the forms of these coefficients can be determined (see Section 3 for details).

Objective Function
For each designated h value, the CNRA and NRR values can be calculated. It is worth noting that a higher h value can lead to a higher CNRA value, but a lower NRR value. As a result, there should be an optimal height of the weir which strikes a balance between the CNRA and NRR values. To this end, the objective function for the optimal design of the weir (i.e., the height of the weir) in this study is expressed by the product of the CNRA and NRR values for each designated h value. Furthermore, to facilitate comparison, the CNRA and NRR values are normalized using the following method: where CNRA j and NRR j are the CNRA and NRR values corresponding to the j-th h value, respectively; NCNRA j and NNRR j are the j-th normalized CNRA and NRR values; M is the data length, which is equal to the total number of h values. Therefore, the objective function is formulated as follows: Max. NCNRA j · NNRR j (12) It indicates that the optimal h of the weir can be obtained when the product of the corresponding NCNRA and NNRR values reaches the maximum.

Validation and Sensitivity Analysis Methods
Since the optimal h of the weir is obtained from the proposed method, it is important and necessary to validate the effectiveness of this method. In this study, two validation methods are used as follows: the first method compares the results from the numerical models for surface and subsurface flow simulation and the empirical equations for calculating the nitrogen removal (both amount and rate) in the HZ. The second method compares the results of denitrification simulation. The denitrification is simulated and is coupled with the hyporheic flow simulation using COMSOL Multiphysics. The following equation is used for denitrification simulation in the HZ: where u is the flow velocity in the streambed (m/s), D is the dispersion coefficient (m 2 ·s), c NO3-is the concentration of the nitrate (mol/m 3 ), and k d is the reaction rate constant of denitrification (s −1 ). The height of the weir is optimal if the nitrogen removal of the whole HZ is the maximum when compared to that of other heights. In this study, two indices are proposed to describe the nitrogen removal performance of the HZ driven by different heights of the weir. One is as follows: (14) where h i is the i-th height (m) and q HZ is the hyporheic flux (m 2 /s). This index considers the HZ as a single reaction tank, and the RT max is the maximum value of the denitrification reaction time. The other one is as follows: The second index regards the HZ as a continuous reaction tank for denitrification simulation. In this equation, T ≥ max(RT max (h i )).
For both of the above two methods, the result using the optimal h value will be obtained first, and then be compared against those obtained using other h values.
In addition, sensitivity analysis of the proposed method to optimize the height of the weir will be conducted from two aspects: the first investigates the sensitivity to the range of the h values used for searching the optimal h, while the second investigates the sensitivity to the variables to be held constant (e.g., river discharge, slope and depth to bedrock in Table 1). The weir is expected to influence the flow conditions of both surface and subsurface water, which can lead to the change of the L max ; moreover, a higher weir should lead to a larger L max . In this study, the simulation results show that the L max increases along with the increase of the height of the weir (h) (Figure 2). In order to evaluate the accuracy of the regression equation, the coefficient of determination (R 2 ) is selected as the assessment criterion. As the R 2 value is very high (i.e., 0.98), the relationship between L max and h can be expressed as follows:

Regression Equations
The positive coefficient of the h in Equation (16) confirms that a higher weir can lead to a larger L max . Moreover, the relationship is simple and linear, which could be partly due to the homogeneity of the streambed.  Figure 3 also proves that the Lmax increases along with the increase of the h. According to the assumptions in subsection 2.2.2, the RT can be expressed as the function of x in the form of a power law using the regression method. Therefore, the coefficients of the regression equations for all h values can be obtained, which are listed in Table 2. The R 2 values of the regression equations are also shown in this table, and they are all higher than 0.91, indicating that the regression equations fit the simulation results very well.    Figure 3 also proves that the L max increases along with the increase of the h. According to the assumptions in Section 2.2.2, the RT can be expressed as the function of x in the form of a power law using the regression method. Therefore, the coefficients of the regression equations for all h values can be obtained, which are listed in Table 2. The R 2 values of the regression equations are also shown in this table, and they are all higher than 0.91, indicating that the regression equations fit the simulation results very well.  Figure 3 also proves that the Lmax increases along with the increase of the h. According to the assumptions in subsection 2.2.2, the RT can be expressed as the function of x in the form of a power law using the regression method. Therefore, the coefficients of the regression equations for all h values can be obtained, which are listed in Table 2. The R 2 values of the regression equations are also shown in this table, and they are all higher than 0.91, indicating that the regression equations fit the simulation results very well.    As mentioned in Section 2.2.2, the coefficients of the regression equations can be related to the h, and Figure 4 shows that a 1 decreases with the increase of the h while b 1 increases with the increase of the h. Thus, the relationships between the two coefficients and h can be expressed as follows: Water 2020, 12, x FOR PEER REVIEW 9 of 16 As mentioned in subsection 2.2.2, the coefficients of the regression equations can be related to the h, and Figure 4 shows that a1 decreases with the increase of the h while b1 increases with the increase of the h. Thus, the relationships between the two coefficients and h can be expressed as follows:  It is worth noting that the logarithmic values of a1 (i.e., Ln(a1)) are used for establishing the regression equation because the range of a1 is large. The R 2 values are high, i.e., 0.89 for Ln(a1) and 0.88 for b1, respectively, which indicate that these two coefficients can be represented well by the h.
Additionally, in Figure 3, for different weir heights, the RT first generally increases along the upstream direction. It then sharply decreases and increases. This is possibly influenced by the backwater from the weir, and the location with sharp decrease and increase of the RT could show the extent of the HZ influenced by different ranges of the backwater from the weir. To be specific, when the weir height is 0.3 m, the influenced range is around 4 m (i.e., = 9-5 m in Figure 3) from the weir. However, when the weir height is 0.7 m, this influenced range could reach around 6.7 m (i.e., = 9-2.3 m in Figure 3). Figure 5 shows the relationship between the Qu and the distance from the left boundary of the domain under the five modeling scenarios with different h values. It can be observed that the maximum Qu increases along with the increase of the h. In order to better match the simulation results, Equation (9) is adjusted as follows:

Hyporheic Flux in the Upstream (Qu)
where the value of 9.5 simply accounts for the location of the weir measured from the left boundary of the domain. Table 3   It is worth noting that the logarithmic values of a 1 (i.e., Ln(a 1 )) are used for establishing the regression equation because the range of a 1 is large. The R 2 values are high, i.e., 0.89 for Ln(a 1 ) and 0.88 for b 1 , respectively, which indicate that these two coefficients can be represented well by the h.
Additionally, in Figure 3, for different weir heights, the RT first generally increases along the upstream direction. It then sharply decreases and increases. This is possibly influenced by the backwater from the weir, and the location with sharp decrease and increase of the RT could show the extent of the HZ influenced by different ranges of the backwater from the weir. To be specific, when the weir height is 0.3 m, the influenced range is around 4 m (i.e., = 9-5 m in Figure 3) from the weir. However, when the weir height is 0.7 m, this influenced range could reach around 6.7 m (i.e., = 9-2.3 m in Figure 3). Figure 5 shows the relationship between the Q u and the distance from the left boundary of the domain under the five modeling scenarios with different h values. It can be observed that the maximum Q u increases along with the increase of the h. In order to better match the simulation results, Equation (9) is adjusted as follows:

Hyporheic Flux in the Upstream (Q u )
where the value of 9.5 simply accounts for the location of the weir measured from the left boundary of the domain. Table 3   Similar to the RT, Figure 6 shows the relationships between the two coefficients of the regression equations for the Qu and h. It is observed that a2 increases with the increase of the h while b2 decreases with the increase of the h. Thus, the relationships between the two coefficients and h can be expressed as follows: It is worth noting that the R 2 value for a2 is equal to 1.00, which indicates the perfect fit; the R 2 value for b2 is also high, i.e., 0.87. As a result, these two coefficients can be closely related to the h.   Similar to the RT, Figure 6 shows the relationships between the two coefficients of the regression equations for the Q u and h. It is observed that a 2 increases with the increase of the h while b 2 decreases with the increase of the h. Thus, the relationships between the two coefficients and h can be expressed as follows:

Cumulative Nitrogen Removal Amount (CNRA) and Nitrogen Removal Ratio (NRR)
It is worth noting that the R 2 value for a 2 is equal to 1.00, which indicates the perfect fit; the R 2 value for b 2 is also high, i.e., 0.87. As a result, these two coefficients can be closely related to the h.  Similar to the RT, Figure 6 shows the relationships between the two coefficients of the regression equations for the Qu and h. It is observed that a2 increases with the increase of the h while b2 decreases with the increase of the h. Thus, the relationships between the two coefficients and h can be expressed as follows: It is worth noting that the R 2 value for a2 is equal to 1.00, which indicates the perfect fit; the R 2 value for b2 is also high, i.e., 0.87. As a result, these two coefficients can be closely related to the h.

Cumulative Nitrogen Removal Amount (CNRA) and Nitrogen Removal Ratio (NRR)
Based on the equations obtained in Section 3.1, the CNRA and NRR values for each designated h value can be calculated. For example, the RT and Q u values for each designated h value can be firstly calculated based on Equations (8) and (17)- (19), and then, the corresponding CNRA values can be calculated based on Equations (3)- (5). Similarly, the NRR values for each designated h value can be calculated based on Equations (6), (8) and (17). Figure 7 shows the relationships of the CNRA and NRR with the h. As h increases, the CNRA increases but the NRR decreases. It is clear that there will be an optimal h which strikes a balance between the two indicators.
Water 2020, 12, x FOR PEER REVIEW 11 of 16 Based on the equations obtained in subsection 3.1, the CNRA and NRR values for each designated h value can be calculated. For example, the RT and Qu values for each designated h value can be firstly calculated based on Equations (8) and (17)- (19), and then, the corresponding CNRA values can be calculated based on Equations (3)- (5). Similarly, the NRR values for each designated h value can be calculated based on Equations (6), (8) and (17). Figure 7 shows the relationships of the CNRA and NRR with the h. As h increases, the CNRA increases but the NRR decreases. It is clear that there will be an optimal h which strikes a balance between the two indicators.

Optimal Height of the Weir
The CNRA and NRR values are firstly normalized using Equations (10) and (11), and then, the objective function, which is the product of the CNRA and NRR values, is calculated. Figure 8 shows the relationship between the objective function and the h. It is observed that, as h increases, the objective function firstly increases and then decreases. The maximum value of the objective function is associated with the optimal h, which is 0.54 m.

Optimal Height of the Weir
The CNRA and NRR values are firstly normalized using Equations (10) and (11), and then, the objective function, which is the product of the CNRA and NRR values, is calculated. Figure 8 shows the relationship between the objective function and the h. It is observed that, as h increases, the objective function firstly increases and then decreases. The maximum value of the objective function is associated with the optimal h, which is 0.54 m.  (8) and (17)- (19), and then, the corresponding CNRA values can be calculated based on Equations (3)- (5). Similarly, the NRR values for each designated h value can be calculated based on Equations (6), (8) and (17). Figure 7 shows the relationships of the CNRA and NRR with the h. As h increases, the CNRA increases but the NRR decreases. It is clear that there will be an optimal h which strikes a balance between the two indicators.

Optimal Height of the Weir
The CNRA and NRR values are firstly normalized using Equations (10) and (11), and then, the objective function, which is the product of the CNRA and NRR values, is calculated. Figure 8 shows the relationship between the objective function and the h. It is observed that, as h increases, the objective function firstly increases and then decreases. The maximum value of the objective function is associated with the optimal h, which is 0.54 m.

Validation
In order to validate the effectiveness of this method, the numerical models of HEC-RAS and COMSOL Multiphysics are run using the optimal h value of 0.54 m, and the results are compared against those obtained with the h values of 0.5 m and 0.6 m. The CNRA and NRR values and their products under these three scenarios are calculated and listed in Table 4. Along with the increase of the h value, the CNRA value increases but the NRR value decreases. This is the same as that shown in Figure 7. The product under the scenario using the optimal h value of 0.54 m (i.e., 0.63) is larger than those under the scenarios using the h values of 0.5 m (i.e., 0.40) and 0.6 m (0.50). As a result, the good performance of the proposed method is confirmed. In addition, the effectiveness of the proposed method is further validated through directly simulating the denitrification in the HZ. The two indices for describing the nitrogen removal performance of the HZ are calculated, and Table 5 lists the results of denitrification simulations with the five h values. It is observed that both indices show that the height of 0.54 m gives the best performance. Therefore, the effectiveness of the proposed method is further validated from the aspect of temporal accumulation of the nitrogen removal in the HZ.

Sensitivity Analysis
In order to investigate the sensitivity of the proposed method to the range of the h values used for searching the optimal h, different ranges of the h values, e.g., [0.1, 0.8] m, [0.2, 0.8] m, and [0.2, 0.9] m, are selected. The relevant results are compared and shown in Figure 8. It is observed that the optimal h (i.e., 0.54 m, the dash line in Figure 8) remains the same, and that it is not sensitive to the range of the h values used for searching the optimal h.
The influence of different river discharges on the optimal height of the weir is also examined. As the river discharge of the base case is 1.0 m 3 /s (Table 1), two other river discharges (i.e., 0.5 m 3 /s and 1.5 m 3 /s) are used to investigate the change of the optimal h when river discharge increases/decreases by 50%. The denitrification simulation method mentioned in Section 2.3 is adopted and M1 den is employed to evaluate the performance of the weir in the nitrogen removal. For each of the three river discharges, simulations under three different heights of the weir (i.e., 0.53 m, 0.54 m and 0.55 m) are conducted in order to detect the generic trends of the optimal h. Table 6 lists the simulation results. If the river discharge decreases from 1.0 m 3 /s to 0.5 m 3 /s, the maximum M1 den value, which indicates the best performance in the nitrogen removal, is associated with the height of 0.53 m (see column 1 in Table 6), and the actual optimal height might be lower for the discharge of 0.5 m 3 /s. Even 0.53 m might not be the optimal h under the river discharge of 0.5 m 3 /s, and it is smaller than the optimal h under the river discharge of 1.0 m 3 /s (i.e., 0.54 m). Therefore, it can be concluded that the optimal h decreases when the river discharge decreases. Conversely, if the river discharge increases from 1.0 m 3 /s to 1.5 m 3 /s, the maximum M1 den value is associated with the height of 0.55 m (see column 3 in Table 6), which indicates that the optimal h increases when the river discharge increases. As a result, the optimal h can be generally sensitive to the river discharge. It is not surprising, as the larger river discharge may bring more nitrogen from surface water to the HZ; if the nitrogen loading is below the HZ nitrate removal capacity, the weir should be higher in order to improve the performance in the nitrogen removal because the higher weir can enhance the interaction between surface water and groundwater and the denitrification in the HZ. However, if the nitrogen loading is above the nitrogen removal capacity of the HZ, the water should be treated first instead of increasing the height of the weir. Finally, sensitivity analyses are also conducted under the case of the optimal height (i.e., 0.54 m) based on different slopes (i.e., 0.001, 0.005 and 0.05) and depths to bedrock (i.e., 3 m, 4 m, 6 m and 7 m), respectively, and hyporheic flux is used for evaluating the different HZ performances. The relevant results, i.e., the influences of different slopes and depths to bedrock on hyporheic flux in the upstream (Q u ), are shown in Figure 9. It is observed that hyporheic flux increases when the slope increases. In contrast, the influence of depth to bedrock on hyporheic flux is not significant, which indicates that the active zone of the HZ can be limited to a narrow range.

Limitations and Future Studies
The limitations of this study mainly include the following five aspects. First, the simulations in this study are all conducted in the steady state, ignoring the transient changes. Thus, the optimal in-stream structure design for purposes with temporal changes (e.g., under the transient conditions such as flooding, changes in carbon availability, microbiological populations, etc.) should be further investigated. Second, only the optimal height of the weir has been discussed in this study; following the proposed method, the optimal design of other variables of different in-stream structures (e.g., shapes, widths and thicknesses) can be further discussed. Third, the results are case-dependent. For other cases (e.g., cases of different slopes and depths to bedrock), the regression equations should be re-established based on the observed data or new simulation results. Fourth, the effects of regional groundwater inflow on the HZ are not considered in this study. However, according to our previous study [29], such effects could be ignored compared to the effects of the surface water. Finally, the results of this study have not been compared to the field data. Although it is in practice very challenging to obtain field data during river restoration projects, laboratory experiments may be conducted to further validate the proposed method.
With awareness of the above limitations, the method proposed in this study can provide a new avenue for investigating optimal in-stream structure design, which would be valuable for hyporheic restoration in river restoration projects. Hyporheic restoration can play its vital role in the ecological functions and services linked to riverine habitats through providing the proper hydro-environment and space. based on different slopes (i.e., 0.001, 0.005 and 0.05) and depths to bedrock (i.e., 3 m, 4 m, 6 m and 7 m), respectively, and hyporheic flux is used for evaluating the different HZ performances. The relevant results, i.e., the influences of different slopes and depths to bedrock on hyporheic flux in the upstream (Qu), are shown in Figure 9. It is observed that hyporheic flux increases when the slope increases. In contrast, the influence of depth to bedrock on hyporheic flux is not significant, which indicates that the active zone of the HZ can be limited to a narrow range.

Conclusions
This study proposes a method for optimizing in-stream structure design based on numerical modelling, through comprehensively considering nitrogen removal in the HZ. The major findings of this study can be summarized as follows: First, a numerical model which couples HEC-RAS and COMSOL Multiphysics by the hydraulic head along the surface of the streambed is built for both surface water and hyporheic flow simulations. Simulations under five different scenarios are conducted, and the model outputs are used to derive an optimal design of an in-stream structure (e.g., the weir in this study).
Second, based on the model outputs under five different scenarios, regression equations for estimating the relevant variables (e.g., the maximum upstream distance influenced by the weir, the RT, and hyporheic flux) are established. Then, the CNRA and NRR can be calculated.
Third, a framework of optimal in-stream structure design by considering the nitrogen removal in the HZ is proposed. Taking the height of the weir as the example, the objective function is formulated as the product of the normalized CNRA and NRR in this study. The optimal height of the weir can be obtained as the height at which the objective function is maximized. Validation shows that this method produces generally good performance. Sensitivity analysis indicates that the optimal height generally can be sensitive to the river discharge, i.e., the optimal height increases when the river discharge increases and vice versa. It is also observed that, in the case of the optimal height, hyporheic flux increases when the slope increases while the influence of depth to bedrock is not significant.
Overall, the modeling and results presented facilitate future studies on optimizing in-stream structure design for the nitrogen removal in the HZ, an important goal of river restoration projects.