A Direct Solution Approach to the Inverse Shallow-Water Problem

The study of open channel flowmodelling often requires an accurate representation of the channel bed topography to accurately predict the flow hydrodynamics. Experimental techniques are the most widely used approaches to measure the bed topographic elevation of open channels. However, they are usually cost and time consuming. Free surface measurement is, on the other hand, relatively easy to obtain using airborne photographic techniques. We present in this work an easy to implement and fast to solve numerical technique to identify the underlying bedrock topography from given free surface elevation data in shallow open channel flows. The main underlying idea is to derive explicit partial differential equations which govern this inverse reconstruction problem. The technique described here is a “one-shot technique” in the sense that the solution of the partial differential equation provides the solution to the inverse problem directly. The idea is tested on a set of artificial data obtained by first solving the forward problem governed by the shallow-water equations. Numerical results show that the channel bed topographic elevation can be reconstructed with a level of accuracy less than 3%. The method is also shown to be robust when noise is present in the input data.


Introduction
The study of open channel flows requires an understanding of the hydrodynamics of the flow in order to accurately capture its main characteristics.Designing and studying engineering structures on rivers, understanding the hydraulic phenomena of the watercourse for water quality control, and prediction of exceptional natural events associated with water flow constitute the three main objectives of the study of open channel flow simulation, 1 .All modelling studies of such applications require well-established governing equations along with the predefined parameters, boundary conditions, and initial conditions.Previous studies show that there has been tremendous progress in the field of hydraulic modelling of open channel flows.However, the representation of accurate channel bed topography is still a challenge, 2-4 .Parameters such as bedrock topographic elevation and roughness coefficients are required prior to the simulation of open channel flows.The level of discrepancy in these data will have a direct effect on the quality of the simulation results.Thus in order to accurately predict the hydrodynamics of open channel flows, precise bedrock topographic elevation data and other hydraulic parameters are required.In 5 and references therein, it is argued that some hydrological models contain conceptual parameters that usually cannot be accessed by in-situ data measuring techniques.On the contrary, some models only use measurable parameters from field surveys.But the procedure is often long and expensive because it is difficult to determine parametric values at each computational grid points.
There have been significant developments in experimental techniques to measure river bathymetry and flow depths.Examples of these techniques for river bathymetry include interferometric synthetic aperture radar SAR digital photogrammetry 4, 6 used for the North Ashburton River, New Zealand; airborne laser altimetry LiDAR ; ground survey 2 used for the Yakima and Trinity River Basins in the USA to assess its ability to map bathymetry.However, most experimental techniques to identify the bed rock elevation are time consuming 4, 7 and expensive 2, 5 especially for wide, gravel bed, and long reach rivers.
On the other hand, many numerical techniques have been developed to simulate open channel flows.The characteristics of the numerical methods differ depending on the particular problem it addresses: methods suitable for studying hydraulic jumps; methods suitable to study dam break problems; methods to study flood plain flows; methods to study steady flows and methods suitable for high Froude numbers are relevant scenario worthy of mention.A detailed review of the relevant numerical techniques on the forward problem is beyond the scope of this study.
Comparatively, very little has been studied on the numerical reconstruction of the bedrock elevation data of open channel flows.Previous studies have shown that there are mainly two approaches that can be implemented to solve such an inverse problem: the direct numerical approach and optimization-based iterative methods.The former method was previously implemented by the authors of this work, 8-10 to infer the bedrock topography from a known free surface data in one-dimensional shallow water flows and by 11-13 to infer the substrate topography and the corresponding flow field 14 from free surface data in thin film flows.This paper extends to three-dimensional shallow-water flows the methodology presented in 8-10 which only dealt with planar flows only dealt with planar flows.In 15 , a numerical bedrock reconstruction approach from known free surface elevation data is presented for zero-inertia two-dimensional shallow-water flows, while the optimization-based reconstruction approach has been implemented by many authors who tried to minimize a cost functional iteratively, see 16, 17 .In the following the forward problem refers to the simulation of open channel flows for the prediction of the water surface elevation, the velocity field, and flood coverage from known hydraulic parameters such as bedrock elevation distribution and the roughness coefficient.Conversely, we call the inverse problem when the bedrock elevation data is inferred from the given free surface elevation data and other hydraulic parameters.Section 2 introduces the shallow water equations and the respective discretization technique for the forward problem.Section 3 presents a set of forward problem benchmark test cases.Following Section 3, the inverse problem governing equations and the corresponding discretization technique and the respective inverse problem test cases are presented in Section 4. In Section 5, we present the concluding remarks.

Governing Equations
Shallow-water flows are three dimensional in nature, but the three-dimensional Navier-Stokes equation can be simplified to a set of two-dimensional equations by considering vertically averaged quantities.This is done by integrating the Navier-Stokes equations over the flow depth for an incompressible fluid.While depth averaging, the two-dimensional shallow-water equations are derived with the following assumptions.a the pressure distribution is hydrostatic over the flow depth, b the angle of inclination of the channel is small such that flow depths measured along the vertical and normal to the channel bed directions are approximately the same, c the length in the vertical direction of the flow is smaller than that in the main flow direction, and d the Coriolis and wind stress effects are neglected, 18 .The hydrostatic pressure distribution assumption is valid for the case of long and shallow waves, wave length is much larger than the depth of the flow , in which the vertical acceleration of fluid elements during the wave passage stays small.The solution of shallow water equations may not accurately represent open channel flows with short or high waves, because the hydrostatic pressure assumption will be violated for such flows, 18 .A general two-dimensional flow orientation is shown in Figure 1.
Thus, the two-dimensional governing equations can be written in conservative form 18 as where u x, y, t is the velocity component in x direction; v x, y, t is the velocity component in y direction; h x, y, t is the depth of the flow; C o is dimensional constant: C o 1 for SI unit system; n is the Manning's roughness coefficient; z x, y is bed topography; g is the acceleration due to gravity.The free surface function, ψ x, y, t h x, y, t z x, y , is the measure of the water surface elevation, and the corresponding depth of the flow is the difference between the water surface and the channel bed elevations.

Discretization of the Governing Equations of the Forward Problem
The discretization of the governing equations is one of the most important steps in the numerical solution of partial differential equations in engineering and science.Many discretization techniques have been developed to suit specific problems.In the context of this work, the MacCormack explicit numerical scheme is implemented for the forward problem.It has been chosen because it is a standard technique which is reasonably easy to implement and which produces reliable results.This discretization technique has a predictor and a corrector stage coupled with sequential solution methodology.Backward differencing numerical technique is implemented in the predictor stage, while forward differencing is used in the corrector stage.The flow variables are known at time step k, and their values at k 1 are to be determined, 18 .Then, for grid point i, j , the approximate equation can be written as the following.

2.4
As mentioned above the scheme first uses backward space differencing ∇ x and ∇ y to predict an intermediate solution from the known information at time step k.Then the forward differences Δ x and Δ y are used in the corrector stage to correct the predicted values.The finite difference operators ∇ and Δ are defined as The subscripts in the operators show the direction of differencing.Substituting the flux terms as U uh and V vh, the final form of the discretized equations can be written as the following.
Predictor Stage: gh k i,j Δt S oy − S fy i,j .

Mathematical Problems in Engineering
The new values at time step k 1 are then calculated from the intermediate values which are determined from the predictor and corrector steps: whereas the primitive variables will be determined from the computed values of U and V in each time step: Reflective wall boundary is implemented at the side walls.This can be implemented by replacing the fictitious point in the solid wall by its mirror point in the flow domain while changing the sign of the normal component of velocity.
However, the inlet and the outlet boundary conditions are handled depending on the type of the flow.For subcritical flow the depth of the flow at the outlet boundary condition is defined with the flow rate at the upstream boundary.On the contrary, the depth and the flow rate at the upstream boundary are defined for supercritical flows.However, supercritical and transcritical flow test cases are not included in this study because of the limitation of the proposed approach.On the other side of the defined boundary conditions, it is assumed that the gradient of the parameters along the flow direction is zero.The transverse velocity component at the inlet is defined also as zero.

The Forward Problem
The governing equations of the forward problem along with its discretization are presented in the above section.The parameters that will be determined from the forward problem are the free surface elevation, the velocity components of the flow, and the depth of the flow.For this purpose some hydraulic parameters are required prior to the computation.These include bed topography elevation, roughness coefficient, flow rate, and the depth of the flow at the boundary depending on the flow regime.In the following, the results of the forward problem test cases are presented, and these results will be used as input data for the validation and assessment of the inverse problem methodology.

Test Case I: One-Dimensional Steady Subcritical Flow over a Frictionless Channel
A 1 m wide 25 m long channel is considered to test subcritical flow over a bump.The channel bed is assumed to have a rectangular cross section and frictionless with a bump.The bed topography is defined by  The depth at downstream boundary condition of the flow is set as h 2 m, and the water inflow condition Q 4.42 m 3 /s is imposed for the computation of results.This test case is the most common canonical test case used by different authors, 19-21 for example, in order to validate their numerical techniques for the convergence of the solution towards steady state and the conservation of discharge along the channel.The Froude numbers for this test case range from 0.496-0.635showing that the flow is subcritical in the entire domain.
The results are in a good agreement with the results presented in the above references and with the respective analytical solution presented in 22 .As can be seen in Figures 2  and 3, the existence of the obstacle creates significant change on the free surface and depth profiles.

Test Case II: Three-Dimensional Flow over a Three-Dimensional Bump
In this test case, an artificial bed topography with a bump having three-dimensional characteristics is considered to simulate three-dimensional shallow water flows.Thus variation of the bed along the longitudinal and the transverse direction is considered in a computational domain of size 100 m by 100 m.The bump is defined by where x max and y max are the length scales of the domain in x and y directions.In this test case x max y max 100 m and α 1 6 α 2 10 are constants.A spatial grid size Δx 1 m and the temporal grid size Δt 0.1 s are implemented.The flow rate Q 2 m 3 /s per unit width is defined at the inlet boundary.A Manning's roughness value n 0.02 is chosen to account for the effect friction.A subcritical flow is considered such that at the outlet a depth is 1.0 m is imposed along the transverse direction.Steady solution is generated from the transient governing equation and steady boundary conditions.The Froude numbers for this test case range from 0.56-0.92showing that the flow is subcritical in the entire domain.In the following, plots of the bedrock elevation along the centre lines and the respective depth of the flow for forward problem are presented.
From the given bedrock topographic elevation, steady flow rate, roughness coefficient, and boundary conditions, the steady flow results are generated.The free surface variation along the centreline is given in Figure 4.As can be seen from Figures 4 and 5, the existence of the bump at the centre of flow domain creates a three-dimensional variation in the free surface elevation.The underlying hypothesis of this work is that the shape of this free surface with known boundary conditions contains sufficient information to reconstruct the underlying bedrock.This will be discussed in detail in the inverse problem analysis section.The free surface elevation data of this test case will be used as a known parameter along with other boundary conditions for the inverse problem analysis.

Test Case III: Three-Dimensional Flow over Number of Three-Dimensional Bumps
In order to investigate the effect of a more complex bedrock topography on the free surface, four bumps are considered in the computational domain.In this test case α 1 6 α 2 10 values in the bed topography function are considered.The flow rate Q 4.42 m 3 /s per unit width defined at the inlet boundary along with a depth of 2.0 m at the outlet boundary is imposed along the transverse direction.Grid sizes and Manning's roughness coefficient values considered here are similar to that of test case II.The free surface elevation is generated for the inverse problem analysis by solving numerically 2.6 , 2.7 , 2.8 , and 2.9 .The bedrock topography is given by  From Figures 6 and 7, it can be seen that each bump introduced on the bed has its own effect on the free surface.The undulations on the free surface are the manifestations of the bump on the bedrock, and these free surface undulations, if accurately captured, can provide sufficient information of the bedrock and other parameters.The Froude numbers for this test case range from 0.49-0.6showing that the flow is subcritical in the entire domain.The free surface generated in this test case will also be used as input data for the respective inverse problem test case.

The Inverse Problem
Unlike the forward problem the inverse problem analysis will be dealing with identifying hydraulic parameters from the known fields and parameters.The shallow-water equations will be rearranged and used for the inverse problem analysis.However, there is a difference between the known and unknown parameters.The governing equation of the inverse problem can be rewritten after substituting a free surface function by ψ x, y z x, y h x, y where ψ x, y is the free surface elevation function.
The discretization of the inverse problem governing equations differs from the forward problem in such a way that only backward differencing explicit numerical scheme is implemented.This approach has stability restriction but is easy to implement and fast to provide a solution: Δt S ψy − S fy i,j .

4.2
The boundary condition is implemented in a similar way as that of the forward problem.The new values at time step k 1 are calculated from the previous time step values.In this analysis, the expression used to find the new values differs from the forward problem approach.This is because, in the inverse problem, unlike to the forward problem analysis, no information is transferred from downstream to upstream direction making the backward difference approach a suitable choice.As the backward differencing scheme is implemented, the definition of the outlet boundary condition does not affect the solution of the problem in the upwind direction as there is no upstream wave propagation.Thus, inlet boundary condition is implemented.The primitive variables will be determined from the computed values of U and V in each time step.Once the values of the primitive variables are determined the bed topography elevation will be evaluated by a simple subtraction for the depth of the flow from the free surface elevation function: The results of the inverse problem will be counterchecked with the forward problem results in order to identify the capability of the algorithm in reconstructing the bed topographic elevation.Thus similar test cases will be considered for the numerical experiment of the algorithm.

Test Case I: Steady One-Dimensional Subcritical Flow over a Bump in a Frictionless Channel
The free surface elevation obtained from the forward problem analysis is used as an input parameter for the inverse problem.This includes the free surface profile and the boundary and initial conditions along with the steady inflow rate.Unlike the depth at the downstream boundary in the case of the forward problem, the inverse problem requires boundary conditions at the inlet boundary.A depth at the upstream boundary h 2 m and the water inflow rate Q 4.42 m 3 /s are imposed.A time step of 0.05 sec and spatial grid size of 0.5 m are used in the computation.Figure 8 shows the comparison of the reconstructed and the actual bed forms.In Figure 9 the reconstructed and the actual depth of the flow are compared.The results show that there is a perfect agreement between the reconstructed and the forward problem parameters.Relative to the results presented in 16 , this approach has a capability to overcome the constant shift of the reconstructed bed topography.This confirms that the algorithm can be used for bedrock reconstructions of open channel flows of one-dimensional nature.

Test Case II: Steady Three-Dimensional Flows over a Bump in a Rectangular Channel
This test case is used to test the algorithm for its capability to reconstruct the unknown bed topographic elevation from known free surface elevation data which was generated in the forward test case II.In the respective forward test case, the effect of bed slope, friction coefficient, and local bed elevation was implemented to generate the free surface profile which we use as input parameter.Along with the free surface profile, depth variation and flow rate at the inlet boundary are used.For the inverse problem analysis, a spatial grid size Δx 1 m and the time step Δt 0.05 sec are implemented.The flow rate of Q 2 m 3 /s per unit width and depth at the upstream boundary which is imported from the forward problem are used.A Manning's roughness value n 0.02 is used as in the case of the forward problem.A total computation time of 100 seconds was sufficient to generate steady reconstructed parameters.
As can be seen from Figure 10, the bedrock is successfully reconstructed with a very good agreement with the actual bed used in the forward problem.The bump is well reconstructed along with the bed slope showing that the methodology used can use the signal from the free surface to successfully identify its cause.Quantitatively the bed topography is reconstructed with 3% maximum deviation.The sensitivity of the algorithm to the introduction of noise in the free surface data was also tested on ranges of noise magnitude with respect to the signal on the free surface.Noise magnitudes ranging from 1-5% of the "peak-to-peak" free surface deviations are considered.In Figure 11, a reconstructed bed topography from a 1% noise is shown, and it is evident that the noise is not amplified in the reconstructed bed.

Test Case III: The Reconstruction of Bed Topography Used in Test Case III of the Forward Problem
Like the above inverse problem test cases, the free surface topographic data is used as input parameter to reconstruct the corresponding bed topographic elevation.Additionally a flow  rate of Q 4.42 m 3 /s per unit width and depth h 2 m at the upstream boundary are imported from the corresponding forward problem.The reconstructed bed along with the free surface profile is shown in Figure 12.
As can be seen from Figure 12, the four bumps in the original bump are successfully reconstructed.However, there is an apparent but small difference, downstream of the bumps,  between the reconstructed and the original bed.The comparison of the bed topography along the centreline is presented in Figure 13.
As can be seen in Figure 13, the reconstructed and original bed forms are in good agreement with each other.The difference between the original and the reconstructed bed lies below 0.5%.
Like test case II, the reconstruction algorithm is tested for its applicability to infer the channel bed topography from noisy free surface elevation data.The magnitude of the noise introduced is 1% of the signal of the free surface variation, and the results Figure 14 show that the bed topography is reconstructed without the amplification of the noise in the reconstructed bed elevation.
The above test cases are all based on numerical results with smooth free surfaces which are used as input variable.However, in practice it is difficult to get a smooth free surface from field measurements.Often the measured free surface elevation includes noise.Thus, to identify the capability of the algorithm to handle noisy input data, we introduce 1, 2.5, and 5% noise in the free surface data of all of the previous test cases.The noise is based on the signal on the free surface.A Savitzky-Golay smoothing filter is used in MATLAB to smooth the 2D noisy free surface data before the computation.The reconstruction of the bed topographies from the noisy data revealed that the level of discrepancy in the reconstructed bedrock elevation is in the range of 1%, 5%, and 9.5%, respectively.
The algorithm used for the inverse reconstruction of the bedrock topography from the free surface elevation data is fast and easy to implement.It has the ability to reconstruct the bed from noisy free surface data with the help of smoothing techniques.The methodology requires the values of the steady flow rate, depth of the flow at inlet boundary, and the roughness coefficient in addition to the free surface data.The flow in all test cases is subcritical as it is shown by the calculated ranges of the Froude number.Similar analysis indicated that the proposed approach is not suitable for transcritical and supercritical flows because hydraulic jumps and surges are not accommodated by the solution approach; thus it requires further developments to perform bedrock reconstruction in the cases of such flows.
The presented methodology works well with steady flows of subcritical nature.The effect of a flux change in the flow will have local effect on the free surface which will send a wrong signal to the algorithm in the reconstruction process.Thus scope of this study is therefore limited to steady shallow open channel flows.

Conclusion
The study of open channel flow modelling calls for input parameters like bed topography and roughness coefficient in order to accurately predict the hydrodynamics of the flow.A methodology based on an explicit finite difference scheme is used to reconstruct the bedrock elevation from the free surface data.The methodology requires a steady flow, the knowledge of the roughness coefficient, and the depth of the flow at the inlet boundary.
The algorithm is tested on a set of benchmark test cases, and encouraging results are found.The bed topography is well reconstructed with a deviation of 3% or less.The numerical methodology is easy to implement and fast to produce a solution but has a CFL restriction because of its explicit nature.This approach is suitable for steady open channel flows for which the shallow-water approximation holds where the signal due to the existence of the underlying bed topography is captured by the free surface measurement technique.
In practice, the measured free surface contains noise.The presented methodology is found to be capable of reconstructing the channel bed topography from noisy data.When tested on a set of noisy numerical data, the methodology was found to introduce no noise amplification in all test cases considered.

Figure 3 :
Figure 3: Depth of the flow.

Figure 4 :
Figure 4: Centreline bed and free surface elevation variation a along the length and b along the transverse direction.

Figure 6 :
Figure 6: Water level on an inclined bed with four bumps.

Figure 7 :
Figure 7: Bed and free surface along the centreline a longitudinal b transverse variation.

Figure 9 :
Figure 9: Comparison of the reconstructed and forward problem depth variation for the case of subcritical flow.

Figure 10 :
Figure 10: Comparison of original and reconstructed bed forms.

Figure 11 :
Figure 11: Reconstructed bed form from noisy free surface data.

Figure 13 :
Figure 13: Comparison of the reconstructed and original bed forms.

Figure 14 :
Figure 14: Comparison of the original and bed reconstructed forms from noisy free surface data.