L Pujol
email: l.pujol@unistra.fr

Monnier, J P-A Garambois
email: pierre-andre.garambois@inrae.fr

towards multisource data fusion from river network to headwaters INTEGRATED HYDRAULIC-HYDROLOGICAL ASSIMILATION CHAIN: TOWARDS MULTISOURCE DATA FUSION FROM RIVER NETWORK TO HEADWATERS

Keywords: Shallow water, Hydrological coupling, 1D-2D seamless model, Variational Assimilation, Meshing strategies

 and references therein) and complex river network portions, forced by spatially distributed inflows (Pujol et al. (2020), Malou et al. (under redaction)). Since hydraulic modeling at the scale of a river basin can be computationally costly, a combination of effective 1D and 2D representations, complemented by hydrological modules, may be useful. Complex river-floodplain interaction zones may be modeled in 2D zooms, while 1D approaches can fit simpler reaches. This contribution presents the development of a complete hydraulichydrological toolchain based on the 2D hydraulic model and variational data assimilation platform DassFlow 2 . A 1D effective modeling approach based on a 2D shallow water model is tested. Then, the implementation of hydrological modules within the DassFlow VDA framework is presented.

INTRODUCTION

In a context of climate change and potential intensification of the hydrological cycle, improving the understanding and representation of water fluxes within river basins is of paramount importance both for hydrological sciences and operational forecasts. It is critical to improve our capacity to model hydrological responses to meteorological variabilities and their hydraulic consequences. Cascading from rainfall to floodplain dynamics is crucial for flood-inundation hazard simulation and several studies have combined hydrological and hydraulic modeling components (e.g. [START_REF] Bonnifait | Distributed hydrologic and hydraulic modelling with radar rainfall input: Reconstruction of the 8--9 September 2002 catastrophic flood event in the Gard region, France[END_REF][START_REF] Pappenberger | Cascading model uncertainty from medium range weather forecasts (10 days) through a rainfall-runoff model to flood inundation predictions within the European Flood Forecasting System (EFFS)[END_REF][START_REF] Uhe | Model cascade from meteorological drivers to river flood hazard: flood-cascade v1.0[END_REF] among others).

Furthermore, the increasing availability of high resolution remotely sensed data, especially global terrain elevation and rivers morphologies, have enabled spatially distributed models to be applied at the regional or global scale (see [START_REF] Paiva | Large scale hydrologic and hydrodynamic modeling using limited data and a GIS based approach[END_REF][START_REF] Yamazaki | A physically based description of floodplain inundation dynamics in a global river routing model[END_REF] among others). Recently, [START_REF] Fleischmann | Modeling the role of reservoirs versus floodplains on large-scale river hydrodynamics[END_REF] proposed a 2D regional model based on simplified channel-floodplain geometry representation on a structured grid and validated using satellite data. They used a hydraulic conceptualization based on the subgrid model developed by [START_REF] Neal | A subgrid channel model for simulating river hydraulics and floodplain inundation over large and data sparse areas[END_REF] for large scale 2D flood modeling. Those 2D approaches, applicable at relatively large scale, consist in more or less simplified hydrological and hydraulic models for the sake of applicability on large computational domains with 1 Corresponding author 2 http://www.math.univ-toulouse.fr/DassFlow SimHydro 2021: Models for complex and global water issues -Practices and expectations [START_REF] Miglio | Model coupling techniques for free-surface flow problems: Part I[END_REF][START_REF] Miglio | Model coupling techniques for free-surface flow problems: Part II[END_REF][START_REF] Monnier | Inverse algorithms for 2D shallow water equations in presence of wet dry fronts. Application to flood plain dynamics[END_REF] June 2021, Sophia Antipolis -Pujol, L., Garambois, P-A., Monnier, J., Finaud-Guyot, P., Larnier, K., Mosé, R. -Integrated hydraulic-hydrological assimilation chain : towards multisource data fusion from river network to headwaters calibration datasets of variable sparsity and information content. Such models may be insufficient for capturing fine scale but crucial hydrodynamic phenomena, like fine streamflow concurrency on a river network and wet-dry fronts propagation on floodplains.

Fine 2D hydraulic models, resolving the full shallow water equations (SWE), on entire river networks and floodplains, i.e. relatively large spatio-temporal domains, can be computationally costly. To overcome this issue, coupled 1D-2D models have been developed, where the 2D model is applied only where complex hydraulic behaviors need to be reproduced, while the less costly but physically sound 1D approach is applied to simpler reaches. Coupling methods usually consist in minimizing the discrepancy between subdomain variables at their, sometimes overlapping, borders. This can be achieved using decomposition methods controls (see [START_REF] Gervasio | Heterogeneous coupling by virtual control methods[END_REF][START_REF] Miglio | Model coupling techniques for free-surface flow problems: Part I[END_REF][START_REF] Miglio | Model coupling techniques for free-surface flow problems: Part II[END_REF]) or superposition methods (see [8,[START_REF] Marin | Superposition of local zoom models and simultaneous calibration for 1D--2D shallow water flows[END_REF]). In both cases, the key is to ensure the hydraulic coherence of the coupled models at their junctions, with regard to the transmitted hydraulic variables, either a priori via a physical analysis or a posteriori using a mathematical analysis of error estimates. 2D zooms are notably used in [8,[START_REF] Marin | Superposition of local zoom models and simultaneous calibration for 1D--2D shallow water flows[END_REF], where a superposed grid is used for coupling the models in view of using data assimilation methods on the 2D grid. Another approach to 1D-2D modeling, following [START_REF] Finaud-Guyot | 1D/2D porosity model for urban flood modeling: case of a dense street networks, 40 : 06010. SimHydro 2021: Models for complex and global water issues -Practices and expectations 16-18[END_REF], consists in an effective "seamless" 1D-2D model obtained directly from simple flux splitting in a finite volume solver. This is achieved here using a finite volume solver, adapted to accurate floodplain flows simulation including wet-dry fronts (see [START_REF] Monnier | Inverse algorithms for 2D shallow water equations in presence of wet dry fronts. Application to flood plain dynamics[END_REF]).

Each component of a hydraulic-hydrological chain is subject to uncertainties (inputs, model structure, parameters). Improving the accuracy of simulations requires to take advantage of information available in multi-source observations (in situ, satellite, opportunistic). Hydrological models are generally calibrated on streamflow data and a few studies started to investigate synergies of those coupled models with remotely sensed data (e.g. [START_REF] Pujol | Estimation of Multiple Inflows and Effective Channel by Assimilation of Multi-satellite Hydraulic Signatures: Case of the Ungauged Anabranching Negro River[END_REF][START_REF] Wongchuig-Correa | Assimilation of future SWOT-based river elevations, surface extent observations and discharge estimations into uncertain global hydrological models[END_REF] among others). They use data assimilation, an adequate technique to produce numerical representation that optimally combine models and data. This technique allows to estimate spatially distributed channel parameters (bathymetry, friction) and/or inflows for example using flow observables (water surface elevation, discharge), see [START_REF] Larnier | On the estimation of river discharges from altimetry[END_REF][START_REF] Pujol | Estimation of Multiple Inflows and Effective Channel by Assimilation of Multi-satellite Hydraulic Signatures: Case of the Ungauged Anabranching Negro River[END_REF]. The increasing observability of river networks from multisensors requires new methods for multi-source data assimilation over full river network models. More precisely, integrated and spatially distributed hydraulic-hydrological modeling approaches that are both accurate and applicable over large domains as well as capable of exploring synergies with multi-sourced data are required. The assimilation of multi-sourced observations in a 1D-2D seamless hydraulic-hydrological model has not been studied yet, to our knowledge. This paper presents the components of an integrated hydraulic-hydrological modeling chain enabling variational data assimilation based on the DassFlow 2D platform (see [START_REF] Monnier | Inverse algorithms for 2D shallow water equations in presence of wet dry fronts. Application to flood plain dynamics[END_REF]) and a simple and widely used conceptual hydrological model GR4H [START_REF] Perrin | A conceptual model of flow and transport in a karst aquifer based on spatial and temporal variations of natural tracers[END_REF] in its state-space version (see [START_REF] Santos | Continuous state-space representation of a bucket-type rainfall-runoff model: a case study with the GR4 model using state-space GR4 (version 1.0)[END_REF]). An effective 1D-2D modeling strategy is implemented as well as hydrological modeling components. Forward and inverse problems are investigated on real and academic cases. In subsection 2.1, a method for the effective modeling of 1D flows using a 2D model, based on [START_REF] Finaud-Guyot | 1D/2D porosity model for urban flood modeling: case of a dense street networks, 40 : 06010. SimHydro 2021: Models for complex and global water issues -Practices and expectations 16-18[END_REF], is presented ; in subsection 2.2, the considered hydrological model is described ; and in subsection 2.3, the inverse method is described. Section 3 presents a series of academic and real cases which study in detail the different parts of the proposed hydraulichydrological modeling method for large scale river networks.

Authors affiliation

Pujol, L.: Laboratoire des Sciences de l'ingenieur, de l'informatique et de l'imagerie (ICUBE), Fluid Mechanics Team, CNRS, Universite de Strasbourg, France Garambois, P-A.: INRAE (Irstea), Aix Marseille Univ, RECOVER, Aix-en-Provence, France Monnier, J.: Institut de Mathematiques de Toulouse (IMT), France Finaud-Guyot, P.: HSM, Univ Montpellier, CNRS, IRD, Montpellier, France Larnier, K. : CS corporation, Business Unit Espace, Toulouse, France Mosé, R.: Laboratoire des Sciences de l'ingenieur, de l'informatique et de l'imagerie (ICUBE), Fluid Mechanics Team, CNRS, Universite de Strasbourg, France SimHydro 2021: Models for complex and global water issues -Practices and expectations [START_REF] Miglio | Model coupling techniques for free-surface flow problems: Part I[END_REF][START_REF] Miglio | Model coupling techniques for free-surface flow problems: Part II[END_REF][START_REF] Monnier | Inverse algorithms for 2D shallow water equations in presence of wet dry fronts. Application to flood plain dynamics[END_REF] June 2021, Sophia Antipolis -Pujol, L., Garambois, P-A., Monnier, J., Finaud-Guyot, P., Larnier, K., Mosé, R. -Integrated hydraulic-hydrological assimilation chain : towards multisource data fusion from river network to headwaters

MODELS AND METHODS

The components of the integrated hydraulic-hydrological modeling chain are described in this section and consist in (see schematic catchment and hydrological-hydraulic domains representations in Fig. 1): (i) a fine 2D hydraulic model adapted to high quality zooms on complex flow zones (ex. confluences, floodplains), (ii) an effective 1D model for stream flows over river network with seamless connexion between effective 1D and 2D zones, (iii) a hydrological model for distributed rainfall-runoff inputs into the hydraulic domain, (iv) a computational inverse method enabling to optimize parameters of the whole modeling chain with multisource data and enable information feedback between flow components. This chain is implemented in the DassFlow platform3 based on existing numerical schemes for solving the 2D shallow water equations (SWE) and corresponding inverse methods based on variational data assimilation (VDA) (see [START_REF] Monnier | Inverse algorithms for 2D shallow water equations in presence of wet dry fronts. Application to flood plain dynamics[END_REF]).

Multi-dimensional hydraulic model

The 2D SWE in conservative form with a Manning-Strickler friction term read as follows:

∂h ∂t + ∂hu ∂ x + ∂hv ∂ y = 0 ∂hu ∂t + ∂ ∂ x ( hu 2 + gh 2 2 )+ ∂huv ∂ y = -gh ∂ z b ∂ x -g n 2 ‖u‖ h 1 /3 u ∂hv ∂ t + ∂huv ∂ x + ∂ ∂ x (hv 2 + gh 2 2 ) = -gh ∂ z b ∂ y -g n 2 ‖v‖ h 1/ 3 v in Ω×] 0 , T ] , (1)
where h denotes the water depth, u and v are the depth averaged flow velocities, g is the gravity and z b is the bottom elevation on a computational domain Ω ∈ℝ 2 and for a time interval [0 ,T ] .

These SWE are numerically solved using finite volume schemes available in DassFlow: first or second order for accurate representation of wet-dry front dynamics (see [START_REF] Monnier | Inverse algorithms for 2D shallow water equations in presence of wet dry fronts. Application to flood plain dynamics[END_REF]).

Effective 1D-2D

A seamless 1D-2D approach, with cheaper effective 1D reaches connected to 2D zooms, is obtained by solving the 2D SWE (eq.( 1)) on a composite mesh including internal boundary conditions between 1D and multiple 2D cells (see [START_REF] Finaud-Guyot | 1D/2D porosity model for urban flood modeling: case of a dense street networks, 40 : 06010. SimHydro 2021: Models for complex and global water issues -Practices and expectations 16-18[END_REF] and references therein). The theoretical equivalence between 1D and 2D models is discussed in [START_REF] Michel-Dansac | Consistent section-averaged shallow water equations with bottom friction[END_REF] for trapezoidal cross-section (XS). In 1D effective reaches, the river channel is meshed using consecutive channel-wide quadrangular cells, each cell is thus connected to a single upstream and a single downstream cell. Cell edges are set perpendicular to the flow direction, thus verifying the 1D hypothesis. In this configuration, the 2D Riemann solver is forced to solve for flows on a single interface and thus is reduced to a 1D solver. Remark that this effective 1D approach implicitly assumes a rectangular channel XS. Details of this solver "degeneration" are not presented here for brieviety and are available on demand.

In reaches where 2D modeling is necessary, classical triangular cell meshes are built. These "1D_eff" and 2D meshes are linked through internal boundary conditions based on classical boundary conditions. Inflows are injected into downstream boundary cells proportionally to the cells water depths. The average WSE over a downstream boundary cell group is imposed to upstream cells. Note that for 1D_eff models presented in this paper, the choice is made of longitudinally large cells (at least 100m). Given the relatively low speeds of signal propagation, the CFL values for these models are low.

Equivalent friction

In view to propose an effective friction term for our 1D_eff modeling approach introduced above and based on the 2D SWE, a simple comparison is made with a 1D model. It solves the classical 1D Saint-Venant SWE in (S,Q) variables:

∂ S ∂ t + ∂ Q ∂ x = 0 ∂Q ∂t + ∂ ∂ x ( Q 2 S )+ gS ∂ H ∂ x = -gS n 2 Q|Q| S 2 R h 4 /3 , (2)
where S denotes the wetted area, Q is the discharge and H is the water surface elevation (WSE).

The projection of eq.( 1) on the x axis leads to an identical formulation to eq.( 2), except for the friction source term on the right hand side of eq.( 2). Assuming the equality of 1D and 1D_eff friction source terms over a cross-section of width W 1 D gives: -g

S 1 D W 1 D n 1 D 2 Q 1 D |Q 1 D | S 1 D 2 R h, 1 D 4 /3 = -g n 2 D 2 ‖h 2 D u 2 D ‖ h 2 D 1 /3 h 2 D u 2 D .
Given our meshing choices, we can assume identical flow velocities ( u 2 D =Q 1 D /S 1 D ) and identical river widths. Let us also assume, as depicted in Fig. 2, identical wet surfaces (

S 2 D =S 1 D , thus h 2 D = S 1 D W 2 D
) in a true 1D model and its equivalent 2D counterpart.

This gives n 2 D =n 1 D √ S 1 D 4 / 3 R h , 1 D 1 /3 W 2 D 1 /3 W 1 D
in the general case or, for a rectangular 1D channel:

n 2 D =n 1 D √ h 4 /3 R h 4 / 3 , (3)
The Manning friction coefficient can be written as a power law of the water depth n=ah b , where a and b are model parameters, to account for the observed depth/friction relation. This formulation may help to effectively reproduce the behavior of a real 1D case with the proposed 1D_eff model. The start of an investigation is found later in this paper.

Hydrological model

A hydrological model is considered to simulate the rainfall-runoff relationship on upstream subcatchments following the schematic catchment discretization introduced in Fig. 1. In this study, we consider the GR4H hourly hydrological model, a parsimonious lumped conceptual model with 4 free parameters (see [START_REF] Perrin | A conceptual model of flow and transport in a karst aquifer based on spatial and temporal variations of natural tracers[END_REF]) in the state-space version with a Nash cascade (differentiable) instead of a unit hydrograph proposed in [START_REF] Santos | Continuous state-space representation of a bucket-type rainfall-runoff model: a case study with the GR4 model using state-space GR4 (version 1.0)[END_REF]. It consists in a series of reservoirs including, in order, a production reservoir, a Nash cascade of 11 reservoirs that simulates a Unit Hydrograph behavior, and a routing reservoir. It sums up to a set of coupled ordinary differential equations (ODE). Their free parameters are x 1 and x 3 , defining the maximum capacity of production and routing reservoirs respectively ; x 4 , the base time of the "Unit Hydrograph" which rules the draining rates of the 11 Nash cascade reservoirs ; and x 2 , an inter-catchment exchange coefficient.

The source code from [START_REF] Santos | Continuous state-space representation of a bucket-type rainfall-runoff model: a case study with the GR4 model using state-space GR4 (version 1.0)[END_REF] has been included in DassFlow software and differentiated using TAPENADE engine (see [START_REF] Hascoët | The Tapenade Automatic Differentiation tool: Principles, Model, and Specification[END_REF]) to enable VDA as explained hereafter. The adjoint of the whole hydrological-hydraulic chain has been successfully validated with classical tests.

Inverse method

DassFlow is a computational software framework enabling Variational Data Assimilation (VDA) (see details in [START_REF]Data assimilation for free surface flows[END_REF][START_REF] Larnier | On the estimation of river discharges from altimetry[END_REF][START_REF] Monnier | Inverse algorithms for 2D shallow water equations in presence of wet dry fronts. Application to flood plain dynamics[END_REF]). The principle of this inverse method is to estimate model parameters, collected in a control vector c , through the minimization of a cost function j that measures the misfit between observations (water levels, discharges, ...) and simulated variables, which depend on the parameters through the hydrodynamic model. The minimizaiton of j is performed with a quasi Newton descent algorithm. It requires the gradient of the cost function computed by solving the adjoint model obtained by numerical differentiation with TAPENADE engine (see [START_REF] Hascoët | The Tapenade Automatic Differentiation tool: Principles, Model, and Specification[END_REF]). Two minimization algorithms are available, either unbounded (M1QN3, see [START_REF] Gilbert | The module M1QN3, INRIA Rep[END_REF]) or bounded (lbfgsb-3.0, see [28]). The control vector can contain riverfloodplain parameters (bathymetry, friction coefficient), boundary conditions (lateral fluxes) and/or the hydrological model parameters on all/each subcatchment(s). The part of the cost function containing the misfit to observations is given for example by j obs = 1 2 ‖ Z (c )-Z obs ‖ 2 , with Z obs is the observed free surface and Z (c) is the simulated free surface.

SimHydro 2021: Models for complex and global water issues -Practices and expectations 16-18 June 2021, Sophia Antipolis -Pujol, L., Garambois, P-A., Monnier, J., Finaud-Guyot, P., Larnier, K., Mosé, R. -Integrated hydraulic-hydrological assimilation chain : towards multisource data fusion from river network to headwaters

RESULTS

A series of simple cases were designed with two goals in mind: (i) evaluate the accuracy of the meshing strategies described in subsection 2.1 to reproduce 1D and "real-like" behaviors, (ii) try the VDA chain capabilities on the newly integrated hydrological module in order to assess ways to improve the parameter search.

Hydraulic cases: 1D_eff modeling

On the 1D-2D interfaces

To evaluate the handling of internal boundaries, we model a steady flow in a simple slopeless prismatic channel using (i) a full 2D reference model and (ii) a 1D-2D model with one 1D-to-2D and one 2D-to-1D internal boundaries. The upstream boundary is an inflow hydrograph, the downstream boundary is a rating curve. In Fig. 3, we see that internal boundaries of the 1D-2D models behave as expected in these simple settings: downstream controls are propagated upstream and discharge is propagated downstream (the expected water depth is enforced at the rating curve). The slight discrepancy between 1D-2D and 2D waterlines can be attributed to mesh resolution variations, which influence the enforcement of the downstream BC. The 1D-2D and 2D models have equal downstream mass fluxes. Discharge is split evenly at the symmetrical diffluence.

Academic cases: 1D to 1D_eff comparison

The following cases allow a comparison of true 1D models to 1D_eff models obtained using our meshing strategy. The 1D code solves the SWE using a semi-implicit finite difference Preissmann scheme (see [START_REF] Larnier | On the estimation of river discharges from altimetry[END_REF]). The 2D code solves an explicit finite volume scheme. Both codes are part of the DassFlow platform. The model comparison is carried out on 3 cases, modeling progressively more complex hydraulic controls. All cases are mild sloped, 10 km long straight channels with constant XS dimensions. Case (a) is a simple, prismatic, rectangular channel. Case (b) is a variation of (a) featuring a slope break at 5 km, the slope remains mild. Case (c) is a variation of (b), where the 1D model has a parabolic XS.

All 1D models have an homogeneous Manning friction of 0.05 m -1/ 3 s . 1D_eff models have friction discretized in 2 patches: from 0 to 5km and from 5 to 10km. Their friction coefficients are either 0.05 m -1/ 3 s or given using the equivalent friction formula (eq.( 3)). A steady state comparison of all models (Fig. 4) allows to discuss the misfit of the effective 1D model to a reference 1D model as well as the pertinence of a theoretical equivalent friction. In case (a), the equivalent friction waterline matches that of the 1D model closely. This good fit is due to the simple structure of the case, where the normal depth is imposed downstream and constitutes the sole hydraulic control of this reach at quasi equilibrium. Equivalent friction is an appropriate effective tool to match a 1D uniform waterline. In case (b), the 1D waterline is gradually varied and features M1 and M2 backwater curves over the bottom slope break. In the downstream reach, the M1 curve is not reproduced by the 1D_eff model, which feature a constant normal depth imposed at the downstream boundary. In the upstream reach , the shape of the M2 backwater curve is matched by the 1D_eff models. Equivalent friction is insufficient to model gradually varied flows in this case. Keep in mind that in real cases, such as the Garonne river presented in this paper, hydraulic controls length may be shorter than in this academic case. In case (c), the shape of the gradually varied waterline, triggered by friction variaiton in space, is matched by the 1D-like models, but with a constant shift. This poor water depth fit is expected, as a shift in bathymetry should be necessary to better represent the reference 1D flow area (parabolic here) and consequently the waterline (see Fig. 2). To reduce the misfit, the bathymetry in the 1D-like model is shifted by the global average misfit between the 1D and 1D-like WSE. For either friction cases, this allows to fit the 1D WSE closely (not shown), although it does not allow to match the 1D water depth. A brief unsteady flow analysis is carried out on case (a). Fig. 5 presents the evolution in water depths at x=2.5km and x = 7.5km. Both stations feature a low water depth misfit between 1D and 2D models, although it is greater at high flows. To reduce this discrepancy, a friction power law n=ah b is introduced and calibrated manually (orange line). The obtained fit on this simple case is encouraging and future work will aim at inferring a and b patterns using data assimilation in order to approximate real flows variabilities in time and space with our 1D effective model.

Real scale real case

In this subsection, the above methodology of effective friction and bathymetry is applied on a large scale case: the Garonne River. This allows to discuss effective modeling using real scale width and bathymetry variations, i.e. hydraulic controls. A 2D model of the river, based on a fine Lidar DEM with inclusion in situ XS bathymetry (see references in [START_REF] Garambois | Inference of effective river properties from remotely sensed observations of water surface[END_REF]), is used as the truth. It has an homogeneous Manning friction coefficient of 0.05 m -1/ 3 s . This model contains complexities in (i) its local geometric hydraulic controls over a varied river bottom and (ii) the propagation of hydraulic signatures upstream and downstream. The goal of our effective modeling is to try to match both the local controls and propagation dynamics, or at least to reach the best compromise of the two. Three 1D-like models are built. A 1D-like model of the main channel is built using bankfull widths, with each quadrangular cell around 100m in length (see Fig. 6, left), an homogeneous friction of 0.05 m -1/ 3 s (equal to the friction of the reference 2D model) and with the 2D bottom bathymetry (minimum elevation of the fine 2D within the 1D like elements) (Model A). Equivalent friction, given by eq.( 3), is calculated from the reference model for each cell (Model B). Equivalent bathymetry is created by shifting the 2D bottom bathymetry by an amount equal to the global average difference in WSE between the 2D model and Model B. This gives model C. This simple approach conserves slopes in all 1D-like models. Recall that the effective model implicitly assumes a rectangular section by cells which tends to overstimate the real flow section and wetted perimeter (cf. Schematic XS view in Fig. 2).

These equivalent models are compared to the 'truth" for a bankfull steady flow (Fig. 6). This allows to focus on local controls representation. Model A features a global WSE misfit of 0.86m (bottom right, in blue). The misfit can stem from from differences in the mesh and in the evaluation and discretization of bathymetry. Model B, with spatialized equivalent friction, reduces the misfits at all points (average misfit of 0.56m), although local misfit variations are mostly conserved. Model C, built using results from Model B, reduces the misfit further, with an average global misfit of under 1mm. Recall that equivalent bathymetry allows to SimHydro 2021: Models for complex and global water issues -Practices and expectations 16-18 June 2021, Sophia Antipolis -Pujol, L., Garambois, P-A., Monnier, J., Finaud-Guyot, P., Larnier, K., Mosé, R. -Integrated hydraulic-hydrological assimilation chain : towards multisource data fusion from river network to headwaters fit the WSE, not the water depth. Moreover, using an effective flow vein representation may affected signal propagation speeds, which motivates the following analysis of the temporal variations of WSE at a station.

A simple propagation analysis is carried out at a station located in the middle of the reach (Fig. 7), with Model C. The discrepancy of relative WSE of the reference 2D model and Model C is constant over the time period. This is due to the simple global method used to determine equivalent bathymetry in Model C. A (further) shift in bathymetry allows to reduce this discrepancy significantly (not shown). No clear temporal misfit is observed during the propagation of the flood hydrograph in this real case.

The 1D_eff model with spatialized equivalent friction and bathymetry allows to model both local geometric hydraulic controls and signal propagation fairly. Future research will focus on the assimilation of spatiotemporal flow observables to improve these results, the determination of appropriate effective friction laws and regularization methodsdepending on inverse problems ill-posedness might be required.

Twin-experiment on hydraulic-hydrological model

In the following subsection, the behavior of the hydrological module integrated in the VDA chain is analyzed and a simple parameter inference is shown.

Direct GR4H run

Parameter values from [START_REF] Lay | MODELISATION HYDROLOGIQUE DANS UN CONTEXTE DE VARIABILITE HYDRO-CLIMATIQUE. Une approche comparative pour l'étude du cycle hydrologique à méso-échelle au Bénin[END_REF][START_REF] Santos | Continuous state-space representation of a bucket-type rainfall-runoff model: a case study with the GR4 model using state-space GR4 (version 1.0)[END_REF] were used to define a range of acceptable values for all parameters. Inputs (rain et evapotranspiration) from [START_REF] Santos | Continuous state-space representation of a bucket-type rainfall-runoff model: a case study with the GR4 model using state-space GR4 (version 1.0)[END_REF] were used. Initial model states (13 reservoir levels) were obtained from a long warm-up run. Fig. 8 shows the evolution of the inputs, outputs and reservoir states during a single rain SimHydro 2021: Models for complex and global water issues -Practices and expectations 16-18 June 2021, Sophia Antipolis -Pujol, L., Garambois, P-A., Monnier, J., Finaud-Guyot, P., Larnier, K., Mosé, R. -Integrated hydraulic-hydrological assimilation chain : towards multisource data fusion from river network to headwaters event. The water level of the production store changes slowly, while the other reservoirs are filled and emptied in a single event. This means that parameters ruling the production store are mostly responsible for low frequency output flow variations, while parameters ruling the other stores create higher frequency variations. Thus, the output flow sensitivity to parameters may be dependent on the observation window and the choice of the assimilation window may be key in determining the parameter search and inference results.

Hydrologic parameter inference

We present a sample hydrological parameter inference as proof of concept of the integrated assimilation chain. The inference is carried out over the sample event presented in the above subsection. The cost function contains the misfit of the hydrological module outputs such as j obs = 1

2 ‖ Q GR 4 (c )-Q GR 4 ,obs ‖ 2 .
Those simulated discharges correspond to one hydrological inflow on an open boundary (a tributary) of the hydraulic domain. The assimilation of the direct hydrological outputs is the first, most simple step towards more complex inferences. The inference is setup as a twin experiment, where a set of hydrological parameters is chosen as the "truth" and used to generate observations (of hydrological output). These observations are then used to infer the "true" parameters starting from a set of incorrect priors.

x 4 is the parameter that rules the 11 reservoirs of the Nash cascade. Given the short time scale of the considered event, we surmise that this parameter should be sensitive. All parameters except x 4 have their prior values equal to the "truth" and are constrained within bounds (Fig. 9, top, dotted colored lines) using a bounded minimization algorithm. The parameter x 4 is unconstrained. All parameters are sought simultaneously. An optimal parameter set is reached after 5 iterations. The cost has been reduced (Fig. 9, top, black line) and the inferred value of x 4 is closer to its "truth" value. x 1 has barely been modified and x 3 is at its upper bound for all iterations. x 2 has been modified but has not reached its bounds. In Fig. 9 (bottom), we can see the reduction of the initial misfit to the observations. The misfit has been decreased during the rain event, but decreased before the rain. Note that without bounds during the parameter search, inferred values are often outside of the accepted range. Once an iteration uses a value far outside of this range, later iterations will not fix the error. This investigation is carried out locally in parameter space as a mean to simply verify the assimilation chain. It shows the capability of VDA to the minimize the cost function in a simple case, while seeking the 4 hydrological parameters simultaneously. It also highlights the difficulty in optimizing hydrological models.

Tayloring optimization method adapted to global optimization over relatively large parameter spaces might be required as well as regularizations adapted to hydrologically coherent optimizations on river networks (e.g. [START_REF] De Lavenne | A regularization approach to improve the sequential calibration of a semidistributed hydrological model[END_REF]).

SimHydro 2021: Models for complex and global water issues -Practices and expectations 16-18 June 2021, Sophia Antipolis -Pujol, L., Garambois, P-A., Monnier, J., Finaud-Guyot, P., Larnier, K., Mosé, R. -Integrated hydraulic-hydrological assimilation chain : towards multisource data fusion from river network to headwaters

DISCUSSION

This paper presents a new integrated hydraulic-hydrological modeling chain based on the DassFlow platform (fine 2D hydraulics and variational data assimilation) and built in view to perform integrated catchmentfloodplains modeling and multi-sourced data assimilation. It is based on a new 1D-2D seamless hydraulic modeling method for representing river networks with 1D effective reaches connected to accurate 2D zooms and solved with the same finite volume scheme. First, the relevance of the implementations have been tested on many numerical cases and successfully compared to full 2D reference solutions, two tests mixing 1D effective and 2D zones are presented. Next, the capabilities of the 1D_eff modeling approach to reproduce water surface signatures is tested on academic and real cases. We show that spatialized friction coefficients/effective laws and effective bathymetry enable to fit reference flow signatures (free surface and propagations), both at the scale of local geometric controls and at the scale of the river reach. Then, we presents the coupled hydrological module GR4H state-space from [START_REF] Santos | Continuous state-space representation of a bucket-type rainfall-runoff model: a case study with the GR4 model using state-space GR4 (version 1.0)[END_REF] and a simple local inference of its parameters from hydrological outputs (boundary flux to the hydraulic domain), using a bounded minimization algorithm. This analysis of simple forward and inverse problems constitutes the first step towards inverse problem resolution on full catchments-river networks using the DassFlow integrated hydraulic-hydrological assimilation chain. Further work on 1D-2D modeling will focus on the fine-tuning of distributed hydraulic controls (channel and floodplain bathymetry-friction, source terms and boundary conditions) using the VDA algorithm and multisourced flow observations. This may require the use of spatialized effective friction laws and bathymetry as well as appropriate cost functions and regularization methods (see [START_REF] Larnier | On the estimation of river discharges from altimetry[END_REF][START_REF] Pujol | Estimation of Multiple Inflows and Effective Channel by Assimilation of Multi-satellite Hydraulic Signatures: Case of the Ungauged Anabranching Negro River[END_REF]). The integrated modeling approach and data assimilation will be applied to real cases, like the Adour case in [START_REF] Barthélémy | On operational flood forecasting system involving 1D/2D coupled hydraulic model and data assimilation[END_REF]. Finally, the link between hydraulic variables and hydrological parameters should be studied, using a setup to test the assimilation of multi-sourced observations of river network and information feedback to hydrological modeling components.

SimHydro 2021 :Figure 1 :

 20211 Models for complex and global water issues -Practices and expectations 16-18 June 2021, Sophia Antipolis -Pujol, L., Garambois, P-A., Monnier, J., Finaud-Guyot, P., Larnier, K., Mosé, R. -Integrated hydraulic-hydrological assimilation chain : towards multisource data fusion from river network to headwaters Design of a 1D-2D hydraulic-hydrological model of a river network. (a) River network conceptual representation with observables. (b) Representation of a river network using our 1D-2D methodology with integrated hydrological module. (c) Concept of a 1D-2D mesh at a confluence with a 2D "zoom" for overbank flows.

Figure 2 :

 2 SimHydro 2021: Models for complex and global water issues -Practices and expectations[START_REF] Miglio | Model coupling techniques for free-surface flow problems: Part I[END_REF][START_REF] Miglio | Model coupling techniques for free-surface flow problems: Part II[END_REF][START_REF] Monnier | Inverse algorithms for 2D shallow water equations in presence of wet dry fronts. Application to flood plain dynamics[END_REF] June 2021, Sophia Antipolis -Pujol, L., Garambois, P-A., Monnier, J., Finaud-Guyot, P., Larnier, K., Mosé, R. -Integrated hydraulic-hydrological assimilation chain : towards multisource data fusion from river network to headwaters Design of an equivalent 1D_eff channel with equivalent friction. Top: Schematic view and hydraulic misfit between a 1D cross section and the rectangular one imposed by our modeling choice in the effective 1D model using a single cell. Bottom: (a) Evolution of the hydraulic radius for 3 sample XS shapes (not represented) and (b) Evolution of the calculated equivalent Manning friction parameter in these XS according to eq.(3).

Figure 3 :

 3 Figure 3: Water depths at steady flow in academic slopeless models using a classical 2D mesh and a 1D-2D approach.

SimHydro 2021 :Figure 4 :

 20214 Models for complex and global water issues -Practices and expectations 16-18 June 2021, Sophia Antipolis -Pujol, L., Garambois, P-A., Monnier, J., Finaud-Guyot, P., Larnier, K., Mosé, R. -Integrated hydraulic-hydrological assimilation chain : towards multisource data fusion from river network to headwaters Effective friction analysis for steady waterline in academic cases. (a) Rectangular prismatic case. (b) Rectangular slope break. (c) Rectangular friction break. Equivalent friction values: (a): 0.0506, (b): (0.05060, 0.05040), (c): (0.05464, 0.05442).

SimHydro 2021 :Figure 5 :

 20215 Figure 5: Effective friction analysis for an unsteady waterline in case (a). The powerlaw is given by n=ah b , with a=0.0506 and b=0.05 . The injected hydrograph is symmetrical. Observations stations at x=2.5km and x = 7.5km.

SimHydro 2021 :Figure 6 :

 20216 Figure 6: Garonne meshes and steady state effective models analysis. Left: 2D and 1D_eff (Model B) meshes simulated water depth. Upper right: simulated WSE elevations in the 2D model, Model B and Model C. Lower right: misfit to 2D WSE for Model A (blue), Model B (equivalent friction, orange) and Model C (equivalent friction and bathymetry, green).

Figure 7 :

 7 Figure 7: Observation station during unsteady flow in the Garonne case near the gauging Verdun station. In red, the true water depth simulated by the 2D model. In blue, the Model C water depth + bathymetric shift (0.562m), .i.e the relative WSE simulated by the 1D_eff model.

Figure 8 :

 8 Figure 8: Hydrological model inputs, outputs and reservoir states evolution during a token event.

SimHydro 2021 :Figure 9 :

 20219 Figure 9: Parameter search and cost evolution. Dashed line are the parameter bounds for the minimization algorithm. The fine dashed line is the "true" value of the x 4 parameter.

http://www.math.univ-toulouse.fr/DassFlow

ACKNOWLEDGMENTS

This work is a part of the PhD thesis of LP. Research plan: LP, PAG, JM, PFG; Computational software DassFlow2D adapted from its previous versions by LP; Numerical investigations by LP with PAG and JM for analysis. JM is the principal designer of the inverse computational method. Fundings: PhD of LP is co-funded by CNES and ICUBE.


Rui Yan

Lanwei Wang

Michel Parrot

Xingguo Zhang

Zhe Hu

Comparison of ULF electric field recorded by ground-based stations in China and a low-altitude satellite

The correlation and synchronization between ground-based and space-based electric field in calm period, is the basis for a comprehensive study of earthquake monitoring. However, it is difficult to study their relevance because of large dissimilarities on data records. The DC-ULF electric field recorded by the DEMETER satellite and the geo-electric field observed in ground-based stations in China were comprehensively analyzed in this paper. The results show that the ionospheric electric field after removing the additional electric field caused by the satellite motion in the magnetic field has a stable annual variation trend; both the annual variation trend and the amplitude of ionospheric electric field above different locations in China are similar. The synchronous, corresponding and similar significant annual variation trends of both ionospheric electric field and geo-electric field in the same direction and with same frequency were found above some locations. The trend displays higher value in summer and lower value in winter. The Sq (Solar quiet) current system is the main source for these phenomena. It is noticed that the shapes and amplitude of geoelectric field for different stations are different, largely due to their individual underground layer conductivity, water level and so on. It is of great significance to explore the correlation and consistency between ground-based and space-based data in order to recognize the anomalies related to the seismic activity and to understand the LAI (Lithosphere-Atmosphere-Ionosphere) coupling mechanism.

Introduction

Electric field perturbations possibly connected with seismic activity is a matter of great concern. There are many reports about this study using the data from groundbased measurements [1][2][3][START_REF] Qian | Seismic electromagnetic observation techniques[END_REF][START_REF] Du | Specification for the construction of seismic station Geoelectrical station part 2:Geoelectrical field observatory (The Seismology Trades Standard of the[END_REF][START_REF] Ma | [END_REF] and data from space-based satellite experiment [7][8][9][10][11]. A joint analysis of variation related to seismic activity using these two kinds of data have also been studied in recent years [12][13][14][15][16]. They found that some electromagnetic anomalous phenomena appeared on the ground as well as in the ionosphere a e-mail: yanxiaoxiao best@163.com at a rough same time. They also pointed out that studying a number of earthquakes using ground-based and satellite data simultaneously can be a crucial and effective way of investigating the Ultra-Low-Frequency/Extremely Low Frequency (ULF/ELF) electromagnetic anomaly triggered by the seismic activity.

These anomalous phenomena are considered as the manifestation for the existence of a lithosphere-atmosphere-ionosphere (LAI) coupling [17,18]. As for this LAI mechanism, two main models are currently considered. One is the Electro Motive Force (EMF) model which is based on the assumption that ionization of lower atmosphere (mainly due to radon) leads to the seismic-related electric field formation in the ionosphere [START_REF] Pulinets | Lithosphere -Atmosphere -Ionosphere -Magnetosphere -Coupling (LAIMC) -a Concept for Pre-earthquake Signals Generation[END_REF]. Another model is related to air ionization at rock surfaces, which is a massive air ionization prior to major earthquakes which increases the electrical conductivity in the air column and may cause ionospheric perturbations because there are conditions under which the air ionization at the ground-to-air interface produces exclusively positive airborne ions. If the atmosphere becomes laden with positive ions, it will introduce an electrostatic (Coulomb) perturbation that can be expected to affect the vertical E field which is a part of the global electric circuit. Then, it will affect the ionosphere [START_REF] Freund | [END_REF]. So there are a lot of works at this time to study the anomalous phenomena related to earthquakes and to try to understand their generation mechanism

In this work, it is very meaningful to study the electric field characteristics during calm period synchronously. The study of space-and ground-based electromagnetic observational data and their correlation can help researchers to understand the LAI coupling. Especially the basic characteristics of space-and ground-based data and their correlation during the period without events such as earthquakes, storms etc. can help us to understand and recognize the characteristics of these data during the period with events. It is also of great scientific significance for comprehensively utilizing the satellite and ground electromagnetic monitoring data for earthquake monitoring research [21]. However, it is difficult to study their relevance because of large dissimilarities on data records. Very few papers have been published on this topic.

Only some studies on the characteristics of these two kinds of data during calm period separately can be found. About the geo-electric field characteristics, Mogi et al. [22] have studied the short-term and the long-time variation of electric field. Shortterm electric field variations were found to correspond mainly to geomagnetic activity, while long-term variation was mostly gradual shift and was clearly uncorrelated to precipitation and to ground water level variations. Ye et al. [START_REF] Ye | Master thesis Lanzhou Institute of Seismology[END_REF] studied the spectrum characteristics of geo-electric diurnal variation based on data from stations in China. This study showed that in geo-electric diurnal variation, the amplitude of the 12 h semidiurnal wave is the largest. Cui et al. [START_REF] Cui | [END_REF] has studied the time and frequency characteristics of the geo-electric field diurnal variation along two longitude chains and two latitude chains and found the same conclusion with Ye [START_REF] Ye | Master thesis Lanzhou Institute of Seismology[END_REF]. Concerning the characteristics of the ionospheric electric field during calm period, Yu et al. [25][26][27] studied the yearly variation of ionospheric electric fields at mid-and low-latitudes under the condition of low solar and quiet geomagnetic activities with a theoretical model of ionospheric electric fields. Results revealed that the coupled electric field would show a pronounced semi-annual variation due to the "coupling effect" along the magnetic field lines between two hemispheres, whereas the uncoupled electric field have an annual variation in respective hemisphere.

The study in this paper closely follows a study by the author in 2013 [21], where research on the characteristics and relationship using geo-electric field observed in two stations and corresponding ionospheric electric field measured by DEMETER were done. Their results showed that during the night, the annual trends of electric field from space and ground-based observational data have exhibited correlation and synchronization every year. In order to research whether or not there is a certain relationship between the ionospheric electric field and geo-electric field in a wider area; we extend the data set to all over China. We also update the data processing method and selected the north-south (NS) direction both in geo-electric field and ionospheric electric field as the main direction in the research. Then we evaluate the correlation between them. Data source and selection is presented in Section 2. We will introduce data processing method in Section 3. Section 4 presents results. They are discussed in Section 5 and summarized in Section 6.

2 Data source and selection

Data source

The ground-based electric field data are recorded by the instrument ZD9A in ground stations. ZD9A is an electric field detector which has been developed in China and installed in many stations [28,29]. The frequency band is 0-0.1 Hz and the sampling rate is 1 Hz. Horizontal component, north-south (NS) and west-east (WE) directions of geo-electric field have been registered, with a set of orthogonal configurations. Each of them is composed of two non-polarized solid state electrodes with spacing of 300 m or 400 m, and buried 3-4 m below the ground surface. At present, geo-electric field observations related to earthquakes in China are mostly fixed-point measurements, with 85 stations throughout the country [START_REF] Du | Specification for the construction of seismic station Geoelectrical station part 2:Geoelectrical field observatory (The Seismology Trades Standard of the[END_REF]16,[START_REF] Ye | Master thesis Lanzhou Institute of Seismology[END_REF][START_REF] Cui | [END_REF].

The space-based electric field data from the French spacecraft DEMETER have been used. DEMETER was a microsatellite launched in June 2004 and stopped in 2010. Its orbit was circular and quasi sun-synchronous (10.30 and 22.30 LT) with an altitude of about 700 km and an inclination of 98 • . The altitude was changed to about 660 km in December, 2005 [30]. We have focused on the analysis of electric field data measured by the ICE instrument. There are four stacer booms used by the ICE experiment, which have a length of 4 m. The probes E1 and E2 are deployed parallel to the satellite Y -axis, i.e. perpendicular to the satellite XZ plane that is the orbital plane for the nominal attitude of the satellite. The probes E3 and E4 are deployed at 45 • from the X-axis, which is the nadir direction, and at distances of about 15 cm on both sides of the satellite XZ plane. All 3 components of electric field in satellite coordinate are Ex (upward), Ey (vertical to satellite orbital plane), and Ez (opposite to the satellite direction) [31,32]. In order to keep the same frequency range for the two sets of electric field data, we select the DEMETER DC/ULF (0-15 Hz) data. The signals are sampled at 39.0625 Hz, and the time resolution is 6.55360 s. The data unit is in mV/m [31,32]. We can get the data of the E1, E2, E3, and E4 probes, and the electric field Ex, Ey, and Ez in the satellite coordinate system.

Data selection

The stations in China are not evenly distributed. There are fewer stations in the west due to geological conditions. In order to study the geo-electric field in China, it is better to obtain data in the stations which have fixed distance intervals and cover the whole country. On this basis, considering the data quality and the distribution of stations available all over China, geo-electric field data observed from 30 stations have been selected among the total 85 stations. The latest digital ZD9A instruments have been deployed at most of stations in China since the second half of 2007. We have firstly analysed the annual variation trend of the observed data. It shows that only data from 5 stations display obvious good regular annual variation trends. For 9 stations, after removing the interference data [START_REF] Qian | Seismic electromagnetic observation techniques[END_REF], a rough annual variation trends can be seen. So for these 30 stations selected there are annual variation trends in 14 stations whereas for other most stations, there are no regular variation trends. Considering the observation periods of these two kinds of data we obtained data of two complete years (from 2008 to 2009).

All ionospheric electric field data in the whole territory of China have similar obvious annual change trend, the correlation between the two kinds of observation data is completely dependent on whether the geo-electric field data has annual change trend or not. If there is no annual trend in the geo-electric field, obviously, it is not meaningful to compute the correlation coefficients of these two kinds of data. It is better to search correlation of data with good and rough annual trend in 14 stations. But for 4 stations. (Guazhou, Songshan, Jiayuguan, Chengdu), the data sets are not complete for 2008 and 2009. So finally, we have selected geo-electric field data in the remaining 10 stations (yellow squares in Fig. 1) with both, good and roughly annual variation trend, and corresponding ionospheric electric field data, to carry out correlation research.

There is no obvious large variation of the electric field in the ionosphere inside a large area, and we will have the same revisited orbits for several nearby stations (revisited orbits are the orbits passing by the same sub-satellite point, which means these revisited orbits are roughly above the same location). So we need not to process ionospheric electric field corresponding to all 30 geo-electric field stations. According to the different annual trends of the geo-electric field in these 30 stations, and considering the station positions for a better distribution of observation stations over China, the night time ionospheric electric field corresponding to a selection of 18 stations (green triangles in Fig. 1) among the 30 stations was processed and analysed in order to compare their different annual variation trends. These 18 stations are corresponding to 5 stations with good annual trend, 6 stations with rough annual trend and 7 stations with no annual trend, including the positions corresponding to the above 10 selected stations. From these data above 18 locations, we can get the characteristics of ionospheric electric field all over China. During the DEMETER operation period from 2004 to 2010, only data observed in 2006-2009 were selected, because the altitude of DEMETER was adjusted from 710 km to 660 km in mid-December 2005 [33] and the sun spot in 2009 is lowest in the solar cycle whereas the one in 2010 begin to ascend. So for each station, all revisited orbits during 2006-2009 are used, similar to the selection described by Yan et al. [21].

We have focused on all data during nighttime with Kp index less than 3 in order to avoid perturbation by the solar activity.

3 Data processing

The processing of electric field in ionosphere

In order to comprehensively compare the space-based and ground-based electric fields, we should try our best to get similar parameters for these two kinds of observations. For the ionospheric electric field registered from DEMETER, there are 4 processing steps as follows.

(a) The data are transformed from the satellite coordinate system to the geographic coordinate system firstly. The data in the Ez direction which is identical to the NS direction of the geo-electric field are selected. Then we designed a low-pass filtering to extract the frequency band of 0-0.1 Hz from the DC-ULF (0-15 Hz) ionospheric electric field so that the space-and ground-based electric fields is consistent as much as possible in direction and in frequency bands.

(b) We need to calculate the additional electric field caused by the motion of the satellite across the Earth's geomagnetic field [31,[START_REF] Pisa | Pre-processing of the ULF Waveform Fluctuations Above Seismic Active Regions Observed by DEMETER, f-2[END_REF]. The additional electric field is in a geographic coordinate system, so data in Ez direction (corresponding to NS direction) is also selected. The time resolution of additional electric field is 30 seconds and then it is interpolated at the resolution of the original data.

(c) The temporal sequences of ionospheric electric field corresponding to each station are obtained. Firstly, the data point A which is nearest to the station position in every corresponding revisited orbit is found. Then we select data within 100 km range around point A and calculate the average value of them. The reasons for selecting 100 km are as follows: DC/ULF electric field signals are sampled at 39.0625 Hz, and there are 256 points per 6.5536 s. We have selected two sets of DC/ULF electric field data of about 13 s, and, as the satellite velocity is 7.8 km/s, this corresponds to a distance of 100 km where we assume that there is no obvious variation for electric field in the ionosphere. The averaged data can remove the small variation and, then can reflect changes above that area.

(d) The average values of all successive revisited orbits are then obtained. Because there is only one revisited orbit every 16 days (2006∼Jun, 2008) or 13 days (July, 2008∼2009), and the revisited orbits were gradually changing from 2008, so the time resolution is not fixed. Only several tens of average values can be acquired during one year, and there are about two data points per month (Fig. 2).

(e) In order to decrease the influence of long-term variation of solar activity, we have eliminated the linear trend of temporal sequences using a least-squares fit (Fig. 3).

(f) Then we got the fitting line of these time series for each year using a Fourier approximation (Fig. 4).

The processing of geo-electric field

Interference data refers to abnormal changes which are caused by the observation process itself, which includes non-natural changes such as single point jump, step change, and periodic change. By removing them, the obvious unreasonable observation data can be distinguished and eliminated. Then the observation data that can truly reflect the actual situation of each station can be obtained [START_REF] Qian | Seismic electromagnetic observation techniques[END_REF]. So the interference data should be removed from the geo-electric field observation data firstly. There are different interference cases for different stations, so no uniform standard can be used to remove interference data for different stations. For different interferences, we adopt different methods. We remove the single jump point directly; for step change, after checking these interference reasons, we first remove them directly, and according to the data before and after them, the data in this range are fitted and then added into the data set; and in case of periodic change, a filter is designed to remove them.

A night time (20:00-24:00) temporal sequence of daily average in each station should be generated, which correspond to the selected revisited orbits time of DEMETER.

Finally, the geo-electric field data will drift unsteadily because of unsteady of electrodes [START_REF] Ma | [END_REF], similar to ionospheric electric field data, the linear trend caused by electrodes drift is also removed from the time series using a least-squares fit.

Correlation calculation

For the correlation between geo-electric field and ionospheric electric field, we focused on the geo-electric field data with obvious annual variation trends in the 10 stations (purple triangles in Fig. 1) and the corresponding ionospheric electric field. The common observation periods of the two kinds of data are 2008 and 2009. The scatter plot and a linear regression model "f (x) = a + bx" are used to find whether there is correlation (Fig. 6). The correlation coefficients (Tab. 1) are calculated by the equation (1).

R(X, Y ) =

Cov(X, Y )

V ar [X] V ar [Y ] (1)
In equation ( 1), for each station, X is the temporal sequence of the geo-electric field in 2008 and 2009, Y is the temporal sequence of the corresponding ionospheric electric field, and the length of X and Y temporal sequences are identical. But for different stations, the length is a little bit different.

Results

Figure 2 shows the temporal sequences of data in geographic NS direction of the ionospheric electric field corresponding to 18 stations selected during 2006-2009, which are the results of the processing explained in Section 3.1d. The points of different colours are related to the 18 stations (green triangles in Fig. 1). As seen in Figure 2, there is similar annual variation trend of ionospheric electric field corresponding to the 18 stations. During one year, the amplitude of ionospheric electric field is higher in summer and lower in winter. For the different years, the annual variation trend is stable and similar, but the amplitude which is between -30 mV/m and 20 mV/m is not the same. The amplitude in 2009 is lower than in the other years. Figure 3 also shows the same temporal sequences as in Figure 2, but the linear trend has been removed. Comparing to Figure 2, characteristics are consistent, however, the amplitude difference between 2009 and other years is weakened, and the amplitude range of each year is reduced. Figure 4 shows the temporal evolution characteristic of the ionospheric electric field year by year which is underlined by a fitting line. The temporal sequences in Figure 4 are similar to the ones in Figure 3, but are displayed in different years. It is obvious that the annual variation characteristic is similar to sinusoidal (Fig. 4). We can also get the evolution model of temporal sequences in one year. The basic model of temporal evolution in each year is given by the equation (2). The coefficients (a, b, n, and c) will change with observations data in different years.

f (X) = a + b * sin(n * X + c) (2)
In order to avoid a possible effect from electrode floating potential [22,36], we study the correlation between ionospheric electric field and geo-electric field using the temporal sequences after removing the linear trend. The temporal sequences of the two kinds of electric fields in NS direction during 2008-2009 are shown in Figure 5 using double coordinates because of the huge difference of magnitude between them. The Y -axis in left side is ionospheric electric field (unit, mV/m) while the Y -axis in the right side is geo-electric field (unit, mV/km). Figure 5a∼j are corresponding to different stations. It can be seen directly that the two kinds of data have very similar annual variation characteristics, which looks like a half sinusoidal curve. Figure 6 indicates the results of scatter plot and regression analysis. Similar to Figures 5 and6a∼j are also corresponding to different stations. The red line is the regression line in each case. The "R" is the correlation coefficient. The function "f (x) = a + bx" is the regression model for each data set. In Figure 6, an obvious linear correlation between two kinds of data is found. For data at different positions, the regression models are nearly similar. Table 1 gives the list of the correlation coefficients of the corresponding data from space-and ground-based observation in the 10 positions selected in this study. We can find that after removing the trends from ground-based and space-based data, the correlation coefficients of data corresponding to the 5 stations with good, annual variation are good, which in 4 stations are more than 0.7. For the 5 stations with annual variation most of the correlation coefficients are more than 0.5 (Tab. 1).

Discussions

From these analyses, the data of ionospheric electric field in different locations all over China showed the nearly same annual variation trends (Figs. 2 and3). The amplitudes in different locations and during different years are also very similar. It should be noted that the amplitude of data in 2009 is slightly lower than data in the other three years. In general, for the ionospheric electric field in China, we can get the fixed annual variation model of approximate a half sinusoidal for each year (Fig. 4), although the coefficients are different in different years. This type of annual variation trend (red lines in Fig. 4) is consistent with the results of the annual variation characteristics of the ionospheric electric field studied by Yu et al. [25][26][27].

According to the statistics results, there are different trends of geo-electric field every year in many stations. Only in some stations, there are similar or same annual variation characteristics with ionospheric electric field. The amplitudes of data both with and without annual variation trends are also different depending on cases. Zhang et al. [16] thinks that the differences in geo-electric field amplitude and shapes between different stations and among ground observing data reflect the different response process at stations with different underground electromagnetic structures and the response mechanism to current system in lithosphere. Li et al. [36] studied the geo-electric field observed at three stations in the Shanghai area and found that data in even very close three stations can vary greatly. They pointed that the huge differences of curve shape, amplitude and phase are caused by different electrical properties under the stations. The different laying way of the electrode, the different distance between the electrodes are other reasons, which is also mentioned by Mogi

(g) (h) (i) (j)
Fig. 6. continued.

et al. [22]. Tan et al. [37] think rock structure, the quaternary sediments, fracture and the bedrock, water content and permeability of rocks, tectonic activity, and other factors also affect the waveform of spontaneous electric field. So the amplitude and shape of data in different stations may be different due to the local Earth surface, the underground layer conductivity, and also the electrode contamination.

For the geo-electric field with annual variation trend in some stations, similar or nearly same annual variation characteristics of electric field both from space-based and ground-based measurement can be found (Fig. 5). A remarkable result is that although the geo-electric field variations are slightly different for these ground-based stations, the ionospheric electric fields measured by DEMETER above these stations well follow same variations. We also calculated the correlation coefficients of the two kinds of electric field and found that the correlation coefficients are good (Fig. 6 and Tab. 1). If we only consider the annual variation trend, there is a good correlation between the geo-electric field and ionospheric electric field for in these stations. This good correlation cannot be explained by the mechanism of LAI coupling. The study in this paper concerns electric field under normal conditions but in case of seismic activity, it could explain why the electric field induced in the lithosphere-atmosphere could propagate up to the ionosphere and trigger ionospheric perturbations [7,[38][39][40][41][42].

In this section, we try to explain why there is a synchronous, corresponding and similar significant annual variation trend between electric field observed both on ground and in space. Sorokin and Hayakawa [17] find that there is no static electric field of lithospheric origin in the ionosphere. So the annual variation trend of these two kinds of electric field is not caused by source from lithosphere. It is probably from outside the Earth. The generator process theory holds that: an induced electric field is generated by a neutral wind dragging plasma and cutting magnetic field; at the same time, due to the spatio-temporal unevenness of ionospheric conductivity, charge accumulation occurs, thus forming a polarized electric field, which produces observable plasma drift in the ionospheric F region. This is the main component of the ionospheric electric field at low and medium latitudes. In the high latitudes, besides electric field caused by generator process, electric fields from the magnetosphere account for a considerable proportion of the electric field. The generator process theory cannot explain ionospheric electric field completely in the high latitudes. However, in the middle and low latitudes far away from the polar region, the influence from the magnetic layer is generally very small, so the generator theory can be applied well [START_REF] Xiong | The Ionosphere Physics Probability[END_REF]. So the widely accepted generator process is a major source of ionospheric electric field in the middle-low latitudes and the equatorial regions. This process has been described in detail in many previous works [START_REF] Chapman | Geomagnetism[END_REF][START_REF] Matsushita | Physics of geomagnetic phenomena 2 volumes (International Geophysics Series)[END_REF][START_REF] Richmond | [END_REF]. The electric field and current systems mainly include the Sq (Solar quiet) current system which is the current system that generates the sun day variation of the quiet solar geomagnetic field and the equatorial current in the middle-low latitudes [25,[START_REF] Kamid | Magnetosphere-ionosphere coupling[END_REF]. Through modelling the Sq current system at solar maximum and minimum, Takeda et al. [START_REF] Takeda | [END_REF] found that the ionospheric conductivity also exhibits a seasonal variation with larger conductivities generally occurring in the summer hemisphere. A comprehensive review of early studies of the Sq current system is provided by Chapman and Bartels [START_REF] Chapman | Geomagnetism[END_REF]. The dominant characteristics of the Sq current system are its seasonal and solar cycle variations. Owing to differences in conductivity and winds, the Sq current system is not hemispherically symmetric and this feature exhibits a marked seasonal variation with the current intensity in the summer hemisphere being greater than the current intensity in the winter hemisphere. Walker et al. [10] found that a predominantly annual variation of Sq current is observed in the Northern Hemisphere with the peak value occurring near the June solstice. All these research results are consistent with the result of the ionospheric electric field in this paper. So the annual trend of ionospheric electric field is caused by Sq current system, which has the same seasonal variation that is high value in summer and low value in winter. Qian and Lin [START_REF] Qian | Seismic electromagnetic observation techniques[END_REF] pointed that the geo-electric field observed by ZD9A in fact has two parts, namely the telluric electric field and spontaneous electric field, based on the different field sources. The telluric electric field is the changing electric field generated by various field sources located outside the Earth and is distributed throughout the entire surface or larger regions. It is associated with the changing magnetic field in the geomagnetic field and has global or large regional characteristics. The spontaneous electric field is the relatively stable electric field formed on the Earth's surface by the subsurface media caused by various physical and chemical actions. It is made of the electric field by distribution of orebody, groundwater and various water systems, which generally has large changes in the horizontal and vertical gradient, in homogeneity, with strong local and even the local characteristics. Now it is difficult to completely separate the two parts. In general, the diurnal variation and its frequency represent the telluric electric field part in the data, which is homogenous with the ionospheric electric field and is correlated with the daily geomagnetic fluctuation [START_REF] Qian | Seismic electromagnetic observation techniques[END_REF]22,[START_REF] Ye | Master thesis Lanzhou Institute of Seismology[END_REF]. As early as 1960, the long-term seasonal and solar cycle variability and significant day-to-day variability in the Sq current system are found [49]. So the annual variation of geo-electric field also reflects the part of telluric electric field, which has both annual and diurnal characteristics. So the annual variation characteristic of geo-electric field is also caused by Sq current system probably.

There is another phenomenon that needs to be discussed. The amplitude of data differs largely between ground-and space-based observation. In this paper, the amplitude of ionospheric electric field concentrated in several tens mV/m, whereas the amplitude of geo-electric field is from a few mV/km to hundreds of mV/km. Vartosos et al. have reported their observations which the SES (seismic electric signals) occurred in certain regions with amplitude of up to 250 mV/km by 7-260 h prior to earthquakes using the VAN method [1][2][3]50]. Chmyrev et al. [7] analyzed the vertical component of quasi-static (DC) electric field Ez, based on the data recorded by INTERCOSMOS-BULGARIA-1300, and reported an anomalous increase of 3-7 mV/m in the quasi-static electric field. Zhang et al. [16] found that ULF electric field showed perturbations a few days before Wenchuan earthquake with amplitude of about 3-5 mV/m, and the amplitude of ground electric field anomalies was from 3 mV/km to 100 mV/km. Xiong et al. [51] and Sorokin et al. [52] pointed out that the anomalous signals of ground observation of electric field, were often in the range of a few tens mV/km to a few hundred mV/km, which means that the amplitude of ground electromagnetic emissions was still far lower than that of space ones, even after consideration of the amplification in the atmosphere of DC/ULF electric field. Yu et al. [25][26][27] simulated ionospheric electric field which varies from less than 1 mV/m to several mV/m. Ye et al. [START_REF] Ye | Master thesis Lanzhou Institute of Seismology[END_REF] and Li et al. [36] researched the geo-electric field observed at stations in China and found that data in even very close stations can vary from several mv/km to more than 100 mv/km. It can be seen from these previous results that the amplitude of the geo-electric and ionospheric electric fields in this paper is in accordance with the observations performed with and without seismic activity. The electrodynamic model of the atmosphere-ionosphere coupling described in Sorokin et al. [40,41] and Sorokin and Chmyrev [START_REF] Sorokin | Atmosphere-ionosphere electrodynamic coupling[END_REF] also cannot be used to explain the different magnitudes of the electric field, because the electric field amplitude is attenuated largely during the propagation from the Earth's surface up to the ionosphere (around several tens of V/km to tens of mV/m), which contradicts with the actual observation data in this paper. According to coupling mechanisms, the anomalous electric fields associated with earthquakes are mostly due to the abnormal electric field in the atmosphere above the Earth's surface. Such electric field is relatively large, especially during thunderstorm activity and other events occurred [START_REF] Sorokin | Atmosphere-ionosphere electrodynamic coupling[END_REF][START_REF] Park | [END_REF][55]. In this paper, the geo-electric field used is the electric field generated by an external field through the underground medium and the electric field induced by the medium itself. Then the amplitude of this electric field is of the order of a few hundred mV/km, which is different from the electric field intensity evaluated in electrodynamic model. In addition, it must be noted that in their theory of peroxy defects, Freund et al. [56] and Scoville et al. [57] claimed that stress activation of positive holes, a prerequisite for pre-earthquake phenomena, is efficient for earthquake hypocentres located in the upper to middle crust. This is the location of the geo-electric field studied here. Shen et al. [58] proposed a new mechanism where the magnetite aggregate in the crust could be a converter from thermal into electric during the process of earthquake preparation, namely, abnormal geo-electricity is caused by a thermoelectric mechanism. Based on the thermoelectric coefficient value tested in the experiments and on the thermoelectric and geo-electric theories the thermoelectric field was calculated. It was showed that the field intensity can reach 10 -6 V/m ∼ 10 -4 V/m. This model result is very consistent with the geo-electric anomalies detected by many earthquake monitoring. At the end, we have not explained the reason of the different amplitudes between ground-and space-based data observed in this paper, but this study shows that these amplitudes are reliable. Further research is needed to investigate this point.

Conclusions

For the ionospheric electric field, at 660 km altitude during night all over China, we have confirmed that there is obvious annual variation trend. We can also get the variation model with approximate sinusoidal.

For geo-electric field in some positions, there are similar or same annual variation trends with ionospheric electric field.

The synchronous, corresponding and similar significant annual variation trends of both ionospheric electric field and geo-electric field are identical at some locations, but the definite correlation between ionospheric electric field and geo-electric field everywhere is not true. Further researches on these phenomena are necessary and meaningful.

This research will hopefully help us to understand the nature of electric field in the ionosphere and on ground, and its penetration mechanism from ground into the ionosphere in case of catastrophic events occurring at the Earth's surface (earthquakes, tsunamis, hurricanes,...).

In this paper, we cannot measure the ionospheric electric field at a fix position and with uniform time intervals, so the results are a little bit biased. But as for the data available at the moment, the results in this paper are objective and effective. When DEMETER was launched in 2006, Zlotnicki et al. [59] had tried to monitor active faults and volcanoes using ground-based electromagnetic data combined with remote sensing based on DEMETER mission, but they did not get ideal results at that time. There are still a lot of work needed to be done, not only in LAI coupling, but also concerning the characteristic in normal conditions. A large-scale normalized electromagnetic observational network for detecting seismic precursors has been built in China. The CSES (China Seismo-Electromagnetic Satellite) had been launched in February 2018, which is proposed to be the first experimental satellite for China earthquake monitoring system from space [60,61]. More studies in this research field are meaningful.

A part of this work was done during the stay of Rui Yan in LPC2E, who expresses sincere thanks to Professor Jean-louis Pinçon for useful and fruitful discussions. The authors also express sincere thanks to Professor Q. Z. Ma for the advices of geo-electric field data processing and D. Pisa for advices about additional electric field calculation. This work is supported by the National Natural Science Foundation of China (41404058) and special fund for public welfare of Scientific Research Project of ICD, CEA (No. ZDJ2016-19). This work is mainly related to data recorded by the ICE of the microsatellite DEMETER that was operated by the French Centre National d'Etudes Spatiales (CNES) and the data measured by ZD9A which was operated by China Earthquake Networks Centre China Earthquake Administration (CEA).

Fig. 1 .

 1 Fig. 1. Distribution map of geo-electric field stations in China.

Fig. 2 .

 2 Fig. 2. Temporal sequences of the ionospheric electric field in NS direction of geographic coordinate system during 2006-2009. The blue vertical lines separate the different years. Data corresponding to different stations are plotted with different colors. The X-axis represents the time (year) whereas the Y -axis shows the amplitude of the ionospheric electric field (mV/m).

Fig. 3 .

 3 Fig. 3. Same temporal sequences as in Figure 2, but after removing the linear trend.

Fig. 4 .

 4 Fig. 4. The same temporal sequences as in Figure 2, but displayed in different years. The red smoothed line is the fitting line of temporal sequences. (a) data in 2006, (b) data in 2007, (c) data in 2008, (d) data in 2009.

Fig. 5 .

 5 Fig. 5. The temporal sequences of electric field in NS direction obtained from space-and ground-based observation during 2008-2009 after removing the linear trend for different stations. (a) Shandan, (b) Wushi, (c) Qianling, (d) Datong, (e) Guyuan, (f) Tengchong.

Fig. 5 .

 5 Fig. 5. continued.

Fig. 6 .

 6 Fig.6. The scatter plot of the same temporal sequences as in Figure5. The red line is the linear regression line. "R" is the correlation coefficient (see text). "F (x) = a + bx" is the regression model for each data set (see text). (a) Shandan, (b) Wushi, (c) Qianling, (d) Datong, (e) Guyuan, (f) Tengchong, (g) Xiaxian, (h) Rushan, (i) Manzhouli, (j) Ganzi.

Table 1 .

 1 The List of the correlation coefficients of the corresponding data from space-and ground-based observations.

		Positions	Correlation coefficient Slope
	1	Shandan	0.79	0.55
	2	Wushi	0.77	0.07
	3	Qianling	0.77	2.1
	4	Datong	0.76	0.98
	5	Guyuan	0.65	0.41
	6	Tengchong 0.64	0.93
	7	Xiaxian	0.60	0.29
	8	Rushan	0.53	0.19
	9	Manzhouli 0.51	1.25
	10 Ganzi	0.41	0.03

Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.