Effects of reconstruction of variables on the accuracy and computational performance of upscaling solutions of the shallow water equations

This paper presents a new sub-grid flood inundation model obtained by upscaling the shallow water equations (SWE) to enhance the model efficiency in large-scale problems. The model discretizes study domains using two nested meshes. The equations are solved at the coarse mesh by a second-order accurate in space (i.e. piecewise linear reconstruction of variables) Godunov-type finite volume (FV) method, while the fine mesh is used to incorporate high-resolution topography and roughness into the solution. The accuracy and performance of the model were compared against a first-order version of the model recently proposed by the authors and a second-order conventional FV model using artificial and real-world test problems. Results showed that improved accuracy is delivered by the proposed model, and that at low-resolution meshes, the spatial reconstruction of variables of the numerical scheme plays a major role in the solution's accuracy.


Introduction
Computational simulation of flood inundation is now a central component of a wide range of critical applications, such as flood risk assessment and mitigation, emergency response, and engineering design. The vast majority of codes developed by industry and research institutions to perform these simulations are based on the numerical solution of the system of shallow water equations (SWE) (e.g. Alcrudo & García-Navarro, 1993;Cea et al., 2010;Deltares, 2019;DHI, 2017;Guinot et al., 2017;Kahl et al., 2022;Kesserwani & Liang, 2015;Liang, 2010;Mignot et al., 2006;Sanders & Schubert, 2019, to cite only a few). While existing SWE models (in particular those implemented in two horizontal dimensions, or simply 2D) are able to accurately and robustly capture the most relevant characteristics of flood inundation, their high computational cost still represents a barrier to a large number of important applications. This computational constraint has restricted the application of these models to even relatively small domains and short duration flood events (Liang & Smith, 2015;Ming et al., 2020), but is even more problematic for large-scale inundation simulations, which are needed to inform important decisions. This is despite the recent progress in the development of high-performance computational modelling using computational parallelization techniques (e.g. Morales-Hernández et al., 2021;Sanders & Schubert, 2019;Vacondio et al., 2014;Xia et al., 2019). The numerical solution procedure has a significant impact on model efficiency (e.g. Liang & Smith, 2015;Yu et al., 2015). Typically, computational methods for the solution of flood inundation offer a compromise between efficiency, accuracy and complexity of a model.
It is widely recognized that high-resolution topography data can substantially improve the accuracy of inundation models. On the other hand, the computational time (C) of a simulation substantially increases as the model resolution is refined (e.g. typically, for an explicit scheme, C ∼ x −3 , where x is the size of a square grid cell). To overcome the challenges associated with the computational cost of high-resolution simulations, a new class of dual-grid models has been developed for flood inundation modelling (e.g. Sanders & Schubert, 2019;Shamkhalchian & de Almeida, 2021;Stelling, 2012;Volp et al., 2013), which has been observed to provide efficient solutions for large-scale simulations (e.g. Sanders et al., 2022;Schubert et al., 2022). This approach may take advantage of Godunov-type finite-volume (FV) methods, which have been the object of extensive research over the last decades. Godunovtype models require assumptions about the spatial structure of the free surface elevation, depth, discharge per unit width, velocity, and topographic heights, which may have an important impact on the accuracy of solutions. For example, Begnudelli et al. (2008) performed an interesting comparison of two Godunov-type (FV) models, namely (i) a first-order accurate (i.e. piecewise constant spatial reconstruction of variables) FV scheme for the solution of the homogeneous SWE combined with a second-order model of topography (i.e. piecewise linear change in bed elevation) and (ii) a second-order accurate scheme for the homogeneous equations (i.e. piecewise linear reconstruction) combined with a first-order representation of topography (i.e. piecewise constant). The results of this comparison indicated that in selected practical cases, the former can be more accurate, efficient and robust than the latter.
Compared to single-grid models, dual-grid models may result in several different types of spatial structure in the solution. Figure 1 compares a FV single-grid against a FV dual-grid. In single-grid models, topography is typically represented as piecewise constant across a domain of study, although a piecewise linear approximation has also been used (e.g. Begnudelli & Sanders, 2006;de Almeida et al., 2018). On the other hand, in FV dual-grid models, topography within the computational cells is not prescribed by mathematical expressions, but instead is represented by the actual data defined at the resolution of the fine grid. Similarly, in FV single-grid models, the conservative variables (such as the water depth and unit width discharge) Figure 1 Bed and water volume representation in dual-grid and single-grid models. In single-grid and dual-grid models, water volume representation is typically based on water depth and water surface level, respectively. Dual-grid models adopt high-resolution bed level, while a single-grid model makes use of the average bed level for the solution of the homogeneous SWE. are typically reconstructed as piecewise linear or piecewise constant. In FV dual-grid models, the solution is typically provided for the non-conservative variable water surface elevation, which can be reconstructed as piecewise constant or piecewise linear. In both cases, the water depth varies within a computational cell. This implies a non-uniform distribution of velocity inside a computational cell even if the unit width discharge is assumed piecewise constant. These differences between FV single-and dual-grid, which are related to spatial reconstruction of topography and flow variables, may lead to very different solutions. The aim of this paper is to assess how accuracy and computational cost vary across different assumptions of spatial structure using a range of grid resolutions. This assessment will be conducted by implementing these assumptions in the dual-grid model by Shamkhalchian and de Almeida (2021). Shamkhalchian and de Almeida (2021) recently proposed and tested a new sub-grid model (SG) that makes use of two nested meshes, an approach that has also been previously used by Volp et al. (2013), Hénonin et al. (2015), and Sanders and Schubert (2019). The main differences among these models relate to the methods and the assumptions adopted to upscale the water depth (or water surface elevation) and velocity (or discharge) from high-to low-resolution, and the numerical techniques used, in addition to other, more specific approaches that are not described here. For example, Volp et al. (2013) used the method of Finite Volume to solve the SWE, while the upscaling of water depths and velocities was performed by spatially averaging the variables over the wet part of coarse cells. Friction was upscaled under the assumptions of a constant friction slope and uniform flow direction. In the model presented by Hénonin et al. (2015), water depths were also averaged over the wet area of the computational cell under the assumption of a piecewise constant free surface elevation, and the SWE were solved using the Alternating Direction Implicit Finite Difference scheme. Mass fluxes and friction were estimated at coarse resolution based on the averaged quantities. Sanders and Schubert (2019) solved the SWE using the Finite Volume method and also defined variables (water depth and unit width discharge) that were averaged over the wet part of computational cells. Fluxes along the edges of the coarse cell were computed at fine resolution by downscaling depths using a piecewise linear reconstruction of the free surface combined with local (i.e. fine resolution) topography, while the values of the velocity on the left and right of the edge are defined by the average over the corresponding coarse grid grid (B. Sanders, private communication). Friction was estimated based on the assumptions of constant friction slope and water surface elevation over the coarse grid. Shamkhalchian and de Almeida (2021) used the FV approach with fluxes along the edges of the coarse cells computed at the fine resolution grid [as in Sanders & Schubert, 2019], but in this case the Riemann problem was solved with both depths and velocities that were downscaled from coarse to fine resolution. Also, friction inside large cells was defined through the integration of the flow resistance equation over the coarse cell under the assumption of a uniform flow direction only. A brief outline of this model is presented in the following section, while further details can be found in Shamkhalchian and de Almeida (2021). The models described above were designed to reduce errors that typically arise with grid coarsening by incorporating topographical information at finer resolution. However, they cannot capture the effects of a complete blockage of the flow by, for example, unsubmerged topography bisecting a coarse cell (e.g. slender elements such as a dyke, or flood defence). New methods have been recently proposed to address this specific condition, (e.g. Ferrari & Viero, 2020;Ferrari et al., 2019;Viero, 2019).
This paper is structured as follows. Section 2 outlines the methodology adopted in the SG model to solve the governing equations. The first-order SG model (piecewise constant reconstruction) was introduced in Shamkhalchian and de Almeida (2021), and is extended to second-order spatial accuracy (piecewise linear reconstruction) in this paper. In Section 3, the artificial and real-world test cases utilized in this paper are described and the results of numerical simulations are presented and analysed. Finally, the main new insights and conclusions of the research are synthesized in Section 4.

The SG model methodology
This section provides a concise description of the SG model, including the extension of the Shamkhalchian and de Almeida (2021) approach to second-order spatial reconstruction of variables. A more detailed description is provided in Shamkhalchian and de Almeida (2021).

Governing equations and nested meshes used in the solution
The SG model solves the following integral form of the 2D shallow water equations to predict the water surface elevation and depth-averaged unit width discharge in time and space.
t is time, η denotes water surface elevation, h represents water depth, denotes the polygonal area (e.g. the quadrilateral area of a fully submerged computational cell) of the domain over which the equations are solved, the boundary of which is represented by , e = [e x , e y ] is the outward unit vector normal to , g is the gravitational acceleration, S 0x , S 0y , S f x , S f y , q x and q y are the x and y components of bed and frictional slopes, and depth averaged unit width discharges, respectively. The frictional slopes are defined based on Manning's expression as where u and v are the x and y components of the velocity vector, the magnitude of which is V = √ u 2 + v 2 . The governing equations are solved using a structured (in the current study, square and rectangular meshes, although the methodology can be easily extended to other mesh types) coarse computational mesh within which another structured, high-resolution mesh is nested (Fig. 2). The coarse and fine meshes form cells that are hereafter also referred to as large cells (computational cells) and small cells, respectively. Topography (i.e. bed elevation) and roughness (i.e. Manning's n M ) are constant within each small cell, and the arrangement of small cells inside a large cell determines a non-uniform distribution of these properties ( Fig. 2 shows an arbitrary rectangular computational domain). The geometrical and other relevant variables are defined as follows. The size of computational domain is The dimensions of the large and small cells along the x and y directions are denoted by x, y, δx, δy, respectively. The edges of the computational cells are denoted by 1 to 4 anticlockwise. The small cells are arranged in J columns and K rows within a large cell. z represents the bed level, which is defined at fine resolution. Large cells are specified by the subscript i, 1 ≤ i ≤ N , where N is the number of large cells in the computational domain. The location of a small cell inside a large cell is defined by the indices j and k (i.e. column and row, respectively). As an example, z i | j ,k specifies the bed elevation at the j th column and kth row of small cells within the ith large cell.

Solution updating procedure
The governing equations are solved using the FV method. In order to improve the accuracy of the solution performed at a coarse resolution (i.e. large cells), each term in the SWE equations is upscaled, allowing for finely resolved information to be included in the solution. The methods for upscaling the equations are described in detail in Shamkhalchian and de 412 A. Shamkhalchian and G. A. M. de Almeida Journal of Hydraulic Research Vol. 61, No. 3 (2023) Figure 2 Nested meshes adopted in the SG model. Almeida (2021) and are only briefly outlined in the following sections. Upscaled properties (i.e. defined at coarse resolution through the integration over a large cell of quantities distributed at fine resolution) are hereafter denoted by the overbar symbol, e.g.
The solution is indirectly updated by one time step t from time level n to n + 1 via an intermediate state or fractional step (LeVeque, 2002), which is denoted by the superscript *. In this intermediate state, the solution is advanced for all the terms of the governing equations except for S f ; the solution at * is then used to integrate the friction term using a semi-implicit method (de Almeida et al., 2018;Liang & Marche, 2009;Sanders & Schubert, 2019) detailed in the next sections. The fractional time-step integration has been widely adopted to prevent numerical instabilities (e.g. Liang & Marche, 2009), and is defined by: where Eq. (5) is obtained from the first-order Euler method (therefore, the methods described in this paper are only first order in time). The time step t is determined dynamically by the Courant-Friedrichs-Lewy (CFL) condition (e.g. LeVeque, 2002;Shamkhalchian & de Almeida, 2021;Toro, 2001) at the level of coarse cells, which are the computational cells.
Fluxes in the first integral of Eq. (5) are computed along each side of large cells by solving the Riemann problem at the scale of small cells, via the efficient Harten, Lax and van Leer Contact (HLLC) Riemann solver (Toro et al., 1994). These fluxes are then integrated along each edge of a large cell to approximate the total flux for the edge. This local Riemann problem requires the reconstruction of the variables in the vector U n i at fine resolution, i.e. variables must be downscaled from the solution obtained at coarse resolution, as described in the next sections. The treatment of Eq. (6) and the second integral in Eq. (5) are detailed in the Section 2.4.

Reconstruction of variables and flux computation
One key difference between the SG model and conventional (i.e. single-mesh) finite volume models based on the SWE is that the SG model requires two reconstructions of variables, while a conventional model only needs one. The vector U n i is first reconstructed at each computational edge (i.e. U n i | m , where m = 1, 2, 3, 4) using conventional methods (e.g. piecewise constant or piecewise linear reconstruction of variables over a large cell). In the second step of the reconstruction (only performed in the SG model), U n i | m is distributed (i.e. downscaled from coarse to fine resolution) over the large cell edges to define values of variables at small cells along the edges (e.g. U n i | j ,1 where j = 1, 2, . . . , J for 1 ). This distribution is implemented based on topography and roughness at the resolution of small cells.
The order of spatial accuracy of the model is hereafter defined in the conventional way, which depends on how the first reconstruction of variables is performed. A piecewise constant reconstruction within a large cell (U n i | m = U n i ) yields a firstorder accurate scheme. A second-order spatial accurate scheme can be obtained by a piecewise linear reconstruction. In this paper, the slope of the linear reconstruction is limited by the minmod method (LeVeque, 2002) to improve numerical stability. Also, for the same reason, the order of accuracy of the scheme is always reduced to first order around the wet-dry fronts. This means that in a wet large cell (fully or partially wet) next to a fully dry cell, the reconstruction is piecewise constant.
For the second step of the reconstruction (i.e. downscaling along the edges), the first component of U (i.e. η) is set constant across each edge. This results in a non-uniform distribution of h = max(η − z, 0) along the edges, which depends exclusively on the fine resolution topography. The second component of U n i (i.e. q n x i ) is distributed along the edge 4 (as shown in Fig. 3) as follows (the process for q n y i distribution is identical, except for the direction considered). The method is based on the assumption of constant friction slope, which is herein defined using Manning's relation (Burguete et al., 2007;Chow, 1959;Cunge et al., 1980;Viero & Valipour, 2017), and yields the following distribution: where q n x i | 1,k represents the unit width discharge (at time level n ) in a small cell neighbouring the edge 4 of a large cell, N n 4 is the number of wet small cells next to the edge 4 (which is a function of η n i ), ψ a is the conveyance of the edge 4 , which is approximated as (see Shamkhalchian & de Almeida, 2021 for wherez n i | 4 is the average of submerged small cells bed levels at the edge 4 at time level n. The parameters, T n 1 i | 4 , T n 2 i | 4 , and T n 3 i | 4 are defined as: Equation (7) has been written for a particular edge as an example; the formulation for the other edges is identical. The parameters T n 1 i to T n 3 i on each edge are dynamic and depend on N and therefore, on η n i . The estimation of these parameters could impose a high computational burden to the model as it requires summations at fine resolution, which would also need to be repeated at each time level. To improve the efficiency of the model, these parameters are computed for all possible values N and stored in sorted vectors at pre-processing. They are retrieved during the simulation as a function of the water surface level at very low computational cost. The reconstruction of the variables h, q x , q y , and z on both sides of a large cell edge at high-resolution defines the states of the Riemann problem at fine resolution, which is solved by the HLLC Riemann solver to approximate the mass and momentum fluxes. Integration of fluxes along each large cell's edge provides the total fluxes, which is substituted into the first integral of Eq. (5) to complete an important step to solve the homogeneous part of the 2D SWE.

Source terms
The bed slope source term, which appears in the second integral of Eq. (5) where K = y/δy, ź x | n k = (ź | 1,k −ź | J ,k ) n i andẑ| n k = 1 2 (ź | 1,k +ź | J ,k ) n i ,ź(z, η) = min(z, η). The estimation of the bed slope source term is based on the assumption that water surface level is constant within each large cell. This means that bed slope source term depends on the water surface level, the mesh resolution and topographic data next to the large cell edges. Topographic data that is not adjacent to the edges of the large cell does not play any role in the estimation of this term, leading to improved computational efficiency.
The friction source term is treated via a widely used implicit scheme, which updates the solution from U * i to U n+1 i via Eq. (6) (e.g. Cea & Bladé, 2015;de Almeida et al., 2018;Kesserwani & Liang, 2012;Liang & Marche, 2009;Sanders & Schubert, 2019). The method is adapted in the SG model in order to include fine resolution data into the solution, and thus to improve the accuracy of the model. To this end, the conservative variables are non-uniformly distributed (i.e. downscaled) over a large cell. Then, Eq. (6) is solved at the high-resolution mesh and finally, the high-resolution solutions are upscaled to the large cell (further details are available in Shamkhalchian & de Almeida, 2021), yielding: where q = [q x , q y ] is the unit width discharge vector, and: , p = 4, 5, 6 (12) T 4 , T 5 and T 6 are computed at pre-processing as a function of N w (i.e. the number of wet small cells within a large cell) and are then retrieved efficiently from the sorted vectors during the simulation. It should be noted that the friction updating only affects the unit width discharge components of U * i and has no effect on η. In the process of updating the solution for the friction source term, the water surface level is treated as piecewise constant over computational cells.

Test cases and numerical results
In the SG model described in this paper, topography is represented at a resolution that is typically much higher than the resolution at which the conservative variables vector is reconstructed. This has two main implications: (i) the model's estimation of terrain elevations is very accurate (and for the purpose of this analysis may be assumed nearly exact) and (ii) this estimation is independent of the resolution of the computational mesh (i.e. large cells). This enables the effect of the computational mesh resolution on the solution of the homogeneous SWE (obtained with first and second-order models) to be separated from the corresponding effect it would have on how bed elevations are represented in the model (e.g. by a piecewise constant or piecewise linear approximation).
In this study, all simulations performed using the SG model are compared to those obtained using the so-called traditional (T) model. The T model is a conventional (i.e. single-mesh) Godunov-type FV model, which also uses the HLLC Riemann solver to solve the 2D SWE over regular quadrilateral cells. It is a particular case of the SG model for the condition x = δx and y = δy (i.e. there is only one small cell within each large cell). In this section, artificial and real-world test cases are employed, which are simulated via four methods: (1) SG1 -the SG (i.e. dual-mesh) model, which implements a first-order spatial reconstruction of flow variables over a detailed terrain model introduced in Fig. 2; (2) SG2 -the second-order accurate SG model (i.e. piecewise linear spatial reconstruction) for the solution of the homogeneous SWE over the detailed terrain model; (3) T -the single-grid model implementing the second order reconstruction of flow variables (i.e. piecewise linear spatial reconstruction of the vector [η, hu, hv]), where the bed levels are assumed piecewise constant within each computational cell; and finally (4) Td -the single-grid numerical solution similar to the T model, but which reconstructs the vector [h, hu, hv] (i.e. the flow depth is reconstructed instead of water surface level) piecewise linearly within a large cell. All four models are only first-order accurate in time.
In the next sections, the following notation is adopted to label the solutions. The labels defined by SG1, SG2, T and Td are used to describe the methods defined in the previous paragraph, while the mesh-resolutions are shown in brackets. For example, SG1 (50/5) represents the solution obtained by the first-order SG1 model in a simulation using 50 m and 5 m coarse and fine resolution meshes, respectively. The second-order T model solution at the mesh resolution 300 m is denoted by T (300).

Test Case 1: 1D steady subcritical flow over bed with large elevation changes
The first test is for steady subcritical flow in a 1D rectangular and prismatic channel with varying bed elevation. This test case is aimed at assessing the performance of the models in problems involving relatively large variation of topography, such as channels with sections of adverse (negative) bed slope. The test case (in particular the longitudinal variation of bed elevations) and the corresponding analytical solution were defined by the inverse method proposed by MacDonald (1996); the water depth profile, flow discharge (Q), channel width, and roughness are assumed and then, using the 1D steady flow form of the SWE, the expression for the bed profile is derived. In this test, a rectangular 3000 m long, 10 m wide channel is used, where the Manning number n M = 0.05 s m −1/3 is constant and the flow discharge is Q = 16 m 3 s −1 . The water depth profile is defined as h(x) = 1.5 + 0.4 sin 5 ( π x 550 ). The water depth and Froude number (F) vary between 1.1 m h 1.9 m and 0.2 F 0.44, respectively. Simulations were performed for computational cells with resolutions of 5, 10, 20, 50, 100, 200 and 300 m. The fine mesh resolution for the SG model was 1 m for all the simulations. Figure 4a shows the water surface elevation predicted by SG1, SG2, T and Td models in Test Case 1. For clarity, only the solutions at low-resolution computational meshes are displayed, as the high-resolution results were all very close to the analytical solution (this is further shown and discussed next). These results show that at low resolution, SG2 yields the best accuracy followed by Td, T and then SG1.
Further information on the accuracy obtained with each method is provided by Fig. 4b, where the root mean square error (RMSE) of the solutions of η relative to the analytical solution is presented. According to Fig. 4b, at high-resolution meshes (e.g. x 20 m), the values of RMSE are similar for all the methods. This would imply that a first-order accurate scheme over detailed terrain model (i.e. SG1) delivers solutions with accuracy similar to those by a second-order accurate solution (SG2, Td and T), which is in agreement with the observation by Begnudelli et al. (2008). On the other hand, the influence of the spatial order of the numerical method used to solve the homogeneous SWE on the solution becomes distinguishable at coarse resolutions (i.e. x > 20 m for this case). This finding was also confirmed in many (but not all) other similar test cases, the results of which are not reported in this paper to avoid repetition. For example, at x = 300 m, the RMSE of SG1 is about 0.6 m, while for SG2, a much smaller value of ∼0.1 m was observed. This has important implications for large-scale flood inundation problems, where simulations at coarse resolution are adopted to reduce the computational cost. Thus, solutions such as SG2, which benefits from a second-order reconstruction of flow variables and the upscaling methods to that make use of high-resolution topography data, may lead to improved accuracy compared to the SG1 model, which uses a first-order reconstruction of flow variables and the same upscaling techniques. In this test, the single-mesh model using a second-order reconstruction of the water depth (Td) provided better accuracy when compared to the corresponding reconstruction of the water surface elevation (T), as was also previously observed by Buttinger-Kreuzhuber et al. (2019).
The above results indicate that at low resolution, the secondorder reconstruction of variables is as important as (or even more important than) the use of detailed topography [e.g. compare the methods SG1 and T (Td) in Fig. 4a]. The reason for this may be explained as follows. Assume for example, the common problem of flow down a positive slope. Figure 5 illustrates this scenario, as well as the corresponding first and second-order accurate reconstruction of variables (herein, the main focus is on water surface) in the computational cells of the SG (Fig. 5a) and traditional models (i.e. T and Td, Fig. 5b). For the T and Td models, the bed level is constant and defined as the averaged bed level within the large cell (Fig. 5b), while the SG model represents the detailed topography using the small cells (Fig.  5a). A first-order reconstruction of η combined with the detailed topography (blue line in Fig. 5a) would imply an increase in depth from left to right (in the direction of the flow) for the SG model's cell, while the depth would be modelled as constant within the cell in the T and Td models (i.e. Fig. 5b). On the other hand, the improved representation of the water surface obtained with a second order reconstruction combined with the detailed topography in the SG model yields a more realistic approximation of depths on the edges of the computational cell (see red line in Fig. 5a). It should be noted that the relative importance of the second-order reconstruction of the variables, compared to the detailed representation of topography, becomes relevant at low resolution, where, η (in Fig. 5) is typically large. These results highlight the need for coherent methods to represent the bed topography and the reconstruction of flow variables.

Test Case 2: real-world flood inundation problem
Test Case 2 is aimed at comparing the performance of the modelling methods described in this paper when simulating a real-world flood inundation problem. The test simulates a flood event that occurred in the River Tiber between 27 November and 1 December 2005 and lasted for 113 h. The River Tiber is an important river in Italy which flows between Apennine Mountains and Tyrrhenian Sea and has a total length of approximately 416 A. Shamkhalchian and G. A. M. de Almeida Journal of Hydraulic Research Vol. 61, No. 3 (2023)    The initial conditions of the model were defined through a steady flow simulation using the boundary conditions immediately before the event (e.g. Q = 374 m 3 s −1 ). The unsteady boundary conditions are presented in Shamkhalchian and de Almeida (2021). The values of Manning's coefficient are 0.035 m s −1/3 and 0.0446 m s −1/3 for the main channel and floodplain, respectively, as previously defined by Morales-Hernández et al. (2016). A 5 m × 5 m resolution digital elevation model (DEM) is used to produce the terrain models required for the T and proposed SG model simulations. The computational mesh resolutions are 20, 40, 50 and 100 m, while the fine mesh used for all the SG model simulations is the 5 m resolution DEM. The solution of the T model at the resolution 4 m is hereafter referred to as the benchmark solution. Water surface levels recorded (field data) during the flood event at two cross-sections C1 and C2 are available, which are shown in Fig. 6. Figure 7 shows the measured and predicted time series of the water surface elevation at cross-sections C1 and C2 during the 113 h flooding event. For the sake of clarity, the results for the mesh with 40 m resolution are not displayed in Fig. 7. Furthermore, results for T (100) are not shown because the model was unable to complete the simulation of this problem at this (and coarser) resolution. This is because at low resolution (e.g. 100 m), the piecewise constant terrain elevation (obtained by averaging) leads to cell bed levels that are higher than the water surface level that is set as the downstream boundary condition. Figure 7 shows that at low resolution, the set of runs obtained by the SG1 and SG2 models provides more accurate results than those from the single-mesh T model. The advantage of the sub-grid approach is particularly clear considering that T (50) benefits from a computational mesh resolution twice finer than that used in SG1 (100/5) and SG2 (100/5). For example, at t = 113 h, the water surface errors (relative to the benchmark solution) of SG1 (100/5), SG2 (100/5) and T (50) are 0.51 m, Figure 7 Test Case 2. Numerical solutions of water surface elevation at cross-sections: (a) C1, (b) C2 obtained by different models, along with the measured field data, and benchmark solution. .04 m at cross-section C2, respectively. The main difference between Test Cases 1 and 2 is that the domain (and therefore the computational cells) are fully wet in the first test, while in test case 2 it is partially wet, and the SG models benefit from the algorithm to account for partially wet cells. The possibility of simulating the effects of partially wet cells, along with the inclusion of high-resolution data at sub-grid scale in the process of solution, leads to the improved accuracy delivered by the SG models. Further analysis of the accuracy of the numerical solutions is presented in Fig. 8, where the water depth RMSE (relative to the benchmark solution) of the results are plotted as a function of the mesh resolution. The results in Fig. 8 are at t = 110 h, when large differences appear between the numerical simulations and the benchmark solution ( Fig. 8). The water depth is obtained by subtracting the 5 m resolution bed elevations from the corresponding water surface level obtained by the simulations. Figure 8 shows that the RMSE from the T simulations is much higher (especially at coarse resolution) compared to the two other sets of runs with the SG models. This shows the importance of high-resolution data in the solution process. At all mesh resolutions, the accuracy of results from the SG2 model is significantly better than those by the SG1 results at the same resolution, which, as previously indicated in test case 1, highlights the importance of higher order numerical schemes at coarser resolution SG simulations. Based on Figs 7 and 8, the Td model provided slightly more accurate results relative to the T model in this test. Table 1 presents the computation time of the models at different resolutions. The simulation cost decreases faster (with increasing cell size) for the traditional models (i.e. T and Td) than for the SG models (i.e. SG1 and SG2). This is because the SG models need to compute the mass and the momentum fluxes across the interface of computational cells by solving  Figure 9 Comparative analysis of the the accuracy (i.e. RMSE of water depths relative to the benchmark solution at t = 110 h) and computational time (C in hours) of the simulations performed with the four models at different resolutions. In (a), as the results for Td and T models are very close, only the regression expression for the T model is shown.
the Riemann problem for all small cells at the interface, while the Riemann problem is only solved once for each cell of the traditional models. The improvement offered by the SG model against the traditional models becomes clearer when the accuracy reported in Fig. 8 is also taken into account. For instance, SG2 (100/5) delivers more accurate results than T (20), while the corresponding runtimes are 0.21 h and 4.41 h, translating into 21.15 times speedup. The results in Table 1 also show that the computational cost of SG2 is only marginally higher than SG1. Figure 9a shows the runtime (computational cost C) of simulations performed with each model as a function of the resolution. The computational cost is commonly expressed as a power function of resolution (i.e. C ∼ x −l ), where l can be used as a measure of rate of a model speed up as a result of grid coarsening. For conventional (single-grid) explicit schemes l = 3, as confirmed (l = 3.08) for the T and Td models. On the other hand, the values of l are lower for SG1 and SG2 (l = 2.4). This may be attributed to the fact that in the SG models, while the number of computational cells and time step are reduced quadratically and linearly, respectively, with increasing x (as it is the case for the T models), the number of flux computations only decreases linearly with increasing x. Figure 9b shows the RMSE of water depth predictions relative to the benchmark solution at t = 110 h (when large differences are observed between the different simulations) as a function of the model runtime. The figure provides a fair comparison of the performance of models, as it helps to identify the trade-off between accuracy and runtime of each model directly, regardless of the resolution at which the simulations were performed. According to this figure, SG2 provides the best performance (e.g. best accuracy for the same computation time), followed by SG1, Td and T.

Conclusions
This research assesses the performance of a new sub-grid model for the solution of the 2D shallow water equations, which was developed by extending the first-order model by Shamkhalchian and de Almeida (2021) to also include a second-order spatial reconstruction of flow variables. The model solves the governing equations using the Godunov finite volume method in a domain discretized by two structured nested meshes. The low-resolution mesh is the computational mesh, over which the governing equations are solved, while the high-resolution mesh provides detailed representation of inputs such as the bed topography and roughness. The proposed model is aimed at upscaling the solution from the high-to low-resolution mesh to enable computationally efficient coarse resolution solutions to be obtained at high accuracy. Four models, named SG1, SG2 (first and second order sub-grid models), T and Td (singlemesh, conventional models that reconstruct the water surface elevation and depth, respectively, and the unit width discharge) were compared via an idealized and a real-world test case. This comparison, which was mainly focused on subcritical flows, showed that SG2 provides the most accurate solutions, and that the improved accuracy is more significant at low-resolution simulations. This makes the SG2 particularly suitable for largescale flood inundation problems, where low-resolution meshes are typically required to reduce the computational cost. Our observations agree with conclusions from previous research (i.e. Begnudelli et al., 2008), and showed that at high to midresolution (e.g. x < 50 m, in our tests), the order of accuracy used to reconstruct the dependent variables is less important than the terrain model (e.g. first versus second order used previously, or the sub-grid mesh in this work). The results in this paper indicate that at coarse resolution, in addition to the terrain model, the order of the reconstruction of flow variables is key to the accuracy of the solutions. = vector of bed slope source term (m s −1 , m 2 s −2 , m 2 s −2 ) S f = vector of friction slope source term (m s −1 , m 2 s −2 , m 2 s −2 ) t = time (s or h) T p = parameters in Eqs (8) and (9), p = 1, 2, 3 (m (p−2/3) s −1 ) T p = parameters in Eqs (11) and (12)