Multi-dimensional hydrological-hydraulic model with variational data assimilation for river networks and floodplains

. This contribution presents a novel multi-dimensional (multi-D) hydraulic-hydrological numerical model with variational data assimilation capabilities. It allows multi-scale modeling over large domains, combining in situ observations with high-resolution hydro-meteorology and satellite data. The multi-D hydraulic model relies on the 2D shallow water equations solved with a 1D2D adapted single finite volume solver. 1Dlike 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 5 as floodplains and confluences, via 1Dlike-2D interfaces. An existing parsimonious hydrological model, GR4H, is implemented and coupled to the hydraulic model. The forward-inverse multi-D computational model is successfully validated on virtual and real cases of increasing complexity, including using the second order scheme version. Assimilating multiple observations of flow signatures leads to accurate inferences of multi-variate and spatially distributed parameters among bathymetry-friction, upstream/lateral hydrographs, and hydrological model parameters. This notably demonstrates the possibility for information 10 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 multi-D model of the complex Adour basin network, inflowed by the semi-distributed hydrological model, is built. High resolution 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 15 chain is obtained by automatic code differentiation. on multi-variate inferability of spatially distributed hydraulic-hydrological parameters.


Introduction
The accurate estimation of storage and fluxes in surface hydrology is an essential scientific question linked to major socioeconomic issues in floods and droughts 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, ad-20 vanced numerical modeling tools are crucially needed to both perform meaningful and detailed representations of basin-scale hydrological processes and provide sensible local forecasts. The quantities of interest range from discharge hydrographs on inertial hydraulic modeling) or at finer scale, e.g. at catchment scale for flash floods in Nguyen et al. (2016); Hocini et al. (2020).
In those studies, conceptual hydrological models of upstream-lateral sub-catchments are used to inflow hydraulic models of river network and floodplains in a weak coupling approach, mostly performed via external coupling of numerical models. In 60 Grimaldi et al. (2018); Fleischmann et al. (2020); Uhe et al. (2020), a simple 2D storage cell inundation model obtained from 1D non-inertial model (Bates et al. (2010) following Hunter et al. (2008), implemented in LISFLOOD-FP model), enables raster based inundation modeling over very large domains at relatively low computational cost (see also Fleischmann et al. (2020) for coupling of this non inertial model with the large scale 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 This work focuses on the following scientific problems: i) the effective representation of river networks flows through multi-D basin-scale hydraulic-hydrological numerical modeling, ii) the inference from multi-source data of spatio-temporally distributed river channel parameters and inflows through a VDA approach, and iii) the upstream informational feedback of river network observations to integrated hydrological models. 95 The present study details upgrades to the DassFlow variational data assimilation framework (DassFlow2D-V2, see Monnier et al. (2016)) in the form of (i) a new multi-D hydraulic computational code relying on an internal mesh coupling and a single finite volume 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 DassFlow2D-V3.
One strength of the DassFlow platform is its variational data assimilation (VDA) algorithm. It can solve high-dimensional 100 optimization problems with descent algorithms using gradients computed with the adjoint model. The adjoint model being 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 multi-variate data assimilation problems, i.e. with composite control vectors (bathymetry, friction, bound-105 ary conditions, hydrological parameters), given multi-source datasets, heterogeneous in nature and spatio-temporal resolutions (see e.g. Brisset et al. (2018); Larnier et al. (2020); Pujol et al. (2020);  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 modelling, by proposing this model as a single hydraulic-hydrological tool/code. The proposed multi-D hydraulic code consists 110 in a single finite volume solver applied to a 2D river network. The network is discretized into "1Dlike" reaches connected to high resolution 2D meshes in a single formulation of the Shallow Water Equations (SWE). 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 modelling of river network is achieved by coupling this hydraulic domain with a well-115 established conceptual hydrological model (GR4H state-space Santos et al. (2018)) in a semi-distributed setup. The resulting integrated tool chain has been differentiated and enables, via the VDA algorithm, 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 1 , and current updated version of the code upon simple request 2 .
The remainder of this article is organized as follows. In Section 2, the modeling hypothesis, the computational resolution and 120 inverse methods are detailed. In Section 3, the multi-D coupling scheme is validated on a series of academic cases and several academic and real-like inference setups are investigated. The study is concluded in Section 4, which also outlines potential applications and improvement perspectives that the proposed method and findings bring.  This section presents the integrated and multi-dimensional hydrological-hydraulic model and the data assimilation approach.

125
The model is designed for simulating spatio-temporal flow variabilities over an entire river network, from upstream hydrological responses to complex flow zones (confluences, multichannels portions, floodplains, ...).
The modeling approach which is detailed below, is based on the following ingredients: -An integrated multi-D hydraulic model: the 2D SWE with finite volume solvers from Monnier et al. (2016) are applied to "1Dlike"-2D composite meshes of river networks using a numerical flux splitting method and an effective friction 130 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.

Multi-D hydraulic-hydrological modeling principle
The flow model consists in a spatially-distributed modeling of hydrological responses coupled to, i.e. inflowing, a seamless multi-scale "1Dlike"-2D hydraulic model of a river network, within a 2D river basin domain denoted Ω ⊂ R 2 . The core idea of this work is to numerically solve the 2D shallow water hydraulic model on a multi-D discretization D hy ⊂ R 2 of 140 the computational hydraulic domain Ω hy ⊂ Ω . This discretization (mesh) D hy being composed of N mixed unstructured triangular/quadrangular cells with internal interfaces between 1Dlike and 2D zones (see Subsection 2.2.3).
Let the basin domain Ω be composed of a hydrological domain Ω rr connected to a hydraulic domain Ω hy (Fig. 1) such that (Ω rr ∪ Ω hy ) ⊆ Ω ; their borders being respectively 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 , 145 contains interfaces with the border Γ rr of the hydrological domain Ω rr . Let us also consider D rr , a 2D discretization of the hydrological domain Ω rr .
Finally, an unstructured lattice D = D rr ∪ D 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 inflowing the hydraulic model at the interface Γ rh , as detailed in this section along with the inverse algorithm applied to the whole 150 hydrological-hydraulic numerical chain.

Hydraulic module
Numerical hydraulic models describing open channel flows generally rely on the resolution of cross sectionally or depth integrated flow equations, respectively the 1D Saint-Venant or 2D SWE (see e.g. Guinot (2010)). While 1D hydraulic models with h the water depth [m] and u = (u, v) T the depth-averaged velocity [m/s] being the flow state variables. g is the gravity 165 magnitude m/s 2 , b the bed elevation [m] and n the Manning-Strickler friction coefficient s/m 1/3 being the flow model parameters. U is the vector of state variables and F(U) (resp. G(U)) is its flux over x (resp. 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. (2016Monnier et al. ( , 2019 for details).
An effective friction law consisting in a simple power-law n = αh β is introduced, as previously done for 1D SWE for 170 effective modeling with simplified multi-channel river geometry in Garambois et al. (2017); Brisset et al. (2018).

Building-up equivalencies between 2D and 1D flow states
1Dlike 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 cross-section (XS) would. This leads to a series of quadrangular cells, each linked to a single upstream and downstream cell.

175
The 1Dlike approach implicitly assumes a rectangular XS shape which potentially impacts the representation of: (i) at a section hydraulic geometry (Leopold and Maddock (1953)), (ii) longitudinal hydraulic controls and flow variabilities.
In view to put the multi-D 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.

180
Let us consider a reference 1D model in (A, Q) variables with the bankfull width value W 1D . The friction term reads In the corresponding 1Dlike model in (h, u, v) variables, see Section 2.2.1, the friction term reads: S f,1Dlike = n 2 ||u|| h 1/3 u Considering 1D flow states over an idealized river section (Fig. 2, left) the hypothesis of local flow equilibrium (uniform, steady-state) with identical wetted areas A and WS widths W , the continuity condition implies that: Variation of the hydraulic radius R h (h) for 3 XS shapes (of similar dimensions). This showcases the potential over-and under-estimation of state variables using "1Dlike equivalent friction" from Eq. (2).
With the additional assumption of a rectangular XS (as it will be assumed in some test cases), we have h 1Dlike = A/W 1D , which leads to n 1Dlike = n 1D 190 This "1Dlike equivalent friction" leads to a perfect fit in WS elevation of 1Dlike model and 1D model in a straight prismatic channel at the given uniform regime (results not shown here for brevity). Fig. 2, right, shows the evolution of the ratio h/R h vs h. For rectangular and parabolic XS, the ratio n 1Dlike /n 1D is expected to increase with h. Thus, one can naturally expect an overestimation (resp. underestimation) of the actual friction coefficient by n 1Dlike at lower flows (resp. greater flows). This would lead to an overestimation (resp. underestimation) of the 1D WS elevation by the "equivalent" 1Dlike model. However, 195 later it will be considered a power-law in h to model the friction coefficient (Subsection 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 Subsection 3.2.2.

Multi-dimensional hydraulic model
Over a given cell K ∈ Ω hy of area m K , the piecewise constant values U K = 1 m K K UdK are approximated. Recall that the 200 finite volume 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 ) writes as follows: where U n K and U n+1 K are the piecewise constant approximations of U = (h, hu, hv) T at time t n and t n+1 (with t n+1 = t n + ∆t n ), F e stands for Riemann fluxes through each edge e of the border ∂K of the cell K, with each adjacent cell Ke. The 205 length of edge e is m e and n e,K is the unit normal to e oriented from K to Ke.
The finite volume schemes are those developped in Couderc et al. (2013); Monnier et al. (2016).The discretization of the friction source term is described in Appendix A. 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 multi-D 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 Subsection 3.3.2 for a real-like 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 215 length m ei . In the end, the flux crossing the interface e = ∪e i is equal to the sum of the fluxes crossing the e i interfaces: 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 second order accuracy in space. This solver, developed in Couderc et al. (2013);Monnier et al. (2016), is accurate and robust for wet-dry front propagations and fully applies in the present context. Note that the lateral 220 distribution of variables accross the 1D2D interface is not constrained. The source code and synthetic cases are available upon simple request.

Hydrological module
In order to simulate the hydrological response of upstream and lateral sub-catchments inflowing the river network within a 225 river basin, a hydrological module is coupled to the 2D SW hydraulic model, in a semi-distributed manner for simplicity here.
.C] ⊂ Ω rr , which outlets coordinates are respectively x i ∈ Ω rr , i ∈ [1..C]. For a given sub-catchment, the hydrological model is a dynamic operator, consisting in coupled state equations (ordinary differential equations, ODE) and output equations, that reads: Where h (x, t) is the state variables vector, P (x, t) and E (x, t) are the observable rainfall and evapotranspiration inputs (here spatially averaged over area of each subcatchment i for semi distributed 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" 235 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 semi-distributed manner in the DassFlow framework.
The model is composed of two non-linear 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 Ordinary Differential Equations (ODE) with explicit dependency to parameters, this hydrological model (Eq. 4, see detailed GR4 equations in 240 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 are the evapotranspiration and precipitation E n and P n , the output is discharge q (t) = Q r + Q d [mm/h]. E n and P n are classically imposed from data time series as piecewise constant on fixed temporal resolution 245 Figure 4. Evolution of GR4 inputs, output and reservoirs states during a sample rain event. Top: temporal forcings (rain and evaporation) and modeled output (discharge). Bottom: hydrological model states (reservoir levels).
(e.g. hourly). The numerical resolution is achieved with an implicit Euler algorithm with an adaptative sub-step algorithm enabling to reduce numerical errors especially for high flows (Santos et al. (2018)). The initial states of the stores is given by a 1 year warm-up run. The discharge q is injected into Ω hy at a sub-interface of Γ hy−rr , either as an upstream or lateral flow.
The calibrated parameters will be the classical 4 parameters θ rr = (c k,i ) k∈1..4,i∈1..C of GR4, i.e the reservoir capacities in mm for each subcatchment i ∈ 1..C. They will constitute the parameter vector θ rr considered in the forthcoming VDA 250 experiments. Other parameters, such as several drainage law exponents or the number of Nash cascade stores are not optimized as classically done with GR4 model Perrin et al. (2003). They are set at values from Santos et al. (2018).

Inverse algorithm: Variational Data Assimilation
Given spatio-temporal 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-255 hydraulic chain composed of a hydraulic model, presented in Section 2.2, and a hydrological model, presented in Section 2.3.
Detailed know-hows on VDA may be found in online courses (see e.g. Bouttier and Courtier (2002); 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 multi-D with GR4 modules is obtained by automatic differentiation with Tapenade engine Hascoet and Pascual (2013). Note that 2D, 1Dlike and 1D2D 260 models are functionally identical from the point of view of the adjoint code and the assimilation process.
We consider the following set of spatio-temporal observations of water surface and discharge over the river domain Ω ⊂ R 2 : with Z o the observed WS elevation [m] above reference elevation, Q o the observed discharge m 3 /s , N o,Z the number of altimetric observations points and N o,Q the number of observed discharges over Ω.

265
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 sub-catchment here. Also note that other water surface observables could be considered, such as water surface velocity observations or dynamic water masks -not the scope of this study.
Given river stage and/or discharge observations, the aim here is to estimate unknown or uncertain quantities of the hydrologicalhydraulic model among: discharge hydrographs Q i (t) , i ∈ [1, N ] on the border of the hydraulic domain (in 1Dlike and/or 2D 270 parts), spatially distributed hydraulic parameters over a river network (bathymetry elevation b or friction n) and hydrological 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 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,  Moreover, we classically enrich the cost function with a regularization term: j (θ) = j obs (θ)+γj reg (θ) with j reg a Tikhonov type regularization term. Here we consider a regularization on the bathymetry only: 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 low pass filtering of the bathymetry shape, see Martin and 290 Monnier (2015); 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 of variables is applied: with B the covariance matrix of the background error.

295
Then by setting J (k) = j (θ), the optimization problem which is solved is actually the following: The first order optimality condition of this optimization problem (Eq. 9) reads B 1/2 ∇j (θ) = 0. The change of variables based on the covariance matrix B acts as a preconditioning of the optimization problem. This optimization problem is solved using a first order gradient-based algorithm, more precisely the classical L-BFGS quasi-Newton algorithm (limited-memory The choice of the covariance matrix B, represents an important a priori 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 B Ω 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 .
The following stopping criterion 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 315 value. The multi-D hydrological-hydraulic tool chain presented above has been implemented in the latest version of Monnier et al. (2019).
3 Results and discussion

Numerical experiments design
Both synthetic and real cases are considered to test the forward and inverse modeling capabilities of the proposed computational 320 chain for river networks and floodplains simulation.
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. inference, is assessed against 325 a reference 2D model built on fine bathymetry of 75 km of the Garonne river and (ii) the whole multi-D hydraulic hydrological tool chain is tested at basin scale on the real complex case of the Adour River network. Then, the model is tested on two real cases. The 1Dlike model is tested against a reference 2D model built on fine bathymetry of 75 km of the Garonne river. The whole multi-D hydraulic-hydrological tool chain is tested on the Adour River network.

Synthetic cases 330
First, the proposed multi-D hydraulic solver (Subsection 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 contain frontal interfaces between 1Dlike and 2D meshes. Next, to investigate the reproductibility 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 i (t) , i ∈ [1..N ] and of a friction power-law n = αh β are carried out using WS observables.

Forward multi-D hydraulic cases
Straight channel case A prismatic rectangular channel and a multi-D 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 multi-D model is compared to a reference 2D model with refined mesh (mesh not shown, 2400 cells, average edge length ∽ 25 m). The multi-D waterline is validated 340 against the 2D model at permanent flow Q = 100 m 3 /s and the modeled downstream discharges are compared for a flood hydrograph ( Fig. 6(b)). Both first (not shown) and second order numerical solver allow a close fit to the target water line (Fig.   6, bottom).
At permanent flow, a slight misfit is observed between the 2D and multi-D WS elevations with the second order scheme ( Fig. 6, top, relative misfit < 0.15% at 1000 m). This is due to approximation error in the multi-D model caused by large spatial 345 steps in 1Dlike reaches (dx = 200 m). Indeed, this misfit is reduced by reducing the 1Dlike cells 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. 6, top). This is due to the second order scheme, which is currently not designed for multi-D interfaces. Recall that no constraint is imposed on the lateral distribution of computed variables. The technical implementation of this reconstruction for 1D2D 350 interfaces will be done in next version of DassFlow.
During a varied flow event, the outflow of both models is close to identical (Fig. 6, middle). The flow is correctly transmitted at multi-D interfaces and, at this scale, the 1Dlike meshes are adequate to model a flood wave propagation.
Simple confluence case   Station locations are noted as black dots. of synthetic cases: (i) a straight rectangular channel, (ii) a straight rectangular channel with a mid-channel slope break and (iii) a straight parabolic channel.

Direct calibration
Recall the equivalent friction formulation (Eq. 2) designed to make match the water line at-a-section and at equilibrium. In a straight rectangular prismatic channel, with the 1D normal water depth imposed downstream ( Fig. 8(a)), a backwater curve is computed by the DassFlow1D model (in red). With a 1D homogeneous friction (n = 0.05 s/m 1/3 ), the 1Dlike approach yields an underestimated waterline (in blue). An homogeneous effective friction allows 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 380 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. 8(b)). 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.

385
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 an homogeneous shift δ b of the reference bathymetry, are needed to match the 1D WS elevation. Equivalent friction only allows to model identical wetted sections at permanent bankfull flow for a given channel width (Section 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.

390
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 (Subsection 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.

395
In the following comparison, power-law 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 1/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 days). 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 400 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.  to infer Q 1 and Q 2 close to perfectly from configurations (a) and (b) (Fig. 10, right).

Hydrological parameters inference
In this second twin experiment, the issue of 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 GR4H module applied on two upstream catchments inflowing the hydraulic module. For each catchment, synthetic rain and evaporation time series and a hydrological parameter set θ rr = (c 1 , ..., c 4 ) are used to generate a discharge Since the observed variables are of different nature and amplitude, we introduce a normalization. The cost function is here j obs = j Z + j Q , with each term being normalized by the number of observations and by their range of variation such that: Those correspond to two separate normalized squared RMSE with N o,Q and N o,Z denoting the number of observations station and T Q and T Z the number of observation time steps.
The control vector c = (θ rr1 , θ rr2 ) contains the two sets of 4 hydrological parameters each. For this synthetic case, an 430 inequality contraint of the control parameters is imposed with the bounded L-BFGS-B 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 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 435 close to the target value starting from an erroneous prior (Fig. 11, right). 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.  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). 1Dlike cell interfaces are perpendicular to the flow direction, as would be XSs in a 1D model. Each cell is about 100 m in longitudinal length. This mesh was generated using the SMS 3 meshing tools.
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 an homogeneous friction parameter of 0.05 s/m 1/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 455 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 (Subsection 2.2.2).

Calibration by hand
We here propose to reduce the misfit using, on one hand, effective friction and, on another 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 1/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 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 480 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 8.6 × 10 −10 ), while for Model B1b it is reached in 85 iterations (for a normalized cost of 5.5 × 10 −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 .

485
Variational calibration of channel parameter has allowed to fit a permanent regime water line. The following paragraphs study the reproduction of propagation in a calibrated models 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 last 10 days, with a peak discharge of 490 702 m 3 /s at the start of day 2, a minimum discharge of 160 m 3 /s and and average discharge of 382 m 3 /s. 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 Using Model A1 as a prior for bathymetry and friction and observations from the 2D model, we find the following optimal parameters: n = 0.054 s/m 1/3 , δ b = 0.669 m . The resulting optimal model is Model C. To allow a better fit during high and 500 low flows, we introduce the friction power-law n = αh β (Subsection 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 1/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, Compared to Model C (average misfit of 0.086 m and 0.20 m at resp. 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 an homogeneous friction value n in a 1D-like model has allowed to closely fit WS elevation observations of a 10 days flood wave over the reference Garonne model. The calibration 510 of depth-dependent 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.

Adour basin: multi-D hydrological-hydraulic model
The whole hydrological-multi-D hydraulic tool chain is now tested in a real context both for forward and inverse problems resolution. A real and complex basin case is now considered: the Adour river basin in the South West of France, with a total 515 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.
In this section, a multi-D model, composed of 1Dlike meshes and 2D zooms over 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 multi-D hydraulic network model and several hydrological models of sub-catchments (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  Table 3. Garonne models parameters and metrics. "PR" stand for Permanent Regime (discharge of 400 m 3 /s). The "Varied" flow event is non-flooding and has a mean discharge of 398 m 3 /s, with a minimum of 183 m 3 /s and a maximum of 702 m 3 /s "High" observation density corresponds to 720 stations, or 1 station per 100 m. "Low" observation density corresponds to 72 stations, or 1 station per 1 km.
45 km), the Oloron and Pau rivers (around 65 km in total). The river networks contains mostly single-branch 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, red points in Fig. 15). Out of the 5 remaining stations kept for data assimilation, 3 are located on river reaches (Peyrehorade, Urt and Villefranque, blue points in Fig. 15) and 2 are located in the Bayonne area (Pont-Blanc and Lesseps, blue points in Fig. 15).

530
At the 4 upstream points of the modeled river network, 4 sub-catchments are modeled with 4 lumped hydrological models (Subsection 2.3). Their drainage areas are respectively about 7811, 842, 2480 and 2464 km 2 . The hourly discharge have been extracted from the HYDRO 4 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 is provided by the SAFRAN reanalysis, Météo France, and then used to calculate the potential evapotranspiration using the Oudin formula 535 Oudin et al. (2005). The rainfall and potential evapotranspiration are at a spatial resolution of 1 km 2 square grid, and processed into hourly time step. Spatial averages of the rainfall and potential evapotranspiration computed using the SMASH distributed hydrological modeling platform (Jay-Allemand et al. (2020); Colleoni et al. (2021)) over the 4 catchments are used as inputs for the lumped GR4H model. Lumped parameter sets for the 4 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 hydro-graphs, not hydrological parameters. Indeed, the study of global calibration and regionalization issues of spatially distributed hydrological models is left for further research.
Our multi-D hydraulic modeling approach (Section 2) is applied to this complex case as follows. First, a "1Dlike-only" model of the whole network is built. Then, a multi-D model is built based on the "1Dlike-only" model, with the addition of a 2D mesh of a floodplain (Fig. 16, left).

545
On the "1Dlike-only" 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 on a large scale river network. This "1Dlike-only" model contains 1409 cells. A 24h simulation runs in 550 around 13s, using uncalibrated parameter estimates.
Then, a multi-D model is obtained: the existing 2D mesh from a Telemac model is coupled to a 1Dlike mesh, similarly to the reference Telemac-Mascaret model from the regional flood forecast center SPC-GAD. The 1Dlike parts of the mesh are kept identical to the "1Dlike-only" model, while the 2D part is the mesh extracted from the Telemac-Mascaret model (Fig. 16 This relatively low computation time, and the potential to decrease it further by using more threads in parallel, indicate that this multi-D method is suited to operational use. The code version this work is based on was proven scalable (Couderc et al. (2013)). Additions made to the current version should not change this, but no numerical testing 565 has been done.

Assimilation of WS observables to infer 4 upstream hydrographs
To investigate the 1Dlike effective modeling on a river network, a twin experiment setup is designed to infer large control vector from a realistic observability on the "1Dlike-only" model.

570
The considered control vector is composed of the 4 upstream hydrographs c = (Q Dax (t) , Q Esc (t) , Q Ort (t) , Q Cam (t)).
Observations of WS elevation are generated at Peyrehorade, Urt and Villefranque stations, and additionally at a virtual station directly downstream from Orthez. These 4 points give theoretically sufficient observability to identify the 4 upstream hydrographs. Indeed, they sample mixed and unmixed signal similarly to the academic setup in Subsection 3.2.1). The observation pattern also corresponds to a reasonable expectation of spatial observability in French river networks. Prior hydrograph values 575 are classically set to a constant average discharge value.
Simultaneous inference of the 4 hydrographs is satisfying. As shown in Fig. 17, left, upstream hydrographs injected at Cambo and Orthez are inferred very accurately (N SE > 0.95). This is due to their WS elevation signature being observed "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 partial error of signature attribution and a less accurate inference 580 (N SE = 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 N SE of 0.50 and closer in shape to the inferred Escos and Orthez hydrographs. This hints at correlated influences of these hydrographs at observation stations and insufficient observability of the Dax hydrograph signal given the model hydrodynamics. Observed signals at the 4 stations are however all very accurate (Fig. 17, right). For upstream stations (Orthez and Villefranque), this is due to the accurate inference of upstream 585 hydrographs. For stations under the influence of the tidal 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 on 1Dlike network models. Note that multi-D models are identical to the investigated 1Dlike model in terms of VDA capabilities. It paves the way towards investigations on multi-variate inferability of spatially distributed hydraulic-hydrological parameters.

590
In the future, investigation on the influence of data sparsity, observation weights and ill-posed problem constraints should be carried out.

Conclusions and perspectives
This article presents a new approach and numerical chain for the multi-D hydrological-hydraulic modeling of complex river networks with variational data assimilation capabilities. It is based on the VDA algorithm and the finite volume solvers (in-595 cluding second order one and accurate treatment of wet/dry front propagation) from Monnier et al. (2016). The resolution of the full 2D shallow water equations (Eq. 1) is performed with a single finite volume solver applied on a multi-D 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 higher resolution unstructured meshes -either triangular or quadrangular. This hydraulic model is inflowed with a hydrological model enabling to describe upstream/lateral catchments inflow hydrographs.

600
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 ). The adjoint of the whole tool chain is generated with Tapenade (Hascoet and Pascual (2013)) and validated.
Forward-inverse capabilities with the new components is assessed on several cases of increasing complexity. RMSE value given are used to calculate cost in the assimilation process (Section 2.4).
The 1Dlike effective hydraulic modeling approach as well as its coupling with higher resolution 2D zooms was validated, 605 against: 1) reference 1D and 2D hydraulic models, on an array of academic; 2) complex real cases featuring 1D and 2D flow variabilities in river networks with confluences and floodplains. Those cases are also employed to test the coupling with the hydrological model implemented in a semi-distributed setup.
Considering single (resp. multiple) type(s) of flow signatures observations in those river networks through a single-(resp. multi)-objective observation cost function, the capabilities of the VDA method for inferring mono/multi-variate control vectors 610 of large dimensions was successfully tested.
From the obtained results, the following conclusions can be raised: 1. A complete integrated multi-D 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.
2. The 1Dlike modeling approach enables to simulate fine physical flow states compared to reference 1D or 2D SW models; 615 hydrograph propagation remains very close also.
3. The inference, from heterogeneous observations in the river network, of multi-variate and high dimension controls among multiple inflow discharge hydrographs, bathymetry, friction of the full shallow water multi-D 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 620 discussion in Brisset et al. (2018);Larnier et al. (2020); Garambois et al. (2020) and references therein).
4. Information feedback from the river network to upstream hydrological models of sub-catchments is shown.
5. 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 geometry-friction patterns. The depth-dependent friction law helps to reduce misfits across flow regimes. To our knowledge, the present numerical tool is the first one proposing, large scale multi-D river network modeling with VDA capabilities. This new toolchain opens the way for the resolution of large scale high dimensional hydrological-hydraulic inverse problems that can be considered given constraints from multi-source datasets. The methods for direct and inverse modeling 630 can be applied at multiple spatial scales including with fine resolution (imposing finer calculation time steps to respect CFL condition in the forward hydraulic model). Building and constraining the models is, however, dependant on data availability and the informative content of observations, which may be linked to the spatial scale.
Short term perspectives will aim to taylor the data assimilation algorithm to perform complex data assimilation experiments at basin scale using various multi-source datasets. To be actually operational, improvements pertaining to the construction of 635 1Dlike models from global public databases are needed to deploy the multi-D approach to a large number of river networks.
Coupling is ongoing with the SMASH spatially distributed hydrological platform (Jay-Allemand et al. (2020); Colleoni et al.
(2021)) on which is based the French flash flood warning system VigiCrues Flash (Garandeau et al. (2018)). Note that the addition of lateral flows (surface and subsurface runoff) would simply consist in imposing additional inflows on the boundary of river/floodplains domain, which is numerically identical to the upstream flows imposed and inferred here (see such lateral cou-640 pling in 1D with inferences from sparse satellite altimetry data in Pujol et al. (2020)). Further work will also test the integrated chain on flash flood Mediterranean basins as well as 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 reaches parameterizations. This is especially true with depth-dependent porosity (Özgen et al. (2017); Guinot et al. (2018) and references therein) applied to complex channel geometries and with spatially distributed calibration. The HLLC approximate Riemann solver is used. This gives the expressions: where the wave speed expression is those proposed in Vila (1986): It has been demonstrated that this insures L ∞ stability, positivity and consistency with entropy condition under a CFL 655 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).

660
The numerical scheme must preserve the fluid at rest property, that is the gradient of bathymetry ∇z b must not provide u n+1 ̸ = 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: ∇ gh The technique which is employed here is those presented in Audusse et al. (2004); Audusse and Bristeau (2005). It is based 665 on the following change of variable.
From now, at left and right sides of edge e, the considered water depth values h * ⊓ are those defined from the "reconstructed" topography z e as: The conservative variable vector U n e,K is now considered with the new variable: Note that this new variable h * e,K depends on the bathymetry values z e,K and z e . The resulting scheme reads: This provides a well-balanced (first order) scheme. Figure A1. Typical situation of desired well-balanced property in presence of a wet/dry front. The WS elevation is here denoted by η.
A3 Prediction-correction time scheme

680
The friction source term is taken into account in the complete SW system by deriving a prediction-correction time scheme, see e.g. Toro (2001).
We denote here byŪ n+1 K the FV solution at time t n+1 of the (well-balanced) 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: U = (h, hu) T .
At each time step, from n to n + 1, 685 -Step 1: computation ofŪ n+1 , solution of the conservative SW system, i.e. the SW system without S f , i.e. the FV solution of the following system: Step 2: given the "predicted value"Ū n+1 , compute U n+1 solution of: General schemes (explicit, implicit or semi-implicit ones) including the friction source term S f in the discretization of the model (Eq. 1) can be written as: for all K, Note that this splitting scheme is consistent at first order in ∆t with the complete SWE. Splitting scheme second order in time is possible; it is not detailed later.

A4 Second order MUSCL scheme
In order to obtain a globally second order scheme, a higher-order time stepping scheme is needed. Let us briefly describe the ingredients of this second order well-balanced positive scheme that is strictly the same as the one proposed in Couderc et al. (2013); Monnier et al. (2016). Actual second order accuracy, considering source terms S g and S f , is achieved through the 720 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 semi-implicitation friction source term S f given by in the subsection above.
It leads to new expressions of U n K and U n Ke . With this linear reconstruction, one can expect a scheme with a second-order 725 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 Finite Volume mass conservation property, a Barth limiter (Barth (2003)) is employed. The friction term S f is classically parameterized with the empirical Manning-Strickler law established for uniform flows |Q| Q K 2 A 2 R 4/3 h where K m 1/3 /s is the Strickler coefficient.

740
The Saint-Venant equations are solved on each segment of the river network and the continuity of the flow between segments is ensured by applying an equality constrain 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; lateral hydrographs q l (t) at in/outflow nodes. The initial condition is set as the steady state backwater curve profile Z 0 (x) = Z (Q in (t 0 ) , q l,1..L (t 0 )) for hot-start. This 1D Saint-Venant 745 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 to deal with flow regimes changes with a one-hour time step ∆t. This is implemented into the computational software DassFlow (Brisset et al. (2018);Larnier (2010)).
The numerical scheme is a semi-implicit finite difference scheme (generalized Preissmann scheme) with a double sweep Local Partial Inertial method to minimize the inertial terms (see documentation in Brisset et al. (2018);Larnier (2010)).
Calibrated parameter for the Adour case were obtained using the airGR global calibration algorithm (Coron et al. (2017)) from the freely available package 5 .  Table C1. Calibrated hydrological parameters of the 4 upstream hydrological catchments from the Adour multi-D hydrographic network model. *Note that the x4 calibrated parameters corresponds to the non-"state-space" GR4H version (not presented, seePerrin et al. (2003)), for which the calibration tool is provided. x4 values in the present run with the "state-space" were set at 0.15.
LP with analyses by all. JM has designed or supervised the 2D numerical scheme and VDA algorithm. PAG and JM have designed and co-supervised the 1D2D algorithm. All authors have participated to the writing. Project management for funding: PAG.
Competing interests. The authors declare that they have no competing interests Acknowledgements. This study was rendered possible by the adaptation of covariance matrices in the inverse algorithm from the DassFlow1D Univ. Montpellier) for fruitful discussions related to hydraulic solvers. They warmly thank Laurent Dieval from SPC-GAD 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, Univ. of Strasbourg), the PhD director of LP.
Financial support. PhD of LP has been co-funded by CNES and ICUBE.