Interpolation of water surface profile in unsteady open channel flow using the adjoint method based on 2 D shallow water equations

Real-time interpolation of high-water stage data obtained at separate water gauges was investigated to provide visual information for evacuation decision-making by residents near a river. The adjoint method was applied to a set of shallow water equations in a generalized coordinate system in which the hydraulic effect of tree communities is included as semipermeable blocks. The cost function to be minimized was assumed as the summation of four misfits for water level and discharge from observation and for Manning’s roughness (n) and tree community permeability (K) from standard values. The Lagrange multipliers method was adopted to determine the uncertain quantities, upstream discharge and spatial distributions of n and K, for minimizing the total cost under the constraint of the shallow water equations. The method was applied to a flood event observed in 2011 over a 20-km reach of the Tone River, where four water gauges are installed. The computation results successfully explained the time series of river discharge observed at stations located in the river reach, as well as the longitudinal profile of the trace water levels measured at 500-m intervals just after the flood event.


Introduction
The increasing frequency of intense rainfall due to recent climate change has increased the risk of levee overtopping.Since overtopping can quickly break down earth levees, residents near the river need real-time information on river flow conditions so they know when to evacuate.At present, however, the major flood information services in Japan handle only data from widely separated gauge stations, and residents often find it difficult to recognize threats.Provision of spatially denser information about the present river water stage in a visual form would be helpful for evacuation decision-making.
In this study, a method for real-time interpolation of water stage data obtained at separate water gauges was investigated to provide an accurate water surface profile in a visual form.The adjoint method [1][2][3] was applied to a set of shallow water equations in a generalized coordinate system in which the hydraulic effect of tree communities is included as semi-permeable blocks.Uncertain variables that needed to be optimized were the river discharge at the upstream end, Manning's bed surface roughness (n), and the permeability constant for immersed bulky vegetation (K).The total cost to be minimized was assumed as the summation of misfits not only for water levels and discharge from observation values but also for n and K from the standard values.
The method was tested for a flood event observed in 2011 over a 20-km reach of the Tone River, where four water gauges are installed at different locations.The river channel has a composite cross section with floodplains covered by grass, but some reed colonies and tree communities remain along the banks of the low-water channel.The channel topography, including tree community distribution with height classification data, was provided by the Tone River administration office.Computation results were compared with observed data for longitudinal water surface profile and river discharge, which were provided by the same office.

Model description
The cost function is expressed as follows: where J T is the total cost to be minimized and is a function of three uncertain control variables.Variable Q bound is the discharge at the upstream end, n is Manning's bed roughness, and K is the permeability constant for tree communities.
Terms J H and J Q on the right side of the equation are the total misfit for the water level and discharge, respectively, and J n_int and J K_ini are the sum of the squared differences between the values obtained from the assimilation analysis and the initially assumed (standard) values of n and K, respectively.The values for n and K were constrained in some ranges to avoid any locally abnormal values.Among the four terms, the most influential one is the first regarding errors in water level.Meanings of the other parameters are as follows: T a and T are initial and last time of calculation, respectively; (x mH , y mH ) are coordinates for the location of the water gauge; t k is the observation time; x mQ is the location of the discharge observation station; H 0 and Q 0 are the observed water level and discharge, respectively; H c and Q c are calculated water level and discharge, respectively;  H0 and  Q0 are the standard deviation of the error of observed water level and discharge, respectively; A is the area of the horizontal plane; num_an and num_ak are numbers representing the classifications of Manning's bed surface roughness and the permeability constant, respectively; n init and K init are the initial values of n and K, respectively, in the assimilation analysis; and  n_init and  K_init are the standard deviations of n and K, respectively.
The continuity equation is written as F h (U, V, h) = 0, where U and V are longitudinal and transverse velocity, respectively, and h is water depth.Then, the longitudinal and transverse momentum equations are written as F U (U, V, h) = 0 and F V (U, V, h) = 0, respectively.Finally, the functional to be minimized in the Lagrange multiplier method under the constraint of the shallow water equations is written as follows:  dAdt The adjoint matrix that Lagrange multipliers satisfy was obtained by taking variation of U, V, and h by Eq. ( 2).The gradient vectors for control variables obtained by differentiating Eq. ( 2) was solved using a quasi-Newton method.The iteration process is summarized in Fig. 1 and listed below as steps: (a) The unsteady shallow water equations are solved with initial values of Q bound , n, and K to obtain U(x, y; t+Δt), V(x, y; t+Δt), and h(x, y; t+Δt) over the entire calculation area for the next calculation time step.(b) Misfits are evaluated using Eq. ( 1).(c) Lagrange multipliers are calculated from the adjoint matrix.(d) The functional L is calculated from Eq. ( 2) and gradient vectors of L relative to Q bound , n, and K are estimated.If convergence is achieved, the process ends.(e) The variables Q bound , n, and K are updated using a quasi-Newton method if L did not converge.Then, the process returns to step (a) for the next iteration.Standard values of Manning's roughness n for each ground cover classification were determined based on Fukuoka's report [4].For permeability constants of tree communities K, we adopted the values recommended by Fukuoka et al. [5].Those parameter values are listed in Table 1.The time increment for unsteady flow calculation was set to 0.5 s in real time.

Results and discussion
Figure 5(a) shows the variation of the cost function relative to the initial value and iteration steps for the computation at T = 6:00 a.m. on September 22 near the flood peak.It decreased sharply through the first 5 iteration steps and became almost constant after the 10th iteration step.Figure 5(b) shows the convergence process of n for each ground cover classification and K.All values became stable by the tenth iteration step.Calculation of 10 iterations for one time increment (0.5 s) took only 0.07 s of CPU time using an IBM POWER8 (10 cores × 2 sockets).Therefore, this data processing method can provide the flow field within a 20km reach fast enough for information alerts of the river condition in a visual form.Colored curves in Fig. 6 show the calculation results for river discharge across gauge station 1 for three cases in which the number of iterations was 1, 5, and 11.The black dotted line shows the reference values estimated from water gauge data using the H-Q relation.With increasing iterations, the calculated water stage became higher than the reference values during the flood peak and lower during the receding phase.This is because the estimation from the H-Q curve was not accurate due to unsteady flow; the water surface slope was steeper in the rising phase than that in the receding phase.
Figures 7(a) and 7(b) compare the calculated and observed water surface levels at gauge stations 1 and 3, respectively, in the same manner as Fig. 6.The calculated result for station 1 was remarkably improved by increasing the iterations, while the results for station 3 agreed with the observed levels even for one iteration because the water level observed at station 4 was given as a boundary condition at the downstream end of the calculation area.The curves in Fig. 8 show the longitudinal water surface profiles averaged for cross sections during the rising phase (green), receding phase (black), and peak for the whole calculation (blue).Colored triangles indicate the observed values at water gauge stations corresponding to the three lines.The red and blue dots are the trace water levels measured at 500-m intervals just after the flood event.The water surface slope in the rising phase was much steeper than in the receding phase, and the calculation showed good agreement with observation for both the rising and receding stages.The curve for the highest water level during the flood generally agrees with the trace water levels.Figure 9 displays a snapshot of the velocity field near water gauge station 1 obtained by our method.Because a flow field picture like this can be obtained for the whole study area at each calculation time step, the output from this method can be presented as visual information helpful to residents near the river through the internet, for example, so they can recognize threats.

Conclusion
The adjoint method was applied to a set of shallow water equations to interpolate the water surface profile from data obtained from water gauges at different locations.The uncertain quantities to be optimized were river discharge and flow resistance factors.The major conclusions obtained from its application to a flood event observed over a 20-km reach of the Tone River are summarized as follows: (1) The computation results successfully replicated the hydrographs observed at the four gauge stations installed in the river reach.In addition, the longitudinal profile of highest water level showed good agreement with the trace water levels measured at 500-m intervals just after the flood event.
(2) Values of uncertain quantities converged to the final values within 10 iterations for the minimization of the total cost function.The CPU time for a 10-iteration calculation was approximately one-eighth of real time, which was fast enough to provide real-time information on the river condition.(3) Because this method provides two-dimensional flow field data as well as the longitudinal water surface profile, the calculation output can be presented as visual information through the internet, which would help residents near the river recognize flooding threats.

Figure 2
Figure 2 shows the location of the study field.The Tone River has its source in the central mountains of Honshu Island, Japan, and flows across the alluvial plain on the north side of Tokyo, eventually discharging into the Pacific Ocean.The model was applied to the flood event observed in September 2011 over a 20-km reach from 152 km to 133 km, where water gauges are installed at four sections.

Fig. 2 .
Fig. 2. Location of study field in the Tone River watershed

Fig. 3 .
Fig. 3. Topography and ground cover of study field

Figure 4 (
Figure 4(a)  shows the flood hydrographs obtained at the four gauge stations during the 2011 flood.A time lag of approximately 3 h was observed between upstream station 1 (150.5 km) and downstream station 4 (133.0km).Figure4(b) shows the river discharge across station 1 estimated from the head-flow (H-Q) curve calibrated with flow data from previous floods.The peak discharge was approximately 4,000 m 3 /s, which is about one-fourth of the design flood discharge at that station.

Fig. 4 .
Fig. 4. Time series of water levels and river discharge

Fig. 6 .
Fig. 6.Dependence of calculated river discharge on iteration number

Fig. 8 .
Fig. 8. Longitudinal water surface profiles in rising and receding phases and highest water level profile

Table 1 .
Standard values of Manning's roughness and permeability constant