Multidimensional hydrological–hydraulic model with variational data assimilation for river networks and floodplains
 ^{1}Laboratoire des sciences de l'ingenieur, de l'informatique et de l'imagerie (ICUBE), Fluid Mechanics Team, CNRS, Université de Strasbourg, France
 ^{2}INRAE (Irstea), Aix Marseille Université, RECOVER, AixenProvence, France
 ^{3}Institut de Mathematiques de Toulouse (IMT), Toulouse, France
 ^{4}INSA Toulouse, Toulouse, France
 ^{1}Laboratoire des sciences de l'ingenieur, de l'informatique et de l'imagerie (ICUBE), Fluid Mechanics Team, CNRS, Université de Strasbourg, France
 ^{2}INRAE (Irstea), Aix Marseille Université, RECOVER, AixenProvence, France
 ^{3}Institut de Mathematiques de Toulouse (IMT), Toulouse, France
 ^{4}INSA Toulouse, Toulouse, France
Correspondence: Léo Pujol (leob.pujol@gmail.com)
Hide author detailsCorrespondence: Léo Pujol (leob.pujol@gmail.com)
This contribution presents a novel multidimensional (multiD) hydraulic–hydrological numerical model with variational data assimilation capabilities. It allows multiscale modeling over large domains, combining in situ observations with highresolution hydrometeorology and satellite data. The multiD hydraulic model relies on the 2D shallowwater equations solved with a 1D–2D adapted single finitevolume solver. Onedimensionallike reaches are built through meshing methods that cause the 2D solver to degenerate into 1D. They are connected to 2D portions that act as local zooms, for modeling complex flow zones such as floodplains and confluences, via 1Dlike–2D interfaces. An existing parsimonious hydrological model, GR4H, is implemented and coupled to the hydraulic model. The forwardinverse multiD computational model is successfully validated on virtual and real cases of increasing complexity, including using the secondorder scheme version. Assimilating multiple observations of flow signatures leads to accurate inferences of multivariate and spatially distributed parameters among bathymetry friction, upstream and lateral hydrographs and hydrological model parameters. This notably demonstrates the possibility for information feedback towards upstream hydrological catchments, that is, backward hydrology. A 1Dlike model of part of the Garonne River is built and accurately reproduces flow lines and propagations of a 2D reference model. A multiD model of the complex Adour basin network, with inflow from the semidistributed hydrological model, is built. Highresolution flow simulations are obtained on a large domain, including fine zooms on floodplains, with a relatively low computational cost since the network contains mostly 1Dlike reaches. The current work constitutes an upgrade of the DassFlow computational platform. The adjoint of the whole tool chain is obtained by automatic code differentiation.
The accurate estimation of storage and fluxes in surface hydrology is an essential scientific question linked to major socioeconomic issues in flood and drought forecasting, particularly with regards to the ongoing climate change and potential intensification of the water cycle and hydrological hazard (Allen et al., 2019; Iturbide et al., 2020). In this context, advanced numerical modeling tools are crucially needed to both perform meaningful and detailed representations of basinscale hydrological processes and provide sensible local forecasts. The quantities of interest range from discharge hydrographs on upstream ungauged parts of the drainage network to their translation into flow depth, velocities and submersion times on downstream floodplains. This information is difficult to access, especially for floods over large territories. Indeed, given the complexity of physical processes involved, their limited observability and the resulting hydrological responses, hydrological modeling remains a hard task and internal state fluxes are generally tinged with uncertainties (Beven, 1993; Schuite et al., 2019; Milly, 1994). Moreover, the accuracy of highresolution hydraulic computations may still be affected by complex dynamics with wet–dry fronts, multiscale and uncertain topography structures and flow model parameters (e.g., friction), uncertain quantities at open boundaries (upstream inflows but also lateral ones due to sudden local runoff, downstream controls and backwater effects), internal inflows or outflows in urban areas, and large computational domains (Monnier et al., 2016). Thus, integrated hydrological–hydraulic approaches are required (e.g., Nguyen et al., 2016; Hocini et al., 2020). Such approaches are now enabled by the increasing informative richness of multisource datasets provided by highresolution hydrometeorology and satellite remote sensing in complement to in situ measurements. Nevertheless, reaching highresolution accuracy and computational efficiency for largescale applications remains a difficult challenge because of multiscale nonlinear hydrodynamic processes over large computational domains and multiple uncertainty sources.
These uncertainties could be reduced by the optimal combination of models and multisource datasets, including highresolution maps and spatially sparse in situ flow measurements but also the growing amount of earth observation data provided by new generations of satellites, drones and sensors (e.g., Biancamaria et al., 2016, 2017; Schumann and Domeneghetti, 2016, among others). Indeed, remote sensing provides very interesting cartographic observations of the variabilities of worldwide catchments characteristics (topography, soil occupation, surface moisture, snow cover, etc.), as well as an unprecedented and increasing hydraulic visibility over river networks (Garambois et al., 2017; Montazem et al., 2019; Rodríguez et al., 2020). This growing wealth of multisensed information is a key to the design and improvement of basinscale models, as shown for accurate river network 1D hydraulic modeling enabled by recent multisource altimetric and optical satellite data in Pujol et al. (2020) and in Malou et al. (2021) (see also references therein) or accurate 2D local floodplain models with radarsensed flooding extent (Hostache et al., 2010).
Still, the underlying issue of the adequate combination of numerical models and multisource data that is heterogeneous in nature and of varied spatiotemporal resolutions remains. In order to exploit this wealth of hydrological and hydraulic information, the complexity of integrated models has to be adapted to observation types and patterns, while retaining a good representation of spatial and temporal river network variabilities and relatively low computational costs over entire catchments, i.e., large computational domains. The challenging search for consistency between such datasets and models, generally involves the calibration of highdimensional multivariate parameter vectors, and corresponding inverse problems are often illposed or even illconditioned – see discussions about “equifinality” (Bertalanffy, 1968; Beven, 2000) in hydraulic inverse problems in Garambois and Monnier (2015), Larnier et al. (2020), Garambois et al. (2020), and Pujol et al. (2020) for river reach parameter inference from satellite altimetric data.
Cascades of 2D hydrological–hydraulic models have been proposed in recent literature, for inundation mapping at large scales using a worldwide digital elevation model (DEM, e.g., Grimaldi et al., 2018; Fleischmann et al., 2020; Uhe et al., 2020) with simplified noninertial hydraulic modeling) or at a finer scale, e.g., at catchment scale for flash floods in Nguyen et al. (2016) and Hocini et al. (2020). In those studies, conceptual hydrological models of upstreamlateral subcatchments are used as the inflow for hydraulic models of river network and floodplains in a weak coupling approach, mostly performed via external coupling of numerical models. In Grimaldi et al. (2018), Fleischmann et al. (2020) and Uhe et al. (2020), a simple 2D storage cell inundation model obtained from a 1D noninertial model (Bates et al., 2010, following Hunter et al., 2008, implemented in the LISFLOODFP model) enables rasterbased inundation modeling over very large domains at relatively low computational cost (see also Fleischmann et al., 2020, for coupling of this noninertial model with the largescale hydrological model MGB; Collischonn et al., 2007, Pontes et al., 2017). In Hocini et al. (2020), an original 2D hydraulic modeling approach is used to compute steady inundation maps of various return periods at high resolution (5 m) for river networks and floodplains at catchment scale of several thousands of square kilometers (up to 5050 km^{2}). It uses “precipiton” for the resolution of the full shallowwater model, which consists in propagating elementary water volumes on the water surface, as proposed by Davy et al. (2017). In Nguyen et al. (2016), an unsteady full 2D shallowwater model (Sanders et al., 2010) is applied at relatively high resolution (10 or 30 m) in the river network and floodplains on a 808 km^{2} catchment. Note that sequential data assimilation methods based on the Kalmann filter have been carried out extensively for monovariate data assimilation with such models (see, e.g., Brêda et al., 2019, with simplified hydraulics in a satellite observability context; references therein and Table 1) at varying spatiotemporal resolutions. Current model development strives to propose combinations of highresolution accuracy and fast computation times over large domains and to incorporate multisource data assimilation methods for statefluxparameter estimation.
Brunner, 1995Sanders et al., 2010Delestre et al., 2017Steinstraesser et al., 2021Guinot et al., 2018Davy et al., 2017Kirstetter et al., 2021Galland et al., 1991Goutal and Maurel, 2002Bates et al., 2013Monnier et al., 2016Brisset et al., 2018In order to combine local accuracy and computational efficiency, the association of 1D and 2D full shallowwater hydraulic models is an appropriate approach for simulating a basinscale network in a way that is both practical and adequately accurate. Methods for coupling models of different dimensions have been developed (Miglio et al., 2005a, b; Amara et al., 2004), classically using domain decomposition (Gervasio et al., 2001) or more recently by coupling 1.5D and 2D equations, so as to model local 2D zooms of overflows overlapping with a 1D domain for inbank flows, in a variational data assimilation framework in Gejadze and Monnier (2007) and Marin and Monnier (2009). An iterative coupling strategy is applied in Barthélémy et al. (2018) between a 1D Mascaret and a 2D Telemac operational numerical model, and a sequential data assimilation technique is performed for correcting water level forecasting. A summary of some established 1D and 2D numerical hydraulic models, external coupling methods and optimization–assimilation methods is presented in Table 1. Note that largescale modeling of river routing–flooding is still generally performed with simplified models and/or external couplings, while enhancing the computational efficiency and realism of multidimensional models is important in the field of hydrological–hydraulic modeling. Furthermore, it is essential to investigate stateparameter estimation, from multisource datasets, for improving the accuracy of hydrodynamic variables estimated with such models. Variational data assimilation (VDA) is an adequate method to tackle such highdimensional and multivariate hydrodynamic inverse problems from heterogeneous datasets (Lai and Monnier, 2009; Monnier et al., 2016; Brisset et al., 2018; Larnier et al., 2020; Garambois et al., 2020; Pujol et al., 2020; Malou et al., 2021). One can identify DassFlow2D (Monnier et al., 2016) as the only 2D full hydraulic model with a secondorder solver with accurate wet–dry front treatment, parallel computation and adjointbased variational data assimilation capabilities.
This work focuses on the following scientific problems: (i) the effective representation of river network flows through multiD basinscale hydraulic–hydrological numerical modeling, (ii) the inference from multisource data of spatiotemporally distributed river channel parameters and inflows through a VDA approach, and (iii) the upstream informational feedback of river network observations to integrated hydrological models.
The present study details upgrades to the DassFlow variational data assimilation framework (DassFlow2DV2, see Monnier et al., 2016) in the form of (i) a new multiD hydraulic computational code relying on an internal mesh coupling and a single finitevolume solver and (ii) the integration and differentiation of a GR4 hydrological module, (iii) with the implementation of a regularization strategy based on Larnier et al. (2020), leading to DassFlow2DV3.
One strength of the DassFlow platform is its variational data assimilation (VDA) algorithm. It can solve highdimensional optimization problems with descent algorithms using gradients computed with the adjoint model. The adjoint model is obtained via the automatic differentiation tool Tapenade (Hascoet and Pascual, 2013). This allows the simultaneous inference of numerous model parameters, which may be needed over large domains. Let us remark that this is not the case for stochastic methods, where computational cost can increase dramatically with the number of sought parameters. The variational method also enables tackling multivariate data assimilation problems, i.e., with composite control vectors (bathymetry, friction, boundary conditions, hydrological parameters), given multisource datasets, heterogeneous in nature, and spatiotemporal resolutions (see, e.g., Brisset et al., 2018, Larnier et al., 2020, Pujol et al., 2020, and Malou et al., 2021).
To apply these unique capabilities at the scale of a river network, there is a need to develop a model with both sufficient complexity and a reasonably low computational cost. In this work, we also focus on tool simplicity, both for direct and inverse modeling, by proposing this model as a single hydraulic–hydrological tool and/or code. The proposed multiD hydraulic code consists in a single finitevolume solver applied to a 2D river network. The network is discretized into “1Dlike” reaches connected to highresolution 2D meshes in a single formulation of the shallowwater equations (SWEs). The resulting product allows building large 1Dlike river networks, connected to fine local zooms. Compared to a full 2D model, this allows computing large networks with low computational costs while maintaining the possibility to model dynamics locally at a fine scale.
The complete hydraulic–hydrological modeling of river network is achieved by coupling this hydraulic domain with a wellestablished conceptual hydrological model (GR4H “state–space” Santos et al., 2018) in a semidistributed setup. The resulting integrated tool chain has been differentiated and via the VDA algorithm enables information feedback within the whole computational domain (basin) and especially from downstream to upstream.
The source code and synthetic cases presented in this study are available online (https://doi.org/10.5281/zenodo.6342723, Pujol et al., 2022), and the current updated version of the code is available upon simple request (http://www.math.univtoulouse.fr/DassFlow, last access: 27 June 2022).
The remainder of this article is organized as follows. In Sect. 2, the modeling hypothesis, the computational resolution and inverse methods are detailed. In Sect. 3, the multiD coupling scheme is validated on a series of academic cases and several academic and reallike inference setups are investigated. The study is concluded in Sect. 4, which also outlines potential applications and improvement perspectives that the proposed method and findings bring.
This section presents the integrated and multidimensional hydrological–hydraulic model and the data assimilation approach. The model is designed for simulating spatiotemporal flow variabilities over an entire river network, from upstream hydrological responses to complex flow zones (confluences, multichannel portions, floodplains, etc.).
The modeling approach which is detailed below, is based on the following ingredients:

an integrated multiD hydraulic model – the 2D SWE with finitevolume solvers from Monnier et al. (2016) is applied to 1Dlike–2D composite meshes of river networks using a numerical flux splitting method and an effective friction power law depending on flow depth;

a numerically coupled hydrological model, the widely used GR4 model from Perrin et al. (2003) in its state–space version (Santos et al., 2018), for the sake of model differentiability;

a computational inverse method based on VDA algorithms from Monnier et al. (2016), Brisset et al. (2018) and Larnier et al. (2020) enabling spatially distributed calibration and variational data assimilation with the whole chain.
2.1 MultiD hydraulic–hydrological modeling principle
The flow model consists in a spatially distributed modeling of hydrological responses coupled to, i.e., inflowing, a seamless multiscale 1Dlike–2D hydraulic model of a river network within a 2D river basin domain denoted Ω⊂ℝ^{2} . The core idea of this work is to numerically solve the 2D shallowwater hydraulic model on a multiD discretization 𝒟_{hy}⊂ℝ^{2} of the computational hydraulic domain Ω_{hy}⊂Ω . This discretization (mesh) 𝒟_{hy} is composed of N mixed unstructured triangular and/or quadrangular cells with internal interfaces between 1Dlike and 2D zones (see Sect. 2.2.3).
Let the basin domain Ω be composed of a hydrological domain Ω_{rr} connected to a hydraulic domain Ω_{hy} (Fig. 1) such that $\left({\mathrm{\Omega}}_{\mathrm{rr}}\cup {\mathrm{\Omega}}_{\mathrm{hy}}\right)\subseteq \mathrm{\Omega}$, their borders being denoted Γ_{rr} and Γ_{hy} and the hydrological–hydraulic coupling interface between Ω_{rr} and Ω_{hy} being denoted Γ_{rh} such that Γ_{rh}⊂Γ_{hy}. In other words, the border Γ_{hy} of the hydraulic domain Ω_{hy} contains interfaces with the border Γ_{rr} of the hydrological domain Ω_{rr}. Let us also consider 𝒟_{rr}, a 2D discretization of the hydrological domain Ω_{rr}.
Finally, an unstructured lattice $\mathcal{D}={\mathcal{D}}_{\mathrm{rr}}\cup {\mathcal{D}}_{\mathrm{hy}}$ covers the basin domain Ω, and a hydrological model and a hydraulic model will be resolved, respectively, on the hydrological and hydraulic subdomains and coupled via hydrological runoff flowing into the hydraulic model at the interface Γ_{rh}, as detailed in this section along with the inverse algorithm applied to the whole hydrological–hydraulic numerical chain.
2.2 Hydraulic module
Numerical hydraulic models describing openchannel flows generally rely on the resolution of crosssectionally or depthintegrated flow equations: respectively, the 1D SaintVenant equations or 2D SWE (see, e.g., Guinot, 2010). While 1D hydraulic models enable a physically sound representation of river flows variabilities in terms of wetted Section A and discharge Q for instance, 2D hydraulic models in flow depth h and depthintegrated velocity $\mathit{u}={\left(u,v\right)}^{T}$ enable us to tackle more complex flow zones such as confluences and/or diffluences and floodplain flows. The 2D shallowwater model used in the proposed approach is presented here with the adaptation of the finitevolume solver from Monnier et al. (2016) for multiD modeling. Note that 1D SaintVenant equations are presented in Appendix B along with their resolution method in DassFlow1D (Brisset et al., 2018; Larnier et al., 2020), which is used for comparison in this study.
2.2.1 Mathematical flow model
On the hydraulic computational domain Ω_{hy}⊂ℝ^{2} and for a time interval [0,T], the 2D SWEs with the Manning–Strickler friction term, in their conservative form, are written as follows:
with h being the water depth [m] and $\mathit{u}={\left(\mathit{u},v\right)}^{T}$, the depthaveraged velocity [m s^{−1}], being the flow state variables. g is the gravity magnitude [m s^{−2}], b the bed elevation [m] and n the Manning–Strickler friction coefficient [s m${}^{\mathrm{1}/\mathrm{3}}$], i.e., the flow model parameters. U is the vector of state variables and F(U) (G(U)) is its flux over the x (y) direction. S_{g} is a gravitational term and S_{f} is a friction term. Classical initial and boundary conditions adapted to real cases are considered (see Monnier et al., 2016, 2019, for details).
An effective friction law consisting in a simple power law n=αh^{β} is introduced, as previously done for 1D SWE for effective modeling with simplified multichannel river geometry in Garambois et al. (2017) and Brisset et al. (2018).
2.2.2 Building up equivalencies between 2D and 1D flow states
Onedimensionallike reaches refer to river reaches where the following meshing strategy has been applied: quadrangular cells are built such that their interfaces are perpendicular to the main flow direction and span the whole river (“bankfull”) width as a traditional 1D crosssection (XS) would. This leads to a series of quadrangular cells, each linked to a single upstream and downstream cell. The 1Dlike approach implicitly assumes a rectangular XS shape which potentially impacts the representation of (i) a section hydraulic geometry (Leopold and Maddock, 1953) and (ii) longitudinal hydraulic controls and flow variabilities.
With a view to putting the multiD model in coherence with real flow physics, a continuity condition between 1D and 1Dlike models states and parameters is required. This continuity condition is enforced at a section, in a prismatic channel such that the uniform permanent flows, that is, equilibrium, are preserved.
Let us consider a reference 1D model in (A,Q) variables with the bankfull width value W_{1D}. The friction term reads ${S}_{\mathrm{f},\mathrm{1}D}=\frac{{n}^{\mathrm{2}}Q\leftQ\right}{{A}^{\mathrm{2}}{R}_{h}^{\mathrm{4}/\mathrm{3}}}$ (Appendix B). In the corresponding 1Dlike model in $\left(h,u,v\right)$ variables (see Sect. 2.2.1), the friction term reads ${S}_{\mathrm{f},\mathrm{1}\mathrm{Dlike}}=\frac{{n}^{\mathrm{2}}\left\right\mathit{u}\left\right}{{h}^{\mathrm{1}/\mathrm{3}}}\mathit{u}$
Considering 1D flow states over an idealized river section (Fig. 2a) and the hypothesis of local flow equilibrium (uniform, steadystate) with identical wetted areas A and water surface (WS) widths W, the continuity condition implies that
where n_{1Dlike} (h_{1Dlike}) is the Manning–Strickler friction coefficient (flow depth) in the 1Dlike model i.e., the coefficient in the 2D SWE (Eq. 1).
With the additional assumption of a rectangular XS (as it will be assumed in some test cases), we have ${h}_{\mathrm{1}\mathrm{Dlike}}=A/{W}_{\mathrm{1}\mathrm{D}}$, which leads to ${n}_{\mathrm{1}\mathrm{Dlike}}={n}_{\mathrm{1}\mathrm{D}}{\left(\frac{{h}_{\mathrm{1}\mathrm{D}}}{{R}_{h,\mathrm{1}\mathrm{D}}}\right)}^{\mathrm{2}/\mathrm{3}}$.
This “1Dlike equivalent friction” leads to a perfect fit in WS elevation of a 1Dlike model and a 1D model in a straight prismatic channel at the given uniform regime (results not shown here for brevity). Figure 2b shows the evolution of the ratio $h/{R}_{h}$ vs. h. For rectangular and parabolic XS, the ratio ${n}_{\mathrm{1}\mathrm{Dlike}}/{n}_{\mathrm{1}\mathrm{D}}$ is expected to increase with h. Thus, one can naturally expect an overestimation (underestimation) of the actual friction coefficient by n_{1Dlike} at lower flows (greater flows). This would lead to an overestimation (underestimation) of the 1D WS elevation by the “equivalent” 1Dlike model. However, later it will be considered a power law in h to model the friction coefficient (Sect. 2.2.1), which provides a better fit to the 1D WS elevations outside of the considered permanent flow.
Note that longitudinal controls and flow variabilities in 1Dlike models are assessed using synthetic cases in Sect. 3.2.4.
2.2.3 Multidimensional hydraulic model
Over a given cell K∈Ω_{hy} of area m_{K}, the piecewise constant values ${\mathbf{U}}_{K}=\frac{\mathrm{1}}{{m}_{K}}{\int}_{K}\mathbf{U}\mathrm{d}K$ are approximated. Recall that the finitevolume approach applied to the homogeneous part of the hyperbolic system of Eq. (1) (that is, without the friction source term S_{f} but including a consistent discretization of the gravitational source term S_{g}) is written as follows:
where ${\mathbf{U}}_{K}^{n}$ and ${\mathbf{U}}_{K}^{n+\mathrm{1}}$ are the piecewise constant approximations of $\mathbf{U}={\left(h,hu,hv\right)}^{T}$ at time t^{n} and t^{n+1} (with ${t}^{n+\mathrm{1}}={t}^{n}+\mathrm{\Delta}{t}^{n}$), and F_{e} stands for Riemann fluxes through each edge e of the border ∂K of the cell K, with each adjacent cell Ke. The length of edge e is m_{e}, and n_{e,K} is the unit normal to e oriented from K to Ke.
The finitevolume schemes are those developed in Couderc et al. (2013) and Monnier et al. (2016). The discretization of the friction source term is described in Appendix A.
Based on the first and secondorder finitevolume solvers of Couderc et al. (2013) and Monnier et al. (2016), a “1D–2D” coupling technique is introduced following a similar concept to FinaudGuyot et al. (2018) (urban geometries and porosity context) to compute numerical fluxes on each interface between a 1Dlike quadrangular mesh cell connected to several 2D cells as schematized in Fig. 3.
At the multiD interfaces, that is, in the case of n>1 cells x_{K,i} i∈[1…n] adjacent to the same interface of another cell Ke (see Fig. 3 for notations and Sect. 3.3.2 for a reallike example), a special treatment is applied. It consists in the Riemann fluxes being calculated for each cell K_{i} using the state from the same corresponding Ke cell over an interface of length ${m}_{{\mathrm{e}}_{\mathrm{i}}}$. In the end, the flux crossing the interface $e=\cup {e}_{\mathrm{i}}$ is equal to the sum of the fluxes crossing the e_{i} interfaces: ${\mathbf{F}}_{\mathrm{e}}={\sum}_{\mathrm{i}=\mathrm{1}\mathrm{\dots}\mathrm{n}}{m}_{{\mathrm{e}}_{\mathrm{i}}}{\mathbf{F}}_{{\mathrm{e}}_{\mathrm{i}}}\left({\mathbf{U}}_{{K}_{\mathrm{i}}}^{n},{\mathbf{U}}_{Ke}^{n},{\mathbf{n}}_{{\mathrm{e}}_{\mathrm{i}},K}\right)$.
This type of internal interface has been implemented in the numerical solvers from Monnier et al. (2016) in the DassFlow platform, which includes a solver with secondorder accuracy in space. This solver, developed in Couderc et al. (2013) and Monnier et al. (2016), is accurate and robust for wet–dry front propagations and fully applies in the present context. Note that the lateral distribution of variables across the 1D–2D interface is not constrained. The source code and synthetic cases are available upon simple request.
2.3 Hydrological module
In order to simulate the hydrological response of upstream and lateral subcatchments flowing into the river network within a river basin, a hydrological module is coupled to the 2D SW hydraulic model, in a semidistributed manner for simplicity here.
Over the hydrological domain Ω_{rr}, we consider the discretization 𝒟_{rr} here into C subcatchments $\left\{{\mathrm{bv}}_{\mathrm{1}},\mathrm{\dots},{\mathrm{bv}}_{C}\right\}$ and corresponding disjointed subdomains ${\mathrm{\Omega}}_{\mathrm{rr},\mathrm{i}\in \left[\mathrm{1}\mathrm{\dots}C\right]}\subset {\mathrm{\Omega}}_{\mathrm{rr}}$, the outlet coordinates of which are ${x}_{\mathrm{i}}\in {\mathrm{\Omega}}_{\mathrm{rr},\mathrm{i}}\in \left[\mathrm{1}\mathrm{\dots}C\right]$. For a given subcatchment, the hydrological model is a dynamic operator, consisting in coupled state equations (ordinary differential equations, ODEs) and output equations, that reads
where h(x,t) is the state variable vector, P(x,t) and E(x,t) is the observable rainfall and evapotranspiration inputs (here spatially averaged over the area of each subcatchment i for semidistributed modeling), Q(x,t) is the “observable” output discharge, and θ(x) is the “unobservable” parameter vector.
The widely used, parsimonious and robust conceptual hydrological model GR4 (Perrin et al., 2003), in its state–space version from Santos et al. (2018), was chosen in this study for parsimony and differentiability reasons. The original lumped hydrological numerical model has been deployed in a semidistributed manner in the DassFlow framework.
The model is composed of two nonlinear stores for production (soil moisture accounting) and routing and a Nash cascade composed of a series of linear stores replacing the unit hydrograph from Perrin et al. (2003). Being a set of ODEs with explicit dependency to parameters, this hydrological model (Eq. 4, see detailed GR4 equations in Appendix C) is differentiable. Moreover, the Fortran code is differentiable with the automatic differentiation tool Tapenade (Hascoet and Pascual, 2013).
The evolution of reservoir states and model outputs and inputs is presented for a sample rain event in Fig. 4.
The input of the hydrological model is the evapotranspiration and precipitation E_{n} and P_{n}; the output is discharge $q\left(t\right)={Q}_{r}+{Q}_{d}$ [mm h^{−1}]. E_{n} and P_{n} are classically imposed from data time series as a piecewise constant at the fixed temporal resolution (e.g., hourly). The numerical resolution is achieved with an implicit Euler algorithm with an adaptive substep algorithm enabling us to reduce numerical errors especially for high flows (Santos et al., 2018). The initial states of the stores is given by a 1year warmup run. The discharge q is injected into Ω_{hy} at a subinterface of Γ_{hyrr}, either as an upstream or lateral flow.
The calibrated parameters will be the classical four parameters ${\mathit{\theta}}_{\mathrm{rr}}={\left({c}_{k,\mathrm{i}}\right)}_{k\in \mathrm{1}\mathrm{\dots}\mathrm{4},\mathrm{i}\in \mathrm{1}\mathrm{\dots}C}$ of GR4, i.e., the reservoir capacities in millimeters for each subcatchment i∈1…C. They will constitute the parameter vector θ_{rr} considered in the forthcoming VDA experiments. Other parameters, such as several drainage law exponents or the number of Nash cascade stores, are not optimized as classically done with the GR4 model (Perrin et al., 2003). They are set at values from Santos et al. (2018).
2.4 Inverse algorithm: variational data assimilation
Given spatiotemporal flow observables, provided by in situ and airborne sensors for instance, the inverse algorithm consisting in variational data assimilation (VDA) aims at estimating the unknown or uncertain “input parameters” of the hydrological–hydraulic chain composed of a hydraulic model, presented in Sect. 2.2, and a hydrological model, presented in Sect. 2.3. Detailed knowhow on VDA may be found in online courses (see, e.g., Bouttier and Courtier, 2002, and Monnier, 2014). This group of VDA algorithms relies on cost gradients, i.e., variations, here computed by the adjoint model, to minimize a misfit to the observed reality (see Fig. 5). The adjoint model of the whole toolchain composed of DassFlow multiD with GR4 modules is obtained by automatic differentiation with the Tapenade engine (Hascoet and Pascual, 2013). Note that 2D, 1Dlike and 1D–2D models are functionally identical from the point of view of the adjoint code and the assimilation process.
We consider the following set of spatiotemporal observations of water surface and discharge over the river domain Ω⊂ℝ^{2}:
with Z_{o} the observed WS elevation [m] above reference elevation, Q_{o} the observed discharge [m^{3} s^{−1}], N_{o,Z} the number of altimetric observation points and N_{o,Q} the number of observed discharges over Ω.
Note that observed discharge Q_{o} may be a value of a hydraulic discharge at a flow XS in Ω_{hy} or within the hydrological domain Ω_{rr} and especially at the outlet of a subcatchment here. Also note that other water surface observables could be considered, such as water surface velocity observations or dynamic water masks, but this is beyond the scope of this study.
Given river stage and/or discharge observations, the aim here is to estimate unknown or uncertain quantities of the hydrological–hydraulic model among discharge hydrographs ${Q}_{\mathrm{i}}\left(t\right),\phantom{\rule{0.125em}{0ex}}i\in \left[\mathrm{1},N\right]$ on the border of the hydraulic domain (in 1Dlike and/or 2D parts), spatially distributed hydraulic parameters over a river network (bathymetry elevation b or friction n) and hydrological parameters ${\mathit{\theta}}_{\mathrm{rr}}={\left({c}_{k,\mathrm{i}}\right)}_{k\in \mathrm{1}\mathrm{\dots}\mathrm{4},\mathrm{i}\in \mathrm{1}\mathrm{\dots}C}$ for all subcatchments ${\cup}_{\mathrm{i}=\mathrm{1}}^{C}{\mathrm{\Omega}}_{\mathrm{rr},\mathrm{i}}={\mathrm{\Omega}}_{\mathrm{rr}}$.
The control vector containing the sought quantities is denoted θ in what follows:
with θ_{hy} and θ_{rr} the control vectors of, respectively, the hydraulic and the hydrological modules, N the number of inflows control points in space and T the number of inflow values in time, M the number of bathymetryfriction control cells, and P the number of hydrological units.
We consider a cost function j_{obs} aiming at measuring the discrepancy between simulated and observed flow quantities on the hydraulic–hydrological computational domain Ω. This cost function is defined as
This cost function contains either misfit to WS elevation, ${j}_{Z}\left(\mathit{\theta}\right)=\frac{\mathrm{1}}{\mathrm{2}}{\u2225{Z}_{\mathrm{o}}\left(t\right)Z\left(\mathit{\theta},t\right)\u2225}_{{\mathcal{O}}_{\mathcal{Z}}}^{\mathrm{2}}$, or misfit to discharge, ${j}_{Q}\left(\mathit{\theta}\right)=\frac{\mathrm{1}}{\mathrm{2}}{\u2225{Q}_{\mathrm{o}}\left(t\right)Q\left(\mathit{\theta},t\right)\u2225}_{{\mathcal{O}}_{\mathcal{Q}}}^{\mathrm{2}}$.
The metrics (symmetric positive definite matrices) 𝒪_{𝒵} and 𝒪_{𝒬} are based on the inverse of the observation error covariances. This enables us to regularize the inverse problem; see, e.g., Bouttier and Courtier, 2002, Monnier, 2021, Asch et al., 2016, and references therein for related discussions.
Moreover, we classically enrich the cost function with a regularization term: $j\left(\mathit{\theta}\right)={j}_{\mathrm{obs}}\left(\mathit{\theta}\right)+\mathit{\gamma}{j}_{\mathrm{reg}}\left(\mathit{\theta}\right)$ with j_{reg} a Tikhonovtype regularization term. Here we consider a regularization on the bathymetry only: ${j}_{\mathrm{reg}}\left(\mathit{\theta}\right)=\frac{\mathrm{1}}{\mathrm{2}}{\sum}_{\mathrm{i}=\mathrm{1}}^{M}\left({\left({\partial}_{x}{b}_{\mathrm{i}}\right)}^{\mathrm{2}}+{\left({\partial}_{y}{b}_{\mathrm{i}}\right)}^{\mathrm{2}}\right)$ with θ defined by Eq. (6).
The regularization term adds convexity to the cost function. Moreover, it here dampens the bathymetry highest frequencies. Recall that low Froude flows; i.e., subcritical flows naturally act as a lowpass filtering of the bathymetry shape; see Martin and Monnier (2015) and Gudmundsson (2003).
The total cost function j is minimized starting from a background value θ^{(0)}. Following Lorenc et al. (2000) (see also Larnier et al., 2020), the following change in variables is applied:
with B being the covariance matrix of the background error.
Then by setting J(k)=j(θ), the optimization problem which is solved is actually the following:
The firstorder optimality condition of this optimization problem (Eq. 9) reads ${\mathbf{B}}^{\mathrm{1}/\mathrm{2}}\mathrm{\nabla}j\left(\mathit{\theta}\right)=\mathrm{0}$. The change in variables based on the covariance matrix B acts as a preconditioning of the optimization problem. This optimization problem is solved using a firstorder gradientbased algorithm, more precisely the classical limitedmemory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) quasiNewton algorithm (Zhu et al., 1997) or, in some cases in this study, its bounded version LBFGSB (Zhu et al., 1997) but without variable change (Eq. 8).
The choice of the covariance matrix B, represents an important a priori information and greatly influences the computed solution of the inverse problem. Assuming the unknown parameters are independent variables, the matrix B is defined as a block diagonal matrix:
Each block matrix of ${\mathbf{B}}_{{\mathrm{\Omega}}_{\mathrm{hy}}}$ is defined as a covariance matrix (positive definite matrix) using 2D kernels for spatial controls. Here following Larnier et al. (2020), we consider
The parameters L_{Q} and (L_{b},L_{n}) act as correlation lengths. These parameters are usually empirically defined. However the expression of B_{b} and the correlation lengths can be derived from physically based estimations following Malou and Monnier (2021).
The following stopping criteria are used to stop the iterative optimization process if (i) the cost function does not decrease over a set number of iterations or (ii) the current cost gradient, normalized by the initial cost gradient, goes under a set objective value. The multiD hydrological–hydraulic tool chain presented above has been implemented in the latest version of Monnier et al. (2019).
3.1 Numerical experiments design
Both synthetic and real cases are considered to test the forward and inverse modeling capabilities of the proposed computational chain for river network and floodplain simulation.
First, the multiD hydraulic model (Sect. 2.2.3) is validated against reference hydraulic models on synthetic cases corresponding to typical hydraulic complexities: (i) simple straight channel, (ii) confluence, and (iii) straight, rectangular and parabolic, channels with effective parameterization of friction and bathymetry. For fluvial regimes in the context of altimetry, these hydraulic complexities generate hydraulic controls. Following Montazem et al. (2019), we define hydraulic controls as characterized by a maximal deviation of the water depth from the normal depth (see, e.g., Chow, 1959, and Dingman, 2009, for definition) at the reach scale.
A series of inference cases are considered in twin experiment setups. In a twin experiment, a reference model acts as a synthetic truth and is used to generate observations of model variables (e.g., b+h in Ω_{hy} or Q in Ω_{rr}). No observation noise is considered in this study. The VDA method is then applied with these observations on an altered version of the reference model. Inferences of temporal forcings, channel parameters and hydrological parameters are carried out.
3.2 Synthetic cases
First, the proposed multiD hydraulic solver (Sect. 2.2.3) is evaluated and compared to a fine 2D reference model on two simple configurations, a straight channel and a confluence, that feature frontal interfaces between 1Dlike and 2D meshes. Next, to investigate the reproducibility of hydraulic controls using a 1Dlike meshing approach and effective modeling, this approach is compared to a 1D reference model in three typical channel configurations. Finally, inferences of inflow hydrographs ${Q}_{\mathrm{i}}\left(t\right),\phantom{\rule{0.125em}{0ex}}i\in \left[\mathrm{1}\mathrm{\dots}N\right]$ and of a friction power law n=αh^{β} are carried out using WS observables.
3.2.1 Forward multiD hydraulic cases
Straight channel case
A prismatic rectangular channel and a multiD mesh are considered (1Dlike to 2D to 1Dlike, see Fig. 6). The channel width is 300 m and its length is 2300 m. A rating curve is imposed downstream. This multiD model is compared to a reference 2D model with a refined mesh (mesh not shown, 2400 cells, average edge length ∽25 m). The multiD waterline is validated against the 2D model at permanent flow $\left(Q=\mathrm{100}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}\right)$ and the modeled downstream discharges are compared for a flood hydrograph (Fig. 6b). Both first (not shown) and secondorder numerical solvers allow a close fit to the target water line (Fig. 6a, bottom diagram).
At permanent flow, a slight misfit is observed between the 2D and multiD WS elevations with the secondorder scheme (Fig. 6a, top diagram, relative misfit <0.15 % at 1000 m). This is due to the approximation error in the multiD model caused by large spatial steps in 1Dlike reaches (dx=200 m). Indeed, this misfit is reduced by reducing the 1Dlike cell length (not shown). This is confirmed and showcased in the next subsection.
At the interface between 1Dlike and 2D meshes, a slight jump in WS elevations can be observed at all 2D cells (Fig. 6a, top diagram). This is due to the secondorder scheme, which is currently not designed for multiD interfaces. Recall that no constraint is imposed on the lateral distribution of computed variables. The technical implementation of this reconstruction for 1D–2D interfaces will be done in the next version of DassFlow.
During a varied flow event, the outflow of both models is close to identical (Fig. 6a, middle diagram). The flow is correctly transmitted at multiD interfaces and, at this scale, the 1Dlike meshes are adequate to model a flood wave propagation.
Simple confluence case
A simple symmetrical confluence is modeled using a multiD mesh. The channel width is 300 m in the downstream reach and 150 m in the upstream reaches. Two inflow hydrographs are imposed at the two upstream interfaces. The maximum abscissa of the mesh points are 0 and 2075 m. The 2D part (average edge length ∽5 m) contains a confluence flow zone (Fig. 6), while the 1Dlike mesh (dx=100 m) covers the upstream and downstream reaches. At permanent flow, the multiD model compares well with the reference 2D fine model (mesh not shown, 15 000 cells, average edge length ∽5 m) and leads to similar conclusions to the above paragraphs. Using a shorter spatial step in the 1Dlike reaches (10 m) reduces the difference between reference and multiD model and allows for a nearly perfect fit to the reference WS elevation at permanent flow (Fig. 6a, middle diagram). During a varied flow event, the outflows modeled with the multiD model are very close to reference ones: slight differences during rising and falling limbs, same peak time, and Nash–Sutcliffe efficiency (NSE) =0.996.
Hydraulic controls and effective friction
To investigate the reproducibility of hydraulic controls for fluvial flows, which are characterized by a maximal deviation of the water depth from the normal depth at the reach scale (see definition in Montazem et al., 2019), using a 1Dlike approach, a set of typical channel variabilities and hydraulic controls are considered. Let us compare 1D and 1Dlike waterlines in a series of synthetic cases: (i) a straight rectangular channel, (ii) a straight rectangular channel with a midchannel slope break and (iii) a straight parabolic channel.
Direct calibration
Recall the equivalent friction formulation (Eq. 2) designed to match the water line at a section and at equilibrium. For each 1Dlike case, waterlines with 1D friction ($n=\mathrm{0.05}\phantom{\rule{0.125em}{0ex}}\mathrm{s}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}/\mathrm{3}}$) and effective equivalent friction are generated. Effective friction values are calculated using Eq. (2) and results from the corresponding 1D case. Reference 1D flow lines are computed with DassFlow solving 1D SaintVenant equations in (A,Q) state variables with a Preissmann scheme (see Appendix B and also Larnier et al., 2020).
In a straight rectangular prismatic channel, with the 1D normal water depth imposed downstream (Fig. 8a), a backwater curve is computed by the DassFlow1D model (in red). With a 1D homogeneous friction ($n=\mathrm{0.05}\phantom{\rule{0.125em}{0ex}}\mathrm{s}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}/\mathrm{3}}$), the 1Dlike approach yields an underestimated waterline (in blue). A homogeneous effective friction allows us to correct the underestimation and to better match the 1D normal water depth over the whole domain (in green). The remaining misfit can be attributed to numerical errors, especially numerical diffusion from the Preissmann scheme (1D model spatial step is dx_{1D}=100 m).
A first complexification of this case consists in the introduction of a local hydraulic control point in the form of a slope break at x=10 km (Fig. 8b). In this setup, both 1D and 1Dlike models generate M2 backwater curves; see, e.g., Dingman (2009). The hydraulic control generated at the slope break is well represented with a 1Dlike approach, given the aforementioned numerical errors due to the coarse grid.
Another complexification of the first case consists in changing from a rectangular XS to a parabolic one (Fig. 2). In this case, both equivalent friction – a single homogeneous patch – and effective bathymetry, in the form of a homogeneous shift δ_{b} of the reference bathymetry, are needed to match the 1D WS elevation. Equivalent friction only allows us to model identical wetted sections at permanent bankfull flow for a given channel width (Sect. 2.2.2). In this case, matching the 1D wetted section does not equate to matching its WS elevation; thus an effective bathymetry is used.
According to the above model comparisons, effective parameterization of channel parameters is sufficient to reproduce 1D behaviors with a 1Dlike approach for permanent bankfull flows in simple geometries. Outside of the permanent bankfull flow, a friction power law n=αh^{β} can be used (Sect. 2.2.1). This friction aims at better 1Dlike representativity when modeled wetted surfaces (and other XS variables) are different in 1D and 1Dlike models for the same WS elevation.
Inverse calibration
In the following comparison, powerlaw parameters α and β are obtained using VDA in a twin experiment setup. Observations are taken at two stations (Fig. 9). Prior values are α=0.05 s m${}^{\mathrm{1}/\mathrm{3}}$ and β=0.
The classic friction law and a power law with calibrated parameters are compared in Fig. 9 during a varied flow event (assimilation time window of 270 d). A range of water depth higher than the bankfull depth used to calculate equivalent friction parameters is simulated. This results in an increase in misfit to the 1D WS elevation at high flow (Fig. 9, in blue and red). The calibrated friction power law (in orange) somewhat reduces the misfit during high flows in this simple synthetic channel and within the simulated water depth range.
3.2.2 MultiD hydraulic–hydrological data assimilation
Temporal forcing inference
Simple inference tests are carried on the confluence case from Sect. 3.2.1, following Pujol et al. (2020), in a twin experiment setup. A control vector $\mathit{c}=\left({Q}_{\mathrm{1}}\left(t\right),{Q}_{\mathrm{2}}\left(t\right)\right)$ is considered, where Q_{1} and Q_{2} are sinusoidal inflow hydrographs injected at the upstream cell of the two upstream reaches. Pujol et al. (2020) shows that the minimal requirement to infer multiple spatially distributed temporal forcings simultaneously is to observe either (i) both of their “unmixed” signatures, at one point each, or (ii) one of their unmixed signatures at one point and their mixed signatures at a second point. We verify this for configuration (ii). Unmixed signature observations are generated at a 1Dlike cell of the upstream reaches; mixed signature observations are generated at a cell in the 2D part (Fig. 10a, red crosses). Prior values for the inflow hydrographs are constant and set to their average value $\left({Q}_{\mathrm{1}}\left(t\right)={Q}_{\mathrm{2}}\left(t\right)=\mathrm{100}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}\right)$. Using sets of two observation points over a 5 h period (see Fig. 10), we are able to reproduce the results of Pujol et al. (2020) but with the proposed multiD hydraulic model: the VDA algorithm enables us to infer Q_{1} and Q_{2} close to perfectly from configurations (i) and (ii) (Fig. 10b).
Hydrological parameter inference
In this second twin experiment, the issue of the spatially distributed calibration of a hydrological model is studied, from multisource observations of the river network. The confluence case above is used and this time, the upstream inflows are generated by the GR4H module applied to two upstream catchments flowing into the hydraulic module. For each catchment, synthetic rain and evaporation time series and a hydrological parameter set ${\mathit{\theta}}_{\mathrm{rr}}=\left({c}_{\mathrm{1}},\mathrm{\dots},{c}_{\mathrm{4}}\right)$ are used to generate a discharge time series over a 1year period (not shown). Mixed observables are used: WS elevation is observed at a single downstream point, and hydrological model discharge is observed at the outlet of one of the catchments, over 360 d. They provide the same signal observability – mixed and unmixed signals – as in the above paragraph but with observations of a different nature. Since the observed variables are of a different nature and amplitude, we introduce a normalization. The cost function is here ${j}_{\mathrm{obs}}={j}_{Z}+{j}_{Q}$, with each term being normalized by the number of observations and by their range of variation such that
and
Those correspond to two separate normalized squared RMSE with N_{o,Q} and N_{o,Z} denoting the number of observation stations and T_{Q} and T_{Z} the number of observation time steps.
The control vector $\mathit{c}=({\mathit{\theta}}_{\mathrm{rr}\mathrm{1}},{\mathit{\theta}}_{\mathrm{rr}\mathrm{2}})$ contains the two sets of four hydrological parameters each. For this synthetic case, an inequality constraint of the control parameters is imposed with the bounded LBFGSB algorithm (Zhu et al., 1997). Indeed, restricted research intervals are considered for the three first parameters of each catchment, namely a 5 % bracket around their target values used as a prior, while c_{4} is sought in its expected variation range from an erroneous prior (Table 2). Expected ranges for GR4H parameters are provided in Le Lay (2006), for the classical GR4 formulation.
In this setup, both j_{Q} and j_{Z} are reduced during the iterative steps and the inferred value for both c_{4} parameters are very close to the target value starting from an erroneous prior (Fig. 11b). Bounded parameters c_{1} and c_{2} vary slightly between their bounds, while c_{3} is locked at its lower bound from the first iteration.
3.3 Real cases
3.3.1 Garonne River: 1Dlike effective model of a 2D reference case
A 1Dlike model of a reach of the Garonne River, in southern France, is built from real data and calibrated here via assimilation of spatially distributed WS observations in a twin experiment setup. A full 2D model of 75 km of the Garonne River (not shown, 867 500 cells including floodplains, average edge length ∽25 m) is used to generate reallike observables on this wellknown case (Garambois and Monnier, 2015; Brisset et al., 2018; Larnier et al., 2020; Monnier et al., 2016). First, river extent is derived from the 2D model output at bankfull flow Q=400 m^{3} s^{−1} with a homogeneous friction n=0.05 s m${}^{\mathrm{1}/\mathrm{3}}$. This extent is in turn used to build a 1Dlike mesh, over the whole reach: single quadrangular cells cover the whole river width and are linked sequentially along the river reach (Fig. 12). Onedimensionallike cell interfaces are perpendicular to the flow direction, as would be cross sections (XSs) in a 1D model. Each cell is about 100 m in longitudinal length. This mesh was generated using the SMS.^{1} Cell bathymetry is first set using the lowest bathymetry point at each corresponding 2D XS.
Permanent bankfull flow calibration
A first expectedly imperfect 1Dlike model (Model A1) is built using the 1Dlike coarse mesh and expectedly underestimated bathymetry elevation and a homogeneous friction parameter of 0.05 s m${}^{\mathrm{1}/\mathrm{3}}$. The simulated steady WS elevation at bankfull flow is lower than that of the 2D model (average misfit of 0.858 m), which is expected since the 1D bathymetry is that of the lowest point of the 2D XS (Fig. 2). Furthermore, the 1Dlike friction is underestimated and leads to an underestimated simulated water depth (and flow surface). The WS elevation misfit does not seem to follow a significant trend from upstream to downstream, although it varies sharply at points of width variation (Fig. 13, e.g., around cell 410). Recall that the goal of this effective modeling approach is to accurately reproduce water surface signatures, including WS elevation and, tangentially, flow section (Sect. 2.2.2).
Calibration by hand. We here propose to reduce the misfit using, on the one hand, effective friction and, on the other hand, bathymetry as follows. In Model A2, the introduction of an equivalent friction parameter (Eq. 2), calculated at each 1D cell using observations from 2D XS, improves the WS elevation misfit (mean friction of n=0.062 s m${}^{\mathrm{1}/\mathrm{3}}$, average WS elevation misfit of 0.562 m). It reduces misfit overall but has no significant local influence (Fig. 13). In Model A3, we use the average WS elevation misfit from Model A2 to create a simple bathymetry shift δ_{b} that helps fit the observations. Equivalent friction parameters from Model A2 are kept and all bathymetry points are shifted by δ_{b}=0.562 m (the bottom slope is conserved). The average WS elevation misfit at bankfull flow reaches 0.003 m for Model A3, which is a very satisfying result.
Calibration by VDA. Now, the same calibration problem is addressed with an inference based on the VDA method. It is applied to the same reference model permanent flow waterline, observed at each cell. The control vector contains a single homogeneous friction parameter, as before, and spatialized bathymetry b(x) for each cell. To constrain the parameter search, two VDA processes are performed separately: a bathymetry regularization and a change in variable. Inference with bathymetry regularization leads to Model B1, with γ=1 (Sect. 2.4). Inference with change in variable (Eq. 8) leads to Model B2, with L_{b}=500 m (Eq. 11). Both models lead to average misfits close to that of Model A3: 0.0839 and 0.0844 m, respectively.
Two other inference setups based on Model B1 are considered. The number of observations is divided by 10: 72 stations, homogeneously distributed at 1 per each 1 km, are considered. In Model B1a, no regularization term is considered (γ=0). In Model B1b, a regularization is added (γ=1, chosen by trial and error). This weight can be optimally determined using iterative regularization (Malou and Monnier, 2022). It is dependent on the spatial scales of observed signals and on the discretization of inferred parameters. As presented in Fig. 13 and in Table 3, both inferences lead to low misfits of the 2D permanent WS elevation. The regularization tends to reduce extreme bathymetric variations that tend to appear far from the observation points. The iterative minimization process can be followed through the cost function value and the parameter gradients (Fig. 14). Values are normalized over their initial (iteration 0) value. For Model B1a, the optimal control is reached after 40 iteration (for a normalized cost of $\mathrm{8.6}\times {\mathrm{10}}^{\mathrm{10}}$), while for Model B1b it is reached in 85 iterations (for a normalized cost of $\mathrm{5.5}\times {\mathrm{10}}^{\mathrm{5}}$). Both models follow the same trajectory up to around iteration 12, where the regularization term j_{reg} reaches the same order of scale as the WS elevation misfit term j_{obs}.
Variational calibration of the channel parameter has allowed us to fit a permanent regime water line. The following paragraphs study the reproduction of propagation in a calibrated model during a varied flow event.
Variational calibration for a flow event
Let us now consider a varied flow event, without flooding in the 2D model. This event lasts 10 d, with a peak discharge of 702 m^{3} s^{−1} at the start of day 2, a minimum discharge of 160 m^{3} s^{−1} and an average discharge of 382 m^{3} s^{−1}. Given our previous inference attempts, we know that we can find a couple (n,b(x)) or (n,δ_{b}) that minimizes the average misfit. For the sake of simplicity, we consider the couple (n,δ_{b}). For a varied flow event, this couple would be an optimal parameterization for the average observed water line. A first inference trial is carried out for a densely observed (in space and time) varied flow. This dense observability is not meant to imitate the actual observability of real rivers but to provide sufficient constraints for the considered inverse problem, with the aim to showcase the capability of the variational toolchain and 1Dlike model to fit finescale reallike water surface elevation (WSE) variations (see hydraulic parameter inference from sparse uneven altimetric observations of river surface in Brisset et al., 2018, Larnier et al., 2020, Garambois et al., 2020, and Pujol et al., 2020).
Using Model A1 as a prior for bathymetry and friction and observations from the 2D model, we find the following optimal parameters: $\left(n=\mathrm{0.054}\phantom{\rule{0.125em}{0ex}}\mathrm{s}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}/\mathrm{3}},{\mathit{\delta}}_{\mathrm{b}}=\mathrm{0.669}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\right)$. The resulting optimal model is Model C. To allow a better fit during high and low flows, we introduce the friction power law n=αh^{β} (Sect. 2.2.1). Using observations from the 2D model during the varied flow event, we infer the couple (β,δ_{b}). The value of α is set to 0.054 s m${}^{\mathrm{1}/\mathrm{3}}$, the inferred value from Model D. Three different β priors are tested: −0.5, 0, 0.5. δ_{b} is included in the control vector to allow modeling of more varied water depths, which may be needed to reach the optimum. All three inferences lead to the very close optimal control. Their averaging, a slight bathymetric shift δ_{b}=0.099 m and $\mathit{\beta}=\mathrm{0.017}$, leads to Model D. The direct simulation results of Model D are compared to the reference model. They correspond to an average WS elevation misfit of 0.026 m at high flow and of 0.11 m at low flow. Compared to Model C (average misfit of 0.086 and 0.20 m at, respectively, high and low flow), this is an improvement. Note that NSE values for both model (Table 3) are both extremely close to 1.
The variational calibration of the global bathymetric shift δ_{b} and of a homogeneous friction value n in a 1Dlike model has allowed us to closely fit WS elevation observations of a 10 d flood wave over the reference Garonne model. The calibration of depthdependent friction, in the form of the β parameter in n=αh^{β}, has allowed an even closer fit to this reference WS elevation over the high and low flows, i.e., a better representation of the observed flood wave propagation.
3.3.2 Adour basin: multiD hydrological–hydraulic model
The whole hydrologicalmultiD hydraulic tool chain is now tested in a real context both for forward and inverse problem resolution. A real and complex basin case is now considered: the Adour River basin in the southwest of France, with a total drained area of 16 890 km^{2} at the estuary. It is probably one of the most difficult basins to model in the country because of contrasted hydrological regimes including nival effects in the south, flash floods on small ungauged catchments, complex river network morphology, anthropized floodplains and tidal effects from downstream areas.
In this section, a multiD model, composed of 1Dlike meshes and 2D zooms over the floodplain area is built from available data (Fig. 15, 2D area in green). Then, forward and inverse flow simulations on the river network are presented.
Model construction and rainfall to inundation simulation
The following model of the Adour River combines a multiD hydraulic network model and several hydrological models of subcatchments (Fig. 15). It encompasses around 180 km of river reaches and includes the Adour River from its tidal boundary downstream of Bayonne up to a gauging station around 70 km upstream and part of its main tributaries: the Nive River (around 45 km) and the Oloron and Pau rivers (around 65 km in total). The river network contains mostly singlebranch reaches, with some notable flood areas around the city of Bayonne. The WS is observed in situ at 10 points, 5 of which are used as boundary conditions (Dax, Orthez, Escos, Cambo and Convergent in Fig. 15). Out of the five remaining stations kept for data assimilation, three are located on river reaches (Peyrehorade, Urt and Villefranque; blue points in Fig. 15) and two are located in the Bayonne area (PontBlanc and Lesseps; blue points in Fig. 15).
At the four upstream points of the modeled river network, four subcatchments are modeled with four lumped hydrological models (Sect. 2.3). Their drainage areas are about 7811, 842, 2480 and 2464 km^{2}. The hourly discharge has been extracted from the HYDRO^{2} database, while the rainfall data from the radar observation reanalysis ANTILOPE J+1, which merges radar and in situ gauge observations, is provided by Météo France. The interannual temperature data are provided by the SAFRAN reanalysis by Météo France and then used to calculate the potential evapotranspiration using the Oudin formula (Oudin et al., 2005). The rainfall and potential evapotranspiration are at a spatial resolution of a 1 km^{2} square grid and are processed at an hourly time step. Spatial averages of the rainfall and potential evapotranspiration computed using the SMASH distributed hydrological modeling platform (JayAllemand et al., 2020; Colleoni et al., 2021) over the four catchments are used as inputs for the lumped GR4H model. Lumped parameter sets for the four GR4H models are simply obtained here using for each the airGR global calibration algorithm (Coron et al., 2017). In this section, inferences are carried out only for upstream inflow hydrographs and not for hydrological parameters. Indeed, the study of global calibration and regionalization issues of spatially distributed hydrological models is left for further research.
Our multiD hydraulic modeling approach (Sect. 2) is applied to this complex case as follows. First, a “1Dlikeonly” model of the whole network is built. Then, a multiD model is built based on the 1Dlikeonly model, with the addition of a 2D mesh of a floodplain (Fig. 16a).
On the 1Dlikeonly model, the goal is to analyze 1Dlike signal propagation representation at a low computational cost. The Adour 1Dlike model is built similarly to the Garonne 1Dlike model, using DEM data to determine minor bed bank line placement and build the quadrangular cells. Bathymetry comes from 25×25 m DEM data, aggregated from fine lidar data from public databases, extracted at each cell. This rough approximation is sufficient to show the potential of the DassFlow assimilation tool chain for a largescale river network. This 1Dlikeonly model contains 1409 cells. A 24 h simulation runs in around 13 s, using uncalibrated parameter estimates.
Then, a multiD model is obtained: the existing 2D mesh from a Telemac model is coupled to a 1Dlike mesh, similarly to the reference TelemacMascaret model from the regional flood forecast center SPCGAD. The 1Dlike parts of the mesh are kept identical to the 1Dlikeonly model, while the 2D part is the mesh extracted from the TelemacMascaret model (Fig. 16a) provided by SPCGAD. For hydraulic coherence, bathymetry at coupled 1Dlike cells is taken as the average bathymetry of the linked 2D cells. This 1D–2D model contains 66 982 cells in the 2D area and 1342 cells in 1Dlike reaches. Results for a flooding event (around 6 h of computation per day, on a single thread ) in the Bayonne 2D area are presented in Fig. 16b. Modeled variables appear coherent over the 1D–2D area and a highresolution flood map in coherence with flow conditions in the whole river network is obtained.
A 1 d (physical time) flooding event is computed in around 20 min (computation time) in the 1D–2D model (with 68 324 total cells, six threads in parallel). The same event leads to a 7 s computation time for the 1Dlikeonly model (with 1409 cells, same computational setup). Remember that computation time may vary depending on the number of wetted cells and on the adaptive time step calculation. This relatively low computation time and the potential to decrease it further by using more threads in parallel indicate that this multiD method is suited to operational use. The code version this work is based on was proven to be scalable (Couderc et al., 2013). Additions made to the current version should not change this, but no numerical testing has been done.
Assimilation of WS observables to infer four upstream hydrographs
To investigate the 1Dlike effective modeling on a river network, a twin experiment setup is designed to infer a large control vector from a realistic observability in the 1Dlikeonly model.
The considered control vector is composed of the four upstream hydrographs $c=\left({Q}_{\mathrm{Dax}}\left(t\right),{Q}_{\mathrm{Esc}}\left(t\right),{Q}_{\mathrm{Ort}}\left(t\right),{Q}_{\mathrm{Cam}}\left(t\right)\right)$. Observations of WS elevation are generated at Peyrehorade, Urt and Villefranque stations and additionally at a virtual station directly downstream from Orthez. These four points give theoretically sufficient observability to identify the four upstream hydrographs. Indeed, they sample the mixed and unmixed signal similarly to the academic setup in Sect. 3.2.3. The observation pattern also corresponds to a reasonable expectation of spatial observability in French river networks. Prior hydrograph values are classically set to a constant average discharge value.
Simultaneous inference of the four hydrographs is satisfying. As shown in Fig. 17, left, upstream hydrographs injected at Cambo and Orthez are inferred very accurately (NSE >0.95). This is due to their WS elevation signature being observed to be unmixed in the respective downstream reaches. The Escos hydrograph WS signature is observed at Peyrehorade, mixed with the Orthez hydrograph WS elevation signature. This leads to a partial error of signature attribution and a less accurate inference (NSE =0.57). The Dax hydrograph WS elevation signature is observed only at the Urt station, where its signal is mixed with that of Escos and Orthez. The resulting inference is the less accurate, with a NSE of 0.50, and closer in shape to the inferred Escos and Orthez hydrographs. This suggests correlated influences of these hydrographs at observation stations and insufficient observability of the Dax hydrograph signal given the model hydrodynamics. Observed signals at the four stations are, however, all very accurate (Fig. 17, right). For upstream stations (Orthez and Villefranque), this is due to the accurate inference of upstream hydrographs. For stations under the influence of the tidal boundary condition (BC) (Peyrehorade and Urt), this is due to the backwater influence of the BC, which compensates for the error in inferred hydrographs (namely at Dax and less so at Escos).
In conclusion, this experiment shows the capability of the VDA tool chain to infer various upstream hydrographs in 1Dlike network models. Note that multiD models are identical to the investigated 1Dlike model in terms of VDA capabilities. This paves the way towards investigations of multivariate inferability of spatially distributed hydraulic–hydrological parameters.
In the future, investigation of the influence of data sparsity, observation weights and illposed problem constraints should be carried out.
This article presents a new approach and numerical chain for the multiD hydrological–hydraulic modeling of complex river networks with variational data assimilation capabilities. It is based on the VDA algorithm and the finitevolume solvers (including a secondorder one and accurate treatment of wet–dry front propagation) from Monnier et al. (2016). The resolution of the full 2D shallowwater equations (Eq. 1) is performed with a single finitevolume solver applied to a multiD discretization of a river network domain. This lattice consists in 1Dlike reaches meshed with irregular quadrangular cells connected via 1Dlike–2D interfaces to 2D zooms consisting in higherresolution unstructured meshes – either triangular or quadrangular. This hydraulic model with inflow from a hydrological model enabling us to describe upstream/lateral catchment inflow hydrographs. In this work, the parsimonious GR4 model is integrated (Perrin et al., 2003) in its state–space version (Santos et al., 2018) for the sake of differentiability (for the VDA computations). This approach is implemented in the platform DassFlow (Monnier et al., 2019). The adjoint of the whole tool chain is generated with Tapenade (Hascoet and Pascual, 2013) and validated. Forwardinverse capabilities with the new components are assessed in several cases of increasing complexity.
The 1Dlike effective hydraulic modeling approach as well as its coupling with higherresolution 2D zooms was validated, against (1) reference 1D and 2D hydraulic models, on an array of academic cases; (2) complex real cases featuring 1D and 2D flow variabilities in river networks with confluences and floodplains. These cases are also employed to test the coupling with the hydrological model implemented in a semidistributed setup.
Considering single (multiple) type(s) of flow signature observations in those river networks through a singleobjective (multiobjective) observation cost function, the capabilities of the VDA method for inferring mono and/or multivariate control vectors of large dimensions was successfully tested.
From the obtained results, the following conclusions can be drawn:

A complete integrated multiD model coupled with an hydrological model has been implemented and validated in a parallel environment. Moreover, this complete tool chain includes VDA capabilities based on the adjoint code.

The 1Dlike modeling approach enables us to simulate fine physical flow states compared to reference 1D or 2D SW models; hydrograph propagation remains very close also.

The inference, from heterogeneous observations in the river network, of multivariate and highdimension controls among multiple inflow discharge hydrographs, bathymetry, and friction of the full shallowwater multiD hydraulic model but also hydrological model parameters is demonstrated. Very accurate inferences are obtained when the available information contained in system observability and priors is sufficient regarding the nature and quantity of unknown parameters (see discussion in Brisset et al., 2018, Larnier et al., 2020, Garambois et al., 2020, and references therein).

Information feedback from the river network to upstream hydrological models of subcatchments is shown.

Real flows on complex channel geometries can be accurately simulated with the 1Dlike model, despite its intrinsic rectangular XS, thanks to the calibration of effective geometryfriction patterns. The depthdependent friction law helps to reduce misfits across flow regimes.

Highresolution simulations of real flows can be obtained on complex river networks including floodplains and confluences with reduced simulation costs.
To our knowledge, the present numerical tool is the first one proposing largescale multiD river network modeling with VDA capabilities. This new toolchain opens the way for the resolution of largescale highdimensional hydrological–hydraulic inverse problems that can be considered given constraints from multisource datasets. The methods for direct and inverse modeling can be applied at multiple spatial scales including with a fine resolution (imposing finer calculation time steps to respect the Courant–Friedrichs–Lewy (CFL) condition in the forward hydraulic model). Building and constraining the models is, however, dependent on data availability and the informative content of observations, which may be linked to the spatial scale.
Shortterm perspectives will aim to tailor the data assimilation algorithm to perform complex data assimilation experiments at basin scale using various multisource datasets. To be actually operational, improvements pertaining to the construction of 1Dlike models from global public databases are needed to deploy the multiD approach to a large number of river networks. Coupling is ongoing with the SMASH spatially distributed hydrological platform (JayAllemand et al., 2020; Colleoni et al., 2021), on which the French flash flood warning system VigiCrues Flash (Garandeau et al., 2018) is based. Note that the addition of lateral flows (surface and subsurface runoff) would simply consist in imposing additional inflows on the boundary of the river/floodplain domain, which is numerically identical to the upstream flows imposed and inferred here (see such lateral coupling in 1D with inferences from sparse satellite altimetry data in Pujol et al., 2020). Further work will also test the integrated chain on Mediterranean basins at risk of flash flooding and on larger basins in a satellite observability context. The implementation of porosity models (Guinot et al., 2018, and references therein) represents a very interesting research direction regarding effective floodplains and 1Dlike reach parameterizations. This is especially true with depthdependent porosity (Özgen et al., 2017, and Guinot et al., 2018, and references therein) applied to complex channel geometries and with spatially distributed calibration.
A1 Twodimensional solver
Recall the rotational invariance property of the SWE (Eq. 1), which simplifies the sum of 2D problems in Eq. (3) to 1D Riemann problems. The fluxes F_{e} are computed using a Riemann solver. Each local Riemann problem depends on left and right states at the interface e.
The Harten–Lax–van Leer condition (HLLC) approximate Riemann solver is used. This gives the following expressions:
where the wave speed expressions are those proposed in Vila (1986):
It has been demonstrated that this insures L^{∞} stability, positivity and consistency with the entropy condition under a CFL condition.
For the intermediate wave speed estimate, following Toro (2001), we set
A Courant–Friedrichs–Levy (CFL) condition for the time step Δt^{n} applies; see, e.g., Vila and Villedieu (2003).
A2 Well balancing
The numerical scheme must preserve the fluid at rest property; that is, the gradient of bathymetry ∇z_{b} must not provide ${\mathbf{u}}^{n+\mathrm{1}}\ne \mathrm{0}$ if u^{n}=0. In the presence of topography gradients (in particular those perpendicular to the streamlines) the basic topography gradient ∇z_{b} discretization in the gravity source term S_{g}(U) generates spurious velocities. There is no discrete balance between the hydrostatic pressure and the gravity source term anymore: $\mathrm{\nabla}\left({\mathrm{gh}}^{{}^{\mathrm{2}}}/\mathrm{2}\right)\ne \mathrm{gh}\mathrm{\nabla}{z}_{\mathrm{b}}$.
The technique which is employed here is that presented in Audusse et al. (2004) and Audusse and Bristeau (2005). It is based on the following change in variable.
From now, at the left and right sides of edge e, the considered water depth values ${h}_{\sqcap}^{*}$ are those defined from the “reconstructed” topography z_{e} as
(see Fig. A1).
The conservative variable vector ${\mathbf{U}}_{\mathrm{e},K}^{n}$ is now considered with the new variable:
Note that this new variable ${h}_{\mathrm{e},K}^{*}$ depends on the bathymetry values z_{e,K} and z_{e}.
The resulting scheme reads
with
This provides a wellbalanced (firstorder) scheme.
A3 Predictioncorrection time scheme
The friction source term is taken into account in the complete SW system by deriving a predictioncorrection time scheme; see, e.g., Toro (2001).
We denote here by ${\overline{\mathbf{U}}}_{K}^{n+\mathrm{1}}$ the finite volume (FV) solution at time t^{n+1} of the (wellbalanced) scheme, either first or second order, of the SW system with the gravitational term S_{g} but without the friction term S_{f}. Recall that we denote $\mathbf{U}={\left(h,h\mathit{u}\right)}^{\mathrm{T}}$.
At each time step, from n to n+1, the following steps apply:

step 1: computation of ${\overline{\mathbf{U}}}^{n+\mathrm{1}}$, solution of the conservative SW system, i.e., the SW system without S_{f}, i.e., the FV solution of the following system:
$$\begin{array}{}\text{(A8)}& {\partial}_{t}\mathbf{U}+{\partial}_{x}\mathbf{F}\left(\mathbf{U}\right)+{\partial}_{y}\mathbf{G}\left(\mathbf{U}\right)={\mathbf{S}}_{\mathrm{g}}\left(\mathbf{U}\right).\end{array}$$ 
step 2: given the ”predicted value” ${\overline{\mathbf{U}}}^{n+\mathrm{1}}$, compute U^{n+1} solution of
$$\begin{array}{}\text{(A9)}& {\partial}_{t}\mathbf{U}={\mathbf{S}}_{\mathrm{f}}\left(\mathbf{U}\right).\end{array}$$
General schemes (explicit, implicit or semiimplicit ones) including the friction source term S_{f} in the discretization of the model (Eq. 1) for all K can be written as
Note that this splitting scheme is consistent at first order in Δt with the complete SWE. A splitting scheme of second order in time is possible; it is not detailed later.
In the case of the Manning–Strickler law, the friction term reads ${\mathbf{S}}_{\mathrm{f}}=g{n}^{\mathrm{2}}\left[\begin{array}{c}\mathrm{0}\\ \frac{\left\overline{\mathit{u}}\right}{{h}^{\frac{\mathrm{1}}{\mathrm{3}}}}\overline{\mathit{u}}\end{array}\right]$.
Therefore the equation to be solved (Eq. A9) reads
Since the friction source term S_{f} is zero in the mass conservation equation, we remark that ${h}^{n+\mathrm{1}}={\overline{h}}^{n+\mathrm{1}}$. As a consequence, we consider the nonzero momentum component only: ${\partial}_{t}\left(h\overline{\mathit{u}}\right)=g{n}^{\mathrm{2}}\frac{\left\overline{\mathit{u}}\right}{{h}^{\frac{\mathrm{1}}{\mathrm{3}}}}\overline{\mathit{u}}$.
A4 Firstorder expression of U^{n+1}
Let us consider the following implicit scheme:
This implies that
Let us set $c=\frac{{\left({h}^{n+\mathrm{1}}\right)}^{\frac{\mathrm{4}}{\mathrm{3}}}}{\mathrm{\Delta}{t}^{n}g{n}^{\mathrm{2}}}$; c≥0. Note that $\left{\mathit{u}}^{n+\mathrm{1}}\right{\mathit{u}}^{n+\mathrm{1}}+c{\mathit{u}}^{n+\mathrm{1}}=c{\overline{\mathit{u}}}^{n+\mathrm{1}}$. Therefore for nonvanishing velocities, there exists $\mathit{\alpha}\in \left]\mathrm{0},\mathrm{1}\right]$ such that ${\mathit{u}}^{n+\mathrm{1}}=\mathit{\alpha}{\overline{\mathit{u}}}^{n+\mathrm{1}}$. Adopting these notations, we obtain $\left{\overline{\mathit{u}}}^{n+\mathrm{1}}\right{\overline{\mathit{u}}}^{n+\mathrm{1}}{\mathit{\alpha}}^{\mathrm{2}}+c{\overline{\mathit{u}}}^{n+\mathrm{1}}\mathit{\alpha}c{\overline{\mathit{u}}}^{n+\mathrm{1}}=\mathrm{0}$. This simplifies to
Since α≥0, the root of this quadratic equation reads
Let us define the function $\mathit{\epsilon}=\frac{\mathrm{1}}{c}\left{\overline{\mathit{u}}}^{n+\mathrm{1}}\right=\mathrm{\Delta}{t}^{n}g{n}^{\mathrm{2}}\frac{\left{\overline{\mathit{u}}}^{n+\mathrm{1}}\right}{{\left({h}^{n+\mathrm{1}}\right)}^{\mathrm{4}/\mathrm{3}}}$. Observe that ε=O(Δt^{n}); also $\mathit{\epsilon}=O\left(\frac{\left{\overline{\mathit{u}}}^{n+\mathrm{1}}\right}{{\left({h}^{n+\mathrm{1}}\right)}^{\mathrm{4}/\mathrm{3}}}\right)$. By adopting this notation, Eq. (A15) reads $\mathit{\alpha}=\frac{\sqrt{\mathrm{1}+\mathrm{4}\mathit{\u03f5}}\mathrm{1}}{\mathrm{2}\mathit{\u03f5}}$. After some rearrangements, we obtain $\mathit{\alpha}=\frac{\mathrm{2}}{\mathrm{1}+\sqrt{\mathrm{1}+\mathrm{4}\mathit{\epsilon}}}$. At first order in ϵ, we get $\mathit{\alpha}\sim \left(\frac{\mathrm{1}}{\mathrm{1}+\mathit{\u03f5}/\mathrm{4}}\right)\sim \mathrm{1}\mathit{\u03f5}/\mathrm{4}$. Finally, we obtain
with $\left{\overline{\mathit{u}}}^{n+\mathrm{1}}\right$ the solution of Eq. (A8).
A5 Secondorder MUSCL scheme
In order to obtain a globally secondorder scheme, a higherorder time stepping scheme is needed. Let us briefly describe the ingredients of this secondorder wellbalanced positive scheme that is strictly the same as the one proposed in Couderc et al. (2013) and Monnier et al. (2016). Actual secondorder accuracy, considering source terms S_{g} and S_{f}, is achieved through the combination of a Monotonic Upwind Scheme for Conservation Laws (MUSCL) spatial reconstruction and an IMEX RK time scheme (see Monnier et al., 2016, and references therein) as well as a spatial discretization of S_{g} and the semiimplicit nature of the friction source term S_{f} given in the subsection above.
A monoslope secondorder MUSCL scheme is adopted; see, e.g., Chévrier and Galley (1993) and Buffard and Clain (2010). It leads to new expressions of ${\mathbf{U}}_{K}^{n}$ and ${\mathbf{U}}_{{K}_{\mathrm{e}}}^{n}$. With this linear reconstruction, one can expect a scheme with a secondorder accuracy in space (for regular solutions only). In order to prevent large numerical dispersive instabilities, the computed vectorial slopes are limited by applying a maximum principle. Furthermore, to handle the presence of wet–dry fronts that can break the finitevolume mass conservation property, a Barth limiter (Barth, 2003) is employed.
A river network Ω^{1D} is described by connected closed line segments. For a flow XS orthogonal to the main (longitudinal) flow direction of curvilinear abscissa x∈Ω^{1D} (distance from downstream areas) at time $t\in \left[\mathrm{0},T\right]$, let A(x,t) be the flow crosssectional area [m^{2}] and Q(x,t) the discharge [m^{3} s^{−1}] such that Q=UA with U(x,t) defined as the longitudinal XS averaged velocity [m s^{−1}]. The 1D SaintVenant equations in the (A,Q) variables at a flow XS are written as follows:
where Z(x,t) is the WS elevation [m] and $Z=\left(b+h\right)$ (with b(x) [m] the riverbed level and h(x,t) [m] the water depth), ${R}_{\mathrm{h}}\left(x,t\right)=A/{P}_{\mathrm{h}}$ [m] the hydraulic radius, P_{h}(x,t) [m] the wetted perimeter, and g [m s^{−2}] is the gravity magnitude. Let us recall the Froude number definition $\mathrm{Fr}=U/c$ comparing the average flow velocity U to pressure wave celerity $c=\sqrt{\frac{gA}{W}}$, where W [m] is the flow top width.
The friction term S_{f} is classically parameterized with the empirical Manning–Strickler law established for uniform flows $\frac{\leftQ\rightQ}{{K}^{\mathrm{2}}{A}^{\mathrm{2}}{R}_{h}^{\mathrm{4}/\mathrm{3}}}$, where K $\left[{\mathrm{m}}^{\mathrm{1}/\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}\right]$ is the Strickler coefficient.
The SaintVenant equations are solved on each segment of the river network and the continuity of the flow between segments is ensured by applying an equality constraint on water levels at the confluence between two segments.
Boundary conditions are classically imposed (subcritical flows here) at boundary nodes with inflow discharges Q(t) at upstream nodes and WS elevation Z(t) at the downstream node; there are lateral hydrographs q_{l}(t) at in/outflow nodes. The initial condition is set as the steadystate backwater curve profile ${Z}_{\mathrm{0}}\left(x\right)=Z\left({Q}_{\mathrm{in}}\left({t}_{\mathrm{0}}\right),{q}_{\mathrm{l},\mathrm{1}\mathrm{\dots}L}\left({t}_{\mathrm{0}}\right)\right)$ for a hot start. This 1D SaintVenant model is discretized using the classical implicit Preissmann scheme (see, e.g., Cunge et al., 1980) on a regular grid of spacing Δx using a double sweep method enabling us to deal with flow regimes changes with a 1 h time step Δt. This is implemented in the computational software DassFlow (Brisset et al., 2018; Larnier, 2010).
The numerical scheme is a semiimplicit finitedifference scheme (generalized Preissmann scheme) with a double sweep Local Partial Inertial method to minimize the inertial terms (see documentation in Brisset et al., 2018, and Larnier, 2010).
The state–space version of the lumped conceptual hydrological model GR4 presented in Santos et al. (2018) consists in the set of ODEs below:
They involve the following fluxes:
The following parameters are set following Perrin et al. (2003) and Santos et al. (2018): α=2, β=5, γ=5, ω=3.5, $\mathit{\nu}=\mathrm{4}/\mathrm{9}$, Φ=0.9, n_{res}=11.
Calibrated parameters for the Adour case were obtained using the airGR global calibration algorithm (Coron et al., 2017) from the freely available package (https://webgr.inrae.fr/logiciels/airgr/, last access: 27 July 2022).
^{∗} Note that the x_{4} calibrated parameter corresponds to the non“state–space” GR4H version (not presented, see Perrin et al., 2003), for which the calibration tool is provided. x_{4} values in the present run with the state–space were set at 0.15.
The current DassFlow code is available on demand at https://www.math.univtoulouse.fr/DassFlow/download.html (last access: 27 July 2022). The exact version of the model used in this article is archived on Zenodo (https://doi.org/10.5281/zenodo.6342723, Pujol et al., 2022), as are input data needed to run the main direct/inverse simulations.
This work is part of the PhD work of LP under the supervision of PAG and JM. The research plan was developed by PAG, LP and JM. The computational software DassFlow2D was adapted from a previous version by LP under the supervision of PAG and JM. Numerical investigations were carried out by LP with analyses by all. JM designed or supervised the 2D numerical scheme and VDA algorithm. PAG and JM designed and cosupervised the 1D–2D algorithm. All authors participated in the writing. PAG was responsible for the project management for funding.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This study was rendered possible by the adaptation of covariance matrices in the inverse algorithm from the DassFlow1D algorithm (Larnier et al., 2020) by Lilian Villenave (INRAE AixenProvence) and the implementation of final tricks in Perl scripts for adjoint postprocessing in DassFlow by Kévin Larnier (CS Group Toulouse). The authors greatly acknowledge Pascal FinaudGuyot (HydroSciences, University of Montpellier) for fruitful discussions related to hydraulic solvers. They warmly thank Laurent Dieval from SPCGAD France, the operational flood forecasting service for the Gaves–Adour–Dordogne river networks, for providing data for the Adour case and also for interesting discussions. They also thank Robert Mosé (ICUBE, University of Strasbourg), the PhD director of Léo Pujol.
Léo Pujol has been cofunded by CNES and ICUBE.
This paper was edited by Charles Onyutha and reviewed by two anonymous referees.
Allen, M., AntwiAgyei, P., AragonDurand, F., Babiker, M., Bertoldi, P., Bind, M., Brown, S., Buckeridge, M., Camilloni, I., Cartwright, A., MassonDelmotte, V., Zhai, P., Pörtner, H.O., Roberts, D., Skea, J., Shukla, P. R., Pirani, A., MoufoumaOkia, W., Péan, C., Pidcock, R., Connors, S., Matthews, J. B. R., Chen, Y., Zhou, X., Gomis, M. I., Lonnoy, E., Maycock, T., Tignor, M., and Waterterfield, T. (Eds.): Technical Summary: Global warming of 1.5 ^{∘}C, An IPCC Special Report on the impacts of global warming of 1.5 ^{∘}C above preindustrial levels and related global greenhouse gas emission pathways, in the context of strengthening the global response to the threat of climate change, sustainable development, and efforts to eradicate poverty, http://pure.iiasa.ac.at/15716 (last access: 27 June 2022), 2019. a
Amara, M., CapatinaPapaghiuc, D., and Trujillo, D.: Hydrodynamical modelling and multidimensional approximation of estuarian river flows, Comput. Vis. Sci., 6, 39–46, https://doi.org/10.1007/s007910030106z, 2004. a
Asch, M., Bocquet, M., and Nodet, M.: Data assimilation: methods, algorithms, and applications, Fundamentals of Algorithms, SIAM, https://hal.inria.fr/hal01402885 (last access: 27 June 2022), 2016. a
Audusse, E. and Bristeau, M.O.: A wellbalanced positivity preserving “secondorder” scheme for shallow water flows on unstructured meshes, J. Comput. Phys., 206, 311–333, https://doi.org/10.1016/j.jcp.2004.12.016, 2005. a
Audusse, E., Bouchut, F., Bristeau, M.O., Klein, R., and Perthame, B.: A Fast and Stable WellBalanced Scheme with Hydrostatic Reconstruction for Shallow Water Flows, SIAM J. Sci. Comput., 25, 2050–2065, https://doi.org/10.1137/S1064827503431090, 2004. a
Barth, T.: Numerical Methods for conservative Laws on Structured and Unstructured Meshes, Tech. rep., VKI Lecture Series, 2003. a
Barthélémy, S., Ricci, S., Morel, T., Goutal, N., Le Pape, E., and Zaoui, F.: On operational flood forecasting system involving 1D$/$2D coupled hydraulic model and data assimilation, J. Hydrol., 562, 623–634, https://doi.org/10.1016/j.jhydrol.2018.05.007, 2018. a
Bates, P., Trigg, M., Neal, J., and Dabrowa, A.: LISFLOODFP, User manual, School of Geographical Sciences, University of Bristol, Bristol, UK, https://vdocument.in/lisfloodfpusermanualuniversityofbristol (last access: 1 August 2022), 2013. a
Bates, P. D., Horritt, M. S., and Fewtrell, T. J.: A simple inertial formulation of the shallow water equations for efficient twodimensional flood inundation modelling, J. Hydrol., 387, 33–45, https://doi.org/10.1016/j.jhydrol.2010.03.027, 2010. a
Bertalanffy, L. v.: General Systems Theory, 1968. a
Beven, K.: Prophecy, reality and uncertainty in distributed hydrological modelling, Adv. Water Resour., 16, 41–51, 1993. a
Beven, K. J.: Uniqueness of place and process representations in hydrological modelling, Hydrol. Earth Syst. Sci., 4, 203–213, https://doi.org/10.5194/hess42032000, 2000. a
Biancamaria, S., Lettenmaier, D., and Pavelsky, T.: The SWOT Mission and Its Capabilities for Land Hydrology, Surv. Geophys., 37, 307–337, https://doi.org/10.1007/s107120159346y, 2016. a
Biancamaria, S., Frappart, F., Leleu, A.S., Marieu, V., Blumstein, D., Desjonquères, J.D., Boy, F., Sottolichio, A., and ValleLevinson, A.: Satellite radar altimetry water elevations performance over a 200 m wide river: Evaluation over the Garonne River, Adv. Space Res., 59, 128–146, https://doi.org/10.1016/j.asr.2016.10.008, 2017. a
Bouttier, F. and Courtier, P.: Data assimilation concepts and methods March 1999, Meteorological training course lecture series, ECMWF, 718, 59, http://msi.ttu.ee/~elken/Assim_concepts.pdf (last access: 27 June 2022), 2002. a, b
Brêda, J., Paiva, R., Bravo, J., Passaia, O., and Moreira, D.: Assimilation of satellite altimetry data for effective river bathymetry, Water Resour. Res., 55, 7441–7463, https://doi.org/10.1029/2018WR024010, 2019. a
Brisset, P., Monnier, J., Garambois, P.A., and Roux, H.: On the assimilation of altimetric data in 1D SaintVenant river flow models, Adv. Water Res., 119, 41–59, https://doi.org/10.1016/j.advwatres.2018.06.004, 2018. a, b, c, d, e, f, g, h, i, j, k
Brunner, G. W.: HECRAS River Analysis System, Hydraulic Reference Manual, Version 1.0., Tech. rep., Hydrologic Engineering Center Davis, CA, https://apps.dtic.mil/sti/pdfs/ADA311952.pdf (last access: 27 June 2022) 1995. a
Buffard, T. and Clain, S.: Monoslope and Multislope MUSCL Methods for unstructured meshes, J. Comput. Phys., 229, 3745–3776, https://doi.org/10.1016/j.jcp.2010.01.026, 2010. a
Chévrier, P. and Galley, H.: A Van Leer finite volume scheme for the Euler equations on unstructured meshes, ESAIMMath. Model. Num., 27, 183–201, 1993. a
Chow, V.: Openchannel Hydraulics, McGrawHill civil engineering series, McGrawHill, NewYork, USA, ISBN 9780070859067, 1959. a
Colleoni, F., Garambois, P.A., Javelle, P., JayAllemand, M., I., G., Organde, D., and Arnaud, P.: SMASH v1.0 platform for spatially distributed hydrological modeling and data assimilation: hypothesis testing and signatures analysis, J. Hydrol., submitted, 2021. a, b
Collischonn, W., Allasia, D., Da Silva, B. C., and Tucci, C. E.: The MGBIPH model for largescale rainfallrunoff modelling, Hydrolog. Sci. J., 52, 878–895, https://doi.org/10.1623/hysj.52.5.878, 2007. a
Coron, L., Thirel, G., Delaigue, O., Perrin, C., and Andréassian, V.: The suite of lumped GR hydrological models in an R package, Environ. Modell. Softw., 94, 166–171, https://doi.org/10.1016/j.envsoft.2017.05.002, 2017. a, b
Couderc, F., Madec, R., Monnier, J., and Vila, J.P.: DassFowShallow, Variational Data Assimilation for ShallowWater Models: Numerical Schemes, User and Developer Guides, University of Toulouse, CNRS, IMT, INSA, ANR, Research report, https://hal.archivesouvertes.fr/hal01120285 (last access: 27 June 2022), 2013. a, b, c, d, e
Cunge, J. A., Holly, F., M., and Verwey, A.: Practical Aspects of Computational River Hydraulics, Pitam Publishing, ISBN 9780273084426, 1980. a
Davy, P., Croissant, T., and Lague, D.: A precipiton method to calculate river hydrodynamics, with applications to flood prediction, landscape evolution models, and braiding instabilities, J. Geophys. Res.Earth, 122, 1491–1512, https://doi.org/10.1002/2016JF004156, 2017. a, b
Delestre, O., Darboux, F., James, F., Lucas, C., Laguerre, C., and Cordier, S.: FullSWOF: Full ShallowWater equations for Overland Flow, J. Open Source Softw., 2, 448, https://doi.org/10.21105/joss.00448, 2017. a
Dingman, S. L.: Fluvial hydraulics, Oxford University Press, ISBN 9780195172867, 2009. a, b
FinaudGuyot, P., Garambois, P.A., Chen, S., Dellinger, G., Ghenaim, A., and Terfous, A.: 1D$/$2D porosity model for urban flood modeling: case of a dense street networks, E3S Web Conf., 40, 06010, https://doi.org/10.1051/e3sconf/20184006010, 2018. a
Fleischmann, A. S., Paiva, R. C. D., Collischonn, W., Siqueira, V. A., Paris, A., Moreira, D. M., Papa, F., Bitar, A. A., Parrens, M., Aires, F., and Garambois, P. A.: Tradeoffs between 1D and 2D regional river hydrodynamic models, Water Resour. Res., 56, e2019WR026812, https://doi.org/10.1029/2019WR026812, 2020. a, b, c
Galland, J.C., Goutal, N., and Hervouet, J.M.: TELEMAC: A new numerical model for solving shallow water equations, Adv. Water Resour., 14, 138–148, 1991. a
Garambois, P., Calmant, S., Roux, H., Paris, A., Monnier, J., FinaudGuyot, P., Montazem, A., and da Silva, J.: Hydraulic visibility: Using satellite altimetry to parameterize a hydraulic model of an ungauged reach of a braided river, Hydrol. Process., 31, 756–767, https://doi.org/10.1002/hyp.11033, 2017. a, b
Garambois, P.A. and Monnier, J.: Inference of effective river properties from remotely sensed observations of water surface, Adv. Water Resour., 79, 103–120, https://doi.org/10.1016/j.advwatres.2015.02.007, 2015. a, b
Garambois, P.A., Larnier, K., Monnier, J., FinaudGuyot, P., Verley, J., Montazem, A., and Calmant, S.: Variational estimation of effective channel and ungauged anabranching river discharge from multisatellite water heights of different spatial sparsity, J. Hydrol., 581, 124409, https://doi.org/10.1016/j.jhydrol.2019.124409, 2020. a, b, c, d
Garandeau, L., Belleudy, A., Javelle, P., Organde, D., Janet, B., Demargne, J., De SaintAubin, C., and Fouchier, C.: Vigicrues Flash, un service automatique d'avertissement pour les crues rapides, De la prévision des crues à la gestion de crise, Société Hydrotechnique de France, Avignon, France, 11, https://hal.inrae.fr/hal02608801 (last access: 27 June 2022), 2018. a
Gejadze, I. Y. and Monnier, J.: On a 2D zoom for the 1D shallow water model: Coupling and data assimilation, Comput. Method. Appl. M., 196, 4628–4643, https://doi.org/10.1016/j.cma.2007.05.026, 2007. a
Gervasio, P., Lions, J.L., and Quarteroni, A.: Heterogeneous coupling by virtual control methods, Numer. Math., 90, 241–264, https://doi.org/10.1007/s002110100303, 2001. a
Goutal, N. and Maurel, F.: A finite volume solver for 1D shallowwater equations applied to an actual river, Int. J. Numer. Meth. Fl., 38, 1–19, https://doi.org/10.1002/fld.201, 2002. a
Grimaldi, S., Li, Y., Walker, J., and Pauwels, V.: Effective representation of river geometry in hydraulic flood forecast models, Water Resour. Res., 54, 1031–1057, https://doi.org/10.1002/2017WR021765, 2018. a, b
Gudmundsson, G. H.: Transmission of basal variability to a glacier surface, J. Geophys. Res.Sol. Ea., 108, B5, https://doi.org/10.1029/2002JB002107, 2003. a
Guinot, V.: Wave propagation in fluids: models and numerical techniques, 2nd edn., vol. 49, edited by: ISTE Ltd., ISBN 9789812707789, 2010. a
Guinot, V., Delenne, C., Rousseau, A., and Boutron, O.: Flux closures and source term models for shallow water models with depthdependent integral porosity, Adv. Water Resour., 122, 1–26, https://doi.org/10.1016/j.advwatres.2018.09.014, 2018. a, b, c
Hascoet, L. and Pascual, V.: The Tapenade automatic differentiation tool: principles, model, and specification, ACM T. Math. Software, 39, 1–43, https://doi.org/10.1145/2450153.2450158, 2013. a, b, c, d
Hocini, N., Payrastre, O., Bourgin, F., Gaume, E., Davy, P., Lague, D., Poinsignon, L., and Pons, F.: Performance of automated methods for flash flood inundation mapping: a comparison of a digital terrain model (DTM) filling and two hydrodynamic methods, Hydrol. Earth Syst. Sci., 25, 2979–2995, https://doi.org/10.5194/hess2529792021, 2021. a, b, c
Hostache, R., Lai, X., Monnier, J., and Puech, C.: Assimilation of spatially distributed water levels into a shallowwater flood model. Part II: Use of a remote sensing image of Mosel River, J. Hydrol., 390, 257–268, https://doi.org/10.1016/j.jhydrol.2010.07.003, 2010. a
Hunter, N. M., Bates, P. D., Neelz, S., Pender, G., Villanueva, I., Wright, N. G., Liang, D., Falconer, R. A., Lin, B., Waller, S., Crossley, A. J., and Mason, D. C.: Benchmarking 2D hydraulic models for urban flooding, in: Proceedings of the Institution of Civil EngineersWater Management, Thomas Telford Ltd, 161, 13–30, https://doi.org/10.1680/wama.2008.161.1.13, 2008. a
Iturbide, M., Gutiérrez, J. M., Alves, L. M., Bedia, J., CerezoMota, R., Cimadevilla, E., Cofiño, A. S., Di Luca, A., Faria, S. H., Gorodetskaya, I. V., Hauser, M., Herrera, S., Hennessy, K., Hewitt, H. T., Jones, R. G., Krakovska, S., Manzanas, R., MartínezCastro, D., Narisma, G. T., Nurhati, I. S., Pinto, I., Seneviratne, S. I., van den Hurk, B., and Vera, C. S.: An update of IPCC climate reference regions for subcontinental analysis of climate model data: definition and aggregated datasets, Earth Syst. Sci. Data, 12, 2959–2970, https://doi.org/10.5194/essd1229592020, 2020. a
JayAllemand, M., Javelle, P., Gejadze, I., Arnaud, P., Malaterre, P.O., Fine, J.A., and Organde, D.: On the potential of variational calibration for a fully distributed hydrological model: application on a Mediterranean catchment, Hydrol. Earth Syst. Sci., 24, 5519–5538, https://doi.org/10.5194/hess2455192020, 2020. a, b
Kirstetter, G., Delestre, O., Lagrée, P.Y., Popinet, S., and Josserand, C.: Bflood 1.0: an opensource SaintVenant model for flashflood simulation using adaptive refinement, Geosci. Model Dev., 14, 7117–7132, https://doi.org/10.5194/gmd1471172021, 2021. a
Lai, X. and Monnier, J.: Assimilation of spatially distributed water levels into a shallowwater flood model. Part I: Mathematical method and test case, J. Hydrol., 377, 1–11, https://doi.org/10.1016/j.jhydrol.2009.07.058, 2009. a
Larnier, K.: Modélisation thermohydraulique d'un troncon de Garonne en lien avec l'habitat piscicole: Approches statistique et déterministe, PhD thesis, Institut National Polytechnique de Toulouse, http://ethesis.inptoulouse.fr/archive/00001263/ (last access: 27 June 2022), 2010. a, b
Larnier, K., Monnier, J., Garambois, P.A., and Verley, J.: River discharge and bathymetry estimation from SWOT altimetry measurements, Inverse Probl. Sci. Eng., 29, 759–789, https://doi.org/10.1080/17415977.2020.1803858, 2020. a, b, c, d, e, f, g, h, i, j, k, l, m
Le Lay, M.: Modélisation hydrologique dans un contexte de variabilité hydroclimatique. Une approche comparative pour l'étude du cycle hydrologique à mésoéchelle au Bénin, PhD thesis, Institut National Polytechnique de Grenoble (INPG), https://tel.archivesouvertes.fr/tel00116912 (last access: 27 June 2022), 2006. a
Leopold, L. B. and Maddock, T.: The hydraulic geometry of stream channels and some physiographic implications, US Government Printing Office, 252, https://doi.org/10.3133/pp252, 1953. a
Lorenc, A. C., Ballard, S. P., Bell, R. S., Ingleby, N. B., Andrews, P. L. F., Barker, D. M., Bray, J. R., Clayton, A. M., Dalby, T., Li, D., Payne, T. J., and Saunders, F. W.: The Met. Office global threedimensional variational data assimilation scheme, Q. J. Roy. Meteor. Soc., 126, 2991–3012, 2000. a
Malou, T. and Monnier, J.: Doublescale diffusive wave model dedicated to spatial river observation and associated covariance kernel for variational data assimilation , EGU General Assembly 2021, online, 19–30 Apr 2021, EGU2110355, https://doi.org/10.5194/egusphereegu2110355, 2021. a
Malou, T. and Monnier, J.: Covariance kernels investigation from diffusive wave equations for data assimilation in hydrology, Inverse Probl., https://doi.org/10.1088/13616420/ac509d, accepted, 2022. a
Malou, T., Garambois, P.A., Paris, A., Monnier, J., and Larnier, K.: Generation and analysis of stagefalldischarge laws from coupled hydrologicalhydraulic river network model integrating sparse multisatellite data, J. Hydrol., 603, 126993, https://doi.org/10.1016/j.jhydrol.2021.126993, 2021. a, b, c
Marin, J. and Monnier, J.: Superposition of local zoom models and simultaneous calibration for 1D–2D shallow water flows, Math. Comput. Simulat., 80, 547–560, https://doi.org/10.1016/j.matcom.2009.09.001, 2009. a
Martin, N. and Monnier, J.: Inverse rheometry and basal properties inference for pseudoplastic geophysical flows, Eur. J. Mech. BFluid., 50, 110–126, https://doi.org/10.1016/j.euromechflu.2014.11.011, 2015. a
Miglio, E., Perotto, S., and Saleri, F.: Model coupling techniques for freesurface flow problems: Part I, Nonlinear Anal.Theor., 63, e1885–e1896, https://doi.org/10.1016/j.na.2005.03.083, 2005a. a
Miglio, E., Perotto, S., and Saleri, F.: Model coupling techniques for freesurface flow problems: Part II, Nonlinear Anal.Theor., 63, e1897–e1908, https://doi.org/10.1016/j.na.2005.03.085, 2005b. a
Milly, P.: Climate, soil water storage, and the average annual water balance, Water Resour. Res., 30, 2143–2156, https://doi.org/10.1029/94WR00586, 1994. a
Monnier, J.: Variational data assimilation: from optimal control to large scale data assimilation, Open Online Course, INSA Toulouse, https://www.math.univtoulouse.fr/~jmonnie/Enseignement/CourseVDA.pdf (last access: 27 June 2022), 2014. a
Monnier, J.: Variational Data Assimilation and Model Learning, https://hal.archivesouvertes.fr/hal03040047 (last access: 1 August 2022), 2021. a
Monnier, J., Couderc, F., Dartus, D., Larnier, K., Madec, R., and Vila, J.P.: Inverse algorithms for 2D shallow water equations in presence of wet dry fronts: Application to flood plain dynamics, Adv. Water Resour., 97, 11–24, https://doi.org/10.1016/j.advwatres.2016.07.005, 2016. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q
Monnier, J., Couderc, F., and Vila, J.P.: Data Assimilation for Free Surface Flows, Mathematics Institute of Toulouse – INSA corp. CNESCNRS, Tech. rep., http://www.math.univtoulouse.fr/DassFlow (last access: 27 June 2022), 2019. a, b, c
Montazem, A., Garambois, P.A., Calmant, S., FinaudGuyot, P., Monnier, J., Moreira, D., Minear, J., and Biancamaria, S.: WaveletBased River Segmentation Using Hydraulic ControlPreserving Water Surface Elevation Profile Properties, Geophys. Res. Lett., 46, 6534–6543, https://doi.org/10.1029/2019GL082986, 2019. a, b, c
Nguyen, P., Thorstensen, A., Sorooshian, S., Hsu, K., AghaKouchak, A., Sanders, B., Koren, V., Cui, Z., and Smith, M.: A high resolution coupled hydrologichydraulic model (HiResFloodUCI) for flash flood modeling, J. Hydrol., 541, 401–420, https://doi.org/10.1016/j.jhydrol.2015.10.047, 2016. a, b, c
Oudin, L., Hervieu, F., Michel, C., Perrin, C., Andréassian, V., Anctil, F., and Loumagne, C.: Which potential evapotranspiration input for a lumped rainfallrunoff model?: Part 2 Towards a simple and efficient potential evapotranspiration model for rainfallrunoff modelling, J. Hydrol., 303, 290–306, https://doi.org/10.1016/j.jhydrol.2004.08.026, 2005. a
Özgen, I., Zhao, J.h., Liang, D.f., and Hinkelmann, R.: Wave propagation speeds and source term influences in single and integral porosity shallow water equations, Water Sci. Eng., 10, 275–286, https://doi.org/10.1016/j.wse.2017.12.003, 2017. a
Perrin, C., Michel, C., and Andréassian, V.: Improvement of a parsimonious model for streamflow simulation, J. Hydrol., 279, 275–289, https://doi.org/10.1016/S00221694(03)002257, 2003. a, b, c, d, e, f, g
Pontes, P., Fan, F., Fleischmann, A., Paiva, R., Buarque, D., Siqueira, V., Jardim, P., Sorribas, M., and Collischonn, W.: MGBIPH model for hydrological and hydraulic simulation of large floodplain river systems coupled with open source GIS, Environ. Modell. Softw., 94, 1–20, https://doi.org/10.1016/j.envsoft.2017.03.029, 2017. a
Pujol, L., Garambois, P.A., FinaudGuyot, P., Monnier, J., Larnier, K., Mosé, R., Biancamaria, S., Yesou, H., Moreira, D., Paris, A., and Calmant, S.: Estimation of multiple inflows and effective channel by assimilation of multisatellite hydraulic signatures: The ungauged anabranching Negro river, J. Hydrol., 591, 125331, https://doi.org/10.1016/j.jhydrol.2020.125331, 2020. a, b, c, d, e, f, g, h, i
Pujol, L., Garambois, P.A., and Monnier, J.: DassFlow2DV3 code and cases, Zenodo [code], https://doi.org/10.5281/zenodo.6342723, 2022. a, b
Rodríguez, E., Durand, M., and Frasson, R. P. d. M.: Observing rivers with varying spatial scales, Water Resour. Res., 56, 9, https://doi.org/10.1029/2019WR026476, 2020. a
Sanders, B. F., Schubert, J. E., and Detwiler, R. L.: ParBreZo: A parallel, unstructured grid, Godunovtype, shallowwater code for highresolution flood inundation modeling at the regional scale, Adv. Water Resour., 33, 1456–1467, https://doi.org/10.1016/j.advwatres.2010.07.007, 2010. a, b
Santos, L., Thirel, G., and Perrin, C.: Continuous statespace representation of a buckettype rainfallrunoff model: a case study with the GR4 model using statespace GR4 (version 1.0), Geosci. Model Dev., 11, 1591–1605, https://doi.org/10.5194/gmd1115912018, 2018. a, b, c, d, e, f, g, h
Schuite, J., Flipo, N., Massei, N., Rivière, A., and Baratelli, F.: Improving the Spectral Analysis of Hydrological Signals to Efficiently Constrain Watershed Properties, Water Resour. Res., 55, 4043–4065, https://doi.org/10.1029/2018WR024579, 2019. a
Schumann, G. J.P. and Domeneghetti, A.: Exploiting the proliferation of current and future satellite observations of rivers, Hydrol. Process., 30, 2891–2896, https://doi.org/10.1002/hyp.10825, 2016. a
Steinstraesser, J. G. C., Delenne, C., FinaudGuyot, P., Guinot, V., Casapia, J. K., and Rousseau, A.: SW2DLEMON: a new software for upscaled shallow water modeling, in: Simhydro 2021 – 6th International Conference Models for complex and global water issuesPractices and expectations, Sophia Antipolis, 16–18 June 2021, https://hal.inria.fr/hal03224050/ (last access: 27 June 2022), 2021. a
Toro, E.: Shockcapturing methods for freesurface shallow flows, Wiley Blackwell, ISBN 9780471987666, 2001. a, b
Uhe, P., Mitchell, D., Bates, P. D., Addor, N., Neal, J., and Beck, H. E.: Model cascade from meteorological drivers to river flood hazard: floodcascade v1.0, Geosci. Model Dev., 14, 4865–4890, https://doi.org/10.5194/gmd1448652021, 2021. a, b
Vila, J.P.: Simplified Godunov schemes for 2×2 systems of conservation laws, SIAM J. Numer. Anal., 23, 1173–1192, https://doi.org/10.1137/0723079, 1986. a
Vila, J.P. and Villedieu, P.: Convergence of an explicit finite volume scheme for first order symmetric systems, Numer. Math., 94, 573–602, https://doi.org/10.1007/s002110020396y, 2003. a
Zhu, C., Byrd, R. H., Lu, P., and Nocedal, J.: Algorithm 778: LBFGSB: Fortran subroutines for largescale boundconstrained optimization, ACM T. Math. Software, 23, 550–560, https://doi.org/10.1145/279232.279236, 1997. a, b, c
https://www.aquaveo.com/software (last access: 27 June 2022); a software including comprehensive meshing tools.
http://www.hydro.eaufrance.fr (last access: 27 June 2022); Office français pour la biodiversité.
 Abstract
 Introduction
 The computational hydrological–hydraulic chain
 Results and discussion
 Conclusions and perspectives
 Appendix A: 2D $\left(h,u,v\right)$ shallowwater scheme
 Appendix B: 1D (A,Q) SaintVenant
 Appendix C: GR4 hydrological model operators
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 The computational hydrological–hydraulic chain
 Results and discussion
 Conclusions and perspectives
 Appendix A: 2D $\left(h,u,v\right)$ shallowwater scheme
 Appendix B: 1D (A,Q) SaintVenant
 Appendix C: GR4 hydrological model operators
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References