A Local Multi-Layer Approach to Modelling Interactions between Shallow Water Flows and Obstructions

The capability to accurately predict flood flows via numerical simulations is a key component of contemporary flood risk management practice. However, modern flood models lack the capacity to accurately model flow interactions with linear features, or hydraulic structures like bridges and gates, which act as partial barriers to flow. Presented within this paper is a new Riemann solver which represents a novel approach to modelling fluid-structure interactions within two-dimensional hydrodynamic models. The solution procedure models obstacles as existing at the interface between neighbouring cells and uses a combination of internal boundary conditions, different forms of the conservation laws and vertical discretisation of the neighbouring cells to resolve numerical fluxes across a partially obstructed interface. The predictive capacity of the solver has been validated through comparisons with experimental data collected from experiments conducted in a state-of-the-art hydraulic flume. Since the solution procedure is local, only applying to the cells within the immediate vicinity of a structure, the method is designed to be compatible with existing two-dimensional hydrodynamic models which use a finite volume scheme to solve the shallow water equations.


Introduction
The ominous threat of anthropogenic climate change is driving the requirement for more effective flood risk management in order to better manage what is already a challenging and costly hazard; models estimate that forecasted average annual flood losses for the United States will increase from US$32 billion to more than US$40 billion by 2050 [1], with similar predictions of increasing flood risk being made on a global scale [2].Hydrodynamic models play a vital role in contemporary flood risk management by providing evidence, via numerical predictions, upon which the quantification of flood risk and consequential future investment is based.It is therefore vital for effective flood risk management that hydrodynamic models produce accurate predictions.
Within catchments, channel structures, such as bridges, weirs and gates, can act as obstacles to flow, significantly influencing the local flow characteristics [3].However, within modern hydrodynamic modelling practice, methods for modelling such features are relatively under-developed, with industry standard models using coarse approximations, empirically based methods or even omitting such features entirely [4,5,6].Within academic literature there have been a number of contributions towards bridging this gap in modelling capacity, such as [7,8,9,10,11], however, none of the published works present an accurate method for the generalised treatment of partial barriers to flow within two-dimensional hydrodynamic models.
Within Mckenna et al. [12], the authors of this paper presented a new Riemann solver capable of resolving numerical fluxes across a partially obstructed interface.The proposed solution procedure represents structures as existing at the interface between neighbouring cells and uses a combination of internal boundary conditions and a different form of the conservation laws in the adjacent cells, to resolve numerical fluxes across the partially obstructed interface.Experimental validation, via experiments conducted in a state of the art research flume, demonstrated the accuracy of the solver for a range of flow conditions and barrier configurations.
Despite the successful validation of the solver, there is opportunity for enhancement of the method via more accurate discretisation of the horizontal velocity in the vertical plane.As such, this paper aims to use the basic conceptual idea underpinning the Riemann solver developed in [12], which is the decomposition of the Riemann problem in the vertical plane, to develop a new, more sophisticated and accurate method for representing structures within two-dimensional hydrodynamic models.As for the development of the previous solver, compatibility of the method with existing flood models utilising two dimensional finite volume schemes to solve the shallow water equations was a key consideration throughout the development of the solver.

Mathematical Model
The proposed solution method divides the computational domain into structure cells, intermediate cells and normal cells with corresponding normal interfaces (NI), intermediate interfaces (II) and structure interfaces (SI) as shown in Figure 1.
Intermediate Cell Normal Cell Intermediate Cell Normal Cell At a structure interface, the adjacent structure cells are vertically discretised into sub-cells with a maximum depth capacity corresponding to the dimensions of the idealised structure represented at the interface as shown in Figure 2.For example, the sub-cells U i,1 and U i+1,1 in Figure 2 have a maximum depth capacity of h 1 = z 3 2 − z b , which represents the difference in elevation between the base of the structure and the bed.
For normal interfaces and the corresponding adjacent normal or intermediate cells, a one-dimensional (1D) Figure 2: Division of structure cells into sub-cells corresponding to the base and cover of the idealised structure modelled at the structure interface.The maximum depth capacity of flow in the layer one is z 3/2 − z 1/2 and the maximum depth capacity of flow in the second layer is equal to z 5/2 − z 3/2 .The uppermost layer has no maximum depth capacity.
FV scheme is used to solve the 1D Shallow Water Equations (1D-SWE) given as: Where U is the vector of conserved variables, F(U) is the vector of fluxes and S(U) is a vector of sources comprising of S 0 , the bed slope source term and S f , the bed friction source term.These terms are given as follows: Whereby h denotes the depth of flow, u denotes the velocity component in the x direction, g is the acceleration due to gravity, z is the elevation of the bed and τ f is the shear stress due to bed friction in accordance with Manning's equation: Where n is Manning's roughness coefficient.
For the structure and intermediate interfaces and corresponding adjacent structure and intermediate cells, a 1D FV scheme is used to solve a multi-layer 1D shallow water system [13]: Where U k is the vector of conserved variables for the layer k, F(U k ) is the vector of fluxes for layer k and S k (U k ) is a vector of sources for layer k comprising of S k,0 , the topographic source terms for layer k and S k,f , the friction source terms for layer k.These terms are given as follows: Where k refers to the index of the layer under consideration, labelled in ascending order from layer 1 at the bed, to layer n at the free surface.k + 1/2 and k − 1/2 refer respectively to the upper and lower interface for layer k.The subscript k (+) refers to the properties of the flow above layer k and the subscript k (−) refers to the properties of the flow below layer k, which are defined respectively as: R k+1/2 and R k−1/2 refer to the reaction forces exerted at the interfaces between the layers, with R k+1/2 denoting the reaction force of layer k onto the fluid above and R k−1/2 denoting the reaction force exerted on layer k by the fluid or bed beneath it.τ k+1/2 and τ k−1/2 represent the interlayer viscous friction effect induced at the upper and lower interfaces of layer k.The interlayer friction terms are derived for a multilayer cell by applying a finite difference approximation, across the depth of the fluid layer k, to the viscous stress component of the incompressible Navier-Stokes system, as proposed by Audusse et al., [14]: For the case where k = 1, considering the layer which flows over the bed, τ k−1/2 = τ 0 which is instead derived from Manning's equation (3), where H is the total depth of flow for the whole structure cell.The particular form of the viscous effect on the base of the fluid layer, τ k−1/2 , is accounted for by Kronecker delta in (8), which is defined as: The Kronecker delta also ensures that the τ k+1/2 term is zero at the free surface for layer n.The source terms for structure cells are also illustrated in Figure 3. Effects relating to stresses as a result of volumetric deformation are not considered necessary to include due to their minor influence [15].For simplicity, wind friction effects on the free surface are also ignored however, wind friction effects can be easily added should the required wind data be available and the effects deemed necessary to include.
The domain is divided into cells (V i ) i∈Z and the discretised first order finite volume scheme is given by: Where the subscript i represents the ith cell, the superscript n represents the nth time level and ∆x and ∆t represent the cell size and time step respectively.F i−1/2 and F i+1/2 represent the numerical fluxes at Figure 3: Annotation of the source terms for example structure cells and their component sub-cells on uneven bed topography.R represents a reaction force induced as a result of the uneven bed topography, τ represents a friction force acting at a layer interface, z denotes the elevation above the bed and h denotes the water depth in the sub-cell.U i is the vector of conserved variables for the ith whole cell, which is equal to the sum of the conserved variables for the component sub cells U i,k .
the i ± 1/2 interfaces respectively.For the structure cells, it is the constituent sub-cells which are updated using the following modification of ( 12): Where U n i,k represents the conserved variables for the kth sub-cell in the ith structure cell at time level n.F i−1/2,k and F i+1/2,k represent the numerical fluxes at the kth layer of the i ± 1/2 interfaces respectively.Although a 1D scheme is implemented in this case, implementation as a 2D scheme requires no fundamental changes to the method.

Numerical Flux Computation
The process for resolving fluxes is dependent on the type of interface (NI, II or SI).For structure and intermediate interfaces Harten-Lax-van Leer (HLL) approximate Riemann solvers [16] are used to resolve the intercell numerical fluxes.For normal interfaces, other suitable approximate Riemann solvers may be used, however, HLL approximate Riemann solvers are recommended for consistency.

Normal Interfaces
A robust algorithm presented by Glenis et al. [17] is used to calculate wave speeds for the Riemann problem, which is outlined in Algorithm 1.Following calculation of the wavespeeds, a standard HLL approximate Riemann solver ( 14) is used to determine numerical fluxes across the normal interface.
N I As discussed prior, other suitable approximate Riemann solvers may also be used however, use of a HLL solver is recommended for consistency.

Structure Interfaces
At a structure interface the layers of flow and can be divided into open and closed as shown in Figure 5.
Open layers are considered as having a transmissive boundary at the structure interface, with the portion Open Open of the structure interface shared by the adjacent sub-cells having no influence on the exchange of conserved variables.Closed layers are considered as having a reflective boundary at the structure interface due to the presence of the structure.For each open layer, a single Riemann problem must be constructed and solved whereas, at each closed layer two Riemann problems must be constructed and solved, as shown in Figure 6.Solution of two Riemann problems for a closed layer is necessary to implement the reflective boundary condition at the structure interface, which reflects the flow in both the left and right sub-cells.This process is based on the assumption that the vertical velocity of the flow is negligible, which is a fundamental assumption for the derivation of the shallow water equations, and therefore the direction of the flow can be considered to be primarily parallel to the bed.
Algorithm 1: Calculation of wavespeeds [17].An initial approximation (h 0 ) of the depth in the star region (h * ) using a two-rarefaction approximate state Riemann solver is used to determine whether a two-rarefaction or two-shock approximation is optimal.For a multi-layer system, the wave celerity is defined as c i,k = g(h i,k + h i,k (+) ), where c i,k is the celerity for cell i layer k, h i,k is the thickness of cell i, layer k and h i,k (+) is the depth of water in cell i above layer k.
Use two-rarefaction approximate state Riemann solver Use two-shock approximate state Riemann solver x t 0  The numerical flux for each layer is determined by applying (4) to each layer, where the numerical flux for a layer is given as: Which can then be used to determine the flux at the interface using a standard HLL approximate Riemann solver (14).The method for determining the fluxes at a structure interface is summarised in Algorithm 2.

Intermediate Interfaces
In order to resolve fluxes with the adjacent sub-cells it is necessary to temporarily define layer properties for the intermediate cell as shown in Figure 7.The properties for the temporary layers in the intermediate interfaces are defined by assuming that the velocity in each layer is equal to the average velocity of the whole intermediate cell and that the depth in each layer is limited to the maximum depth capacity of the adjacent sub-cell.The fluxes for each layer can then be found using the process outlined for the open layers in Algorithm 2.
Algorithm 2: Calculation of fluxes for an example structure interface as shown in Figure 6.k is the index of the layer under consideration, n is the total number of layers at the structure interface.
Calculate wavespeeds Calculate layer flux For the closed layer h i+1,ghost ← h i,k Right ghost cell water depth u i+1,ghost ← −u i,k Right ghost cell water velocity calculate S − k L , S + k L using Algorithm (1) Calculate wavespeeds calculate F i+ 1  2 ,k L using (6) Flux for the left side of the structure Flux for the right side of the structure where u i−1 represents the average velocity for the whole intermediate cell.

Conservative Updating of Conserved Variables
Once numerical fluxes have been resolved across all interfaces within the computational domain, the final procedure for each timestep is to update the conserved variables contained within each cell and sub-cell.Normal cells are updated using equation (12), which is standard for a one-dimensional Godunov type scheme.

Normal Cells
For cases involving variable bed topography, a well-balanced treatment of the topographic source terms can be achieved via the hydrostatic reconstruction method [18] or via upwinding of the source terms [19].Suitable explicit or implicit treatment of the remaining source terms are both viable depending on the desired stability and admissible constraint of the stable timestep.For strong stability and the flexibility of advancing the solution at the timestep for the advection problem, the splitting method proposed by Liang and Marche [20] is recommended: The following simple limiter is also recommended to ensure stability in regions where the water depth approaches zero: Structure Cell Intermediate Cell Normal Cell The same procedure for updating a normal cell is applied to an intermediate cell however, due to the fact that fluxes at a intermediate interface are calculated on a sub-cell basis (Figure 9), they must first be summated.For this case illustrated in Figure 9 this is equal to:

Structure Cells
Structure Cells Intermediate Cell Since structure cells are divided into sub-cells, it is necessary to update each individual sub-cell using the respective left and right fluxes as per: Where U n i,k represents the vector of conserved variables for the kth sub-cell contained within the ith cell at time level n.F i−1/2,k and F i+1/2,k represent the left and right fluxes for the kth layer of the ith cell.As for the normal cells, a well-balanced treatment of the topographic source terms may be achieved via the hydrostatic reconstruction method or via upwinding of the source terms.The remaining source terms may be treated using suitable explicit or implicit methods depending on the desired stability and constraint of the timestep.For strong stability and the convenience of advancing the solution at the timestep for the advection problem, a point implicit scheme is recommended for the friction source terms: At the sub-cell interfaces containing structures there are two numerical fluxes as illustrated in Figure 10, as a consequence of the two reflective boundaries implemented at each side of the structure.Since not all of the external forces are accounted for, these fluxes may be unequal, with the difference in the sum of the fluxes at the left face of the structure interface (F (−) ) and the right face of the structure interface (F (−) ) equal to the resultant hydrostatic pressure force exerted on the structure multiplied by the ratio of the timestep to the cell width (∆t\∆x(F (+) − F (−) )).
Once the sub cells have been updated, their updated depth may exceed the maximum depth capacity for the layer and it is therefore necessary to re-define the layer properties of the structure cells in order to maintain alignment of the layers with the obstructions modelled at the interface.The process for redefining the layer properties is outlined in Algorithm 3, for which an illustrative example is also provided via Figure 11.Algorithm 3: Redefinition of the sub-cell properties based on the maximum depth capacity of the layers defined at a structure interface, post updating of the conserved variables.h and q represent the redefined depth and momentum.j refers to the index of the redefined layers and k refers to the index of the updated layer properties pre-redefinition.n is the maximum number of layers defined at a structure interface.
for each structure cell do

Model Validation
Previously published validation data [12], collected from experiments conducted in Newcastle University's Armfield S100 Research Flume, is used to validate the accuracy of the proposed Riemann solver.The S100 Research Flume is a 12.5m long, 1m wide, 0.8m deep flume capable of producing flow rates up to 400ls −1 .Using the control panel, shown in Figure 12, the user can select a desired flow rate which is then produced by the two pumps which draw water from the sump.The flow rate is maintained and corrected via a proportional-integral-derivative control loop, which uses a electromagnetic flow meter (Euromag Model MUT2200EL) to ensure that flow rate within the inflow pipe matches the desired flow rate.According to Euromag technical sheet [21], each sensor is calibrated on a hydraulic test rig equipped with an ISO17025 traceable weighing system, which ensures that the accuracy of the sensor is equal to 0.2% ± 2mms −1 with a repeatability of approximately 0.1%.A summary of the maximum permissible error limits for the instrument, provided by the manufacturer, is presented in Table 1.

Maximum Permissible Error limits for Euromag Model MUT2200EL DN 350 PN 10 EN 1092-1
Flow Rate Table 1: Maximum permissible error limits for the electromagnetic flow meter for a range of flow rates within the inflow pipe (adapted from [21] p.4).
The validation experiments consisted of running the flume at a range of flow rates, with a range of different barrier geometries placed within the flume cross-section, at a distance of 5m downstream.The flume tilt was set to 0% for all validation experiments in order to eliminate any potential numerical errors introduced as a result of topographic source terms.Once steady state flow conditions were achieved for each experiment, depth measurements were obtained using vernier point gauges.The full validation dataset is available as supplementary material from the referenced publication.

Numerical Setup
All numerical simulations were conducted on a 12.5m 1D spatial domain, discretised into a structured grid comprised of 0.1m cells (∆x = 0.1m).In order to ensure satisfaction of the Courant-Friedrichs-Lewy condition, a Courant number of C = (0.95∆x)/(S n max ) was used to determine a stable timestep, where S n max is the maximum absolute wave speed at time level n.Since the bed slope is set to 0% this has the intended effect of simplifying the source terms, only requiring the friction source term to be resolved, facilitating clearer analysis of the accuracy of the Riemann solver.The friction source terms for normal and intermediate cells are resolved using (16).The friction source terms for the structure cells are resolved using (20).A Manning's n equal to 0.012 and a kinematic viscosity of 1.0034 × 10 −6 m 2 s −1 is assumed for all numerical simulations.
The upstream and downstream boundary conditions are both implemented using exterior ghost cells.In order to replicate the constant inflow produced by the S100 flume, an inflow boundary condition is defined at the upstream end utilising relationships derived from the Riemann invariants across a rarefaction wave.At the downstream boundary a critical depth boundary condition is imposed.Full details for the implementation of the boundary conditions are presented in [12].

Results
The following validation test cases can be categorised into three primary flow configurations: • Flow under a barrier.
• Flow under a barrier, producing a downstream stationary hydraulic jump.
• Flow over and under a barrier.
Through comparisons between the experimental and numerical data for the six presented validation test cases, the suitability and accuracy of the proposed solver is demonstrated.

Flow Under a Barrier
For test case one and test case two, the solver produced accurate predictions for the upstream and downstream depth, capturing the interaction of the flow with the obstruction.In both test cases there is a slight overestimation of the upstream depth which equated to an error of 0.7−8.3% for test case one and 0.1−12.6%for test case two.In contrast, the downstream depth was slightly overestimated in both test cases, with a greater error for test case two due to the numerical prediction of a hydraulic jump at approximately x = 10m downstream.This is potentially a consequence of greater uncertainty in the measurement of the downstream depth, due to the presence of turbulent and unsteady flow at the outfall, which contributed to difficulty implementing the correct downstream boundary condition.Moreover, the location of a stationary hydraulic jump was determined to be extremely sensitive to small deviations in the flow during the execution of the lab experiments.
The velocity upstream of the barrier is predicted accurately for both test cases with errors in the region of 7 − 11%.The numerical estimation of the velocity at the upstream face of the barrier has a larger error however, this is a localised error, constrained only to the structure cell immediately upstream of the barrier.Since the discharge predictions are accurate otherwise, the overestimation of the downstream depth corresponds to an underestimation of the downstream velocity equating to an error of 2.4 − 18.8% for test case one and 0.8 − 20.6% for test case two.With the larger errors for test case two arising as a result of the incorrect prediction of the hydraulic jump.

Stationary Hydraulic Jump
Test case five and test case six showcase the capacity of the solver to accurately resolve stationary hydraulic jumps.The two presented test cases use the same barrier configurations with different flow rates, resulting in the formation of different stationary hydraulic jumps for each scenario.In both cases, the numerical results correctly predicted the formation of a stationary hydraulic jump downstream of the barrier.For test case six, the position and height of the jump was accurately captured.For test case five, the height of the jump was accurately captured, however, the formation of the jump was premature occurring at approximately x = 0.5m upstream of the actual location.The robust wave estimation algorithm (Algorithm 1) was determined to be crucial for accurately capturing and maintaining the stationary hydraulic jumps for the relevant numerical simulations.
In both cases the numerical results predict jumps with a zero length roller, characterised by a sharp discontinuity in the depth of flow at the toe of the jump, which is a feature of the classical shallow water equations; since there is no internal energy transcribed within the classical shallow water energy loss through a shock discontinuity is instead captured via Rankine-Hugoniot relations arising from the conservation of mass and momentum [22].This is insufficient to capture the complex behaviour which occurs within the transition region of turbulent hydraulic jumps with a Froude number of greater than 1.5.Methods to overcome the shortcomings of the classical shallow water equations, such as the work of Richard and Gavrilyuk [23], are not appropriate nor necessary for the desired application of flood risk modelling.
More generally, the predictions of the upstream and downstream depth and velocity proved to be accurate for both test cases, outside of the early prediction of the hydraulic jump for test case five.For test case five, there was a slight over estimation of the upstream depth corresponding to an error in the region of 1.4 − 12.3%.Ignoring the region of the domain occupied by the hydraulic jump (6 − 7m), downstream depth predictions were also found to be accurate with errors in the region of 0.7 13.2%.For test case six, the accuracy of the predictions starts to degrade towards the downstream boundary suggesting that the boundary condition may not be optimal.However, despite the increasing errors towards the downstream boundary, the solver still contributed to accurate results overall with depth errors from 0.2 − 19.2% and velocity errors from 0.2 − 23.6% for the data points between x = 0 − 8m, with errors increasing to 27.3% and 36.7% respectively at the boundary.

Flow Over and Under a Barrier
For test case eight and test case nine, depth predictions proved to be accurate with errors increasing towards the downstream boundary in both cases.For test case eight, upstream depth predictions were extremely accurate (0.1 − 3.3%).The upstream depth was overestimated for test case nine but remained accurate with errors in the region of 5.9 − 10.5%.For both test cases, there was a slight overestimation of the downstream depth with errors in the range of 6.5 − 24.7%.Figure 17 demonstrates that the water was observed as vertically flowing over the barrier for test case eight which cannot be captured by the numerical model, due to the nature of the fundamental equations and the structure of the finite volume scheme.Although this behaviour isn't captured by the model, the overall results remain accurate and the general behaviour is well captured.Certainly, for applications concerning flood risk modelling, the key quantities are the upstream and downstream depths which are observed to be consistent with the validation data.
The velocity predictions are similarly accurate with errors in the region of 1.2 − 19.8% for both of the presented test cases.As for the previous test cases, there is also a local error in the prediction of the discharge in the cells proceeding the barrier with discharge predictions otherwise proving accurate.Table 2: Comparison between the experimental data and simulated data for Test Case 6, for a range of mesh resolutions.
The data in Table 2 demonstrates negligible differences in the results for the tested mesh resolutions which range between 1 − 20cm (∆x = 1 − 20cm).The relevant plots illustrating the results can be found in the Appendices.The primary difference between the meshes is in the sharpness of the depth discontinuity at the toe of the stationary hydraulic jump, which becomes steeper as the mesh is refined.

Comparison
In order to demonstrate the comparative value of the solver presented within this paper, designated as solver 2, a comparison is presented with the solver presented in [12], designated as solver 1.A comparison of the solvers for test case one, shown in Figure 19, demonstrates accurate results for both solvers.Solver 2 has a marked increase in accuracy for the depth upstream of the barrier, whereas, the depth and velocity downstream of the barrier is captured slightly more accurately by Solver 1. Similarly for test case 8, shown in Figure 22, there is a improvement in the prediction of the upstream depth for Solver 2, with comparatively accurate results for both solvers downstream of the barrier.The benefits of Solver 2 are however, showcased best via Figure 20 and Figure 21.Whilst Solver 1 is able to broadly capture the upstream and downstream depths, which is of primary concern for flood risk management applications, Solver 2 captures the upstream and downstream depths more accurately including the formation of a stationary hydraulic jump.The difference in the accuracy of the velocity predictions downstream of the barrier is stark and demonstrates that the superiority of the solution procedure utilised by Solver 2 with regards to capturing the horizontal velocity profile in the vertical plane at the structure interface.
Although it is clear that Solver 2 is capable of producing superior results, there is a clear increase in complexity and computational burden in comparison with Solver 1.However, since the implementation is local, on a sufficiently large domain the difference in computational efficiency of the two solvers is unlikely to be a limiting factor since structure cells are likely to comprise a very small percentage of all cells within the computational domain.As such, the primary grounds for the use of Solver 1 over Solver 2 should be limited to scenarios in which the simplicity of implementation is important and for use cases which predominantly involve supercritical downstream flow regimes.Otherwise, Solver 2 proves to be the optimal solution.Furthermore, the capacity for Solver 2 to accurately resolve stationary hydraulic jumps and more accurately capture the velocity at the barrier presents further opportunities such as the modelling of the transport of water soluble contaminants.Passive scalars, such as water soluble contaminants, are assumed to be passively advected with the fluid and via the reintroduction of the contact discontinuity wave, via switching from a HLL to a HLLC (Harten-Lax-van Leer contact) approximate Riemann solver [24], their transport can be modelled.Since this process is highly dependent on the accurate determination of the velocity, this is only possible for Solver 2. This has important applications in terms of modelling water quality, especially since the combination of flows around obstacles and species equations is seldom explored and is therefore to be the subject of further work.Figure 20: Comparison between solver 1 and solver 2 with respect to the experimental results for test case 5. Details for solver 1 can be found in [12].Details for the numerical setup can be found in Section 3.1.

Conclusion
A new Riemann solver, capable of resolving numerical fluxes across a partially obstructed interface, has been presented.Via the validation process, it has been demonstrated that the solver is able to adequately capture fluid-structure interactions for a range of barrier configurations and flow rates.Furthermore, via the comparison process, it has been demonstrated that the solver represents a significant improvement on the previously published solver [12].It is clear that the new solution procedure addresses the identified weakness of the previous solver by more accurately capturing the vertical variation in the horizontal velocity profile at a structure interface.This results in the more accurate determination of the flow characteristics including the ability to resolve the location and jump height of stationary hydraulic jumps.However, this does come at the cost of increased computational demands and implementation complexity although, due to the local nature of the solution procedure and the proportionally small number of structure cells within a computational domain, the increase in computational expense in unlikely to be prohibitive.As for the previous solver, the biggest barrier to implementation is the scarcity of the required data for structures and the availability of suitable meshing algorithms, which remains the subject of further work.
The capability of the solver to resolve numerical fluxes across a partially obstructed interface has significant implications for modelling a variety of structures within two-dimensional hydrodynamic models.This has applications in terms of improving flood inundation modelling capabilities as well as enabling the modelling of infrastructure resilience modelling and the structural health monitoring of hydraulic structures.Moreover, the more accurate determination of the fluid velocity in comparison with the previously presented solver presents new opportunities such as the capacity to model hydraulic jumps and the ability to integrate species equations enabling the modelling of water soluble contaminants in conjunction with flows around obstacles.
As for all models, the underlying assumptions must be considered in order to ascertain the limitations of the model and as such the solver should be considered appropriate for modelling structures at a spatial scale whereby approximating the structure as a partial obstruction existing at a cell interface is appropriate.Although, the proposed model does not capture all of the energy losses which occur as a result of the fluidstructure interaction, such effects are insignificant at this spatial scale in comparison with the effect induced by the blockage of the flow by the structure, which is well captured as shown by the validation results.For detailed analyses of individual structures 3D CFD analyses are recommended.
Avenues for further development of the solver are limited without compromising on the compatibility of the solver with standard numerical schemes utilised in contemporary hydrodynamic models.The layer redefinition process, particularly where layers are shifted downwards a significant distance, presents the greatest weakness of the method.However, such cases involve flow which is inherently vertical in nature and it is difficult to consolidate this with the fundamental nature of the shallow water equations which underpin modern hydrodynamic modelling.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Figure 1 :
Figure 1: A simple computational domain [a, b] illustrating the designation of structure, intermediate and normal cells with their corresponding interfaces.z 1/2 and z 3/2 represent the height above the bed of the base and cover of the idealised structure represented at the structure interface.

Figure 4 : 5 2
Figure 4: (a) Example normal interface with adjacent normal cells and (b) the general structure of the general solution of the Riemann problem for a normal interface.S − is the left wave speed and S + is the right wave speed, as defined in Algorithm 1. h * and u * denote the conserved variables in the star region.F i− 5 2

Figure 5 :
Figure 5: Designation of open and closed layers at a structure interface.
Division of the structure cells into sub-cells and their respective properties.The subscripts L and R are used to differentiate between the left and right face of the structure interface.
b) The general structure of the general solution of the Riemann problems for an example structure interface shown in (a).The introduction of fictitious ghost cells for the purpose of implementing reflective boundary conditions are denoted by grey shading.

Figure 6 :
Figure 6: Method for resolving fluxes for the sub-cells adjacent to a structure interface.

Figure 7 :
Figure 7: Temporary division of an intermediate cell into layers in order to resolve fluxes at a intermediate interface.ui−1,1 = u i−1,2 = u i−1,3 = u i−1where u i−1 represents the average velocity for the whole intermediate cell.

Figure 8 :
Figure 8: Illustration of the numerical fluxes at the normal interfaces bordering a normal cell.

Figure 9 :
Figure 9: Illustration of the numerical fluxes used to update a intermediate cell.

Figure 10 :
Figure 10: Illustration of the numerical fluxes used for updating the sub-cells of which a structure cells is comprised.

Figure 11 :
Figure 11: Illustration of the layer redefinition process post updating of the conserved variables.The redefinition process is required to re-align the updated properties of the sub-cells with the respective boundary conditions implemented at the structure interface.

Figure 12 :
Figure 12: Integrated control panel for the S100 Research Flume, including a schematic of the flume.Two pumps, which draw water from a recirculating sump, supply water to the flume via a pipe connected to the upstream (right) end.At the left end of the flume, the water exits the flume via a sloped free outfall into the the sump.

Figure 14 :
Figure 14: Comparison between numerical and experimental results for test case 2. Details of the numerical setup can be found in Section 3.1.

Figure 15 :
Figure 15: Comparison between numerical and experimental results for test case 5. Details of the numerical setup can be found in Section 3.1.

Figure 16 :
Figure 16: Comparison between numerical and experimental results for test case 6.Details of the numerical setup can be found in Section 3.1.
Figure 17: Comparison between numerical and experimental results for test case 8. Details of the numerical setup can be found in Section 3.1.
Figure 19: Comparison between solver 1 and solver 2 with respect to the experimental results for test case 1.Details for solver 1 can be found in Details for the numerical setup can be found in Section 3.1.