Abderrahim Jardani
email: abderrahim.jardani@univ-rouen.fr

T M Vu

Pierre Fischer

Use of convolutional neural networks with encoder-decoder structure for predicting the inverse operator in hydraulic tomography

Introduction

Since the 1980s, the hydraulic tomography has been adopted in hydrology as an effective technique for mapping the heterogeneity of the hydraulic properties of aquifers [START_REF] Neuman | Stochastic continuum representation of fractured rock permeability as an alternative to the R.E.V. and fracture network concepts[END_REF][START_REF] Gottlieb | Identification of the permeability distribution in soil by hydraulic tomography[END_REF][START_REF] Bohling | Steady shape analysis of tomographic pumping tests for characterization of aquifer heterogeneities[END_REF][START_REF] Yeh | Hydraulic tomography: Development of a new aquifer test method[END_REF]. Indeed, this approach can effectively provide insights into the spatial variability of the hydraulic transmissivity and storage coefficient of a porous or fractured aquifer through a combined interpretation of hydraulic data obtained from several pumping tests [START_REF] Bohling | Steady shape analysis of tomographic pumping tests for characterization of aquifer heterogeneities[END_REF][START_REF] Berg | Three-dimensional transient hydraulic tomography in a highly heterogeneous glaciofluvial aquifer-aquitard system[END_REF][START_REF] Cardiff | Hydraulic conductivity imaging from 3-D transient hydraulic tomography at several pumping/observation densities[END_REF][START_REF] Fischer | Hydraulic tomography of discrete networks of conduits and fractures in a karstic aquifer by using a deterministic inversion algorithm[END_REF]. The approach consists of using an inversion algorithm, either deterministic or stochastic, to retrieve the best models that could match the observed piezometric data (Fu & Gómez-Hernández, 2009a;[START_REF] Kitanidis | Introduction to geostatistics: applications in hydrogeology[END_REF]. In general, the inverse solution is not unique for that reason the inverse formulation usually incorporates a regularization term in the optimization of the objective function ( ) m Ψ that is defined as the sum of data misfit and regularization terms, as below [START_REF] Tarantola | Generalized nonlinear inverse problems solved using the least squares criterion[END_REF]:

[ ] [ ] T T 1 1 obs d obs m 0 0 Data misfit term constraint term ( ) F( ) F( ) , - - - - - -     = +     m m m m m m m
1444444444442 444444444443 144444444 42 4444444443

Ψ h C h C (1)
where the data misfit term assesses the pertinence of the model in term of matching the data and is expressed as a sum of square differences between the observed hydraulic data obs h and the numerical hydraulic data derived by solving a forward operator F( ) m . In hydraulic tomography, the forward problem involves a numerical discretization for solving the groundwater flow equation subject to Darcy's law in transient or steady modes.

[ ] log = -T m is the spatial distribution of the negative logarithm of the hydraulic transmissivity T , which is regarded as an unknown field and is expressed logarithmically to ensure its positivity during the inverse process. d C is a diagonal matrix for incorporating data uncertainty in the optimization. On the other side, the regularization term is expressed with a prior model 0 m and its covariance matrix m C . This term can be obtained from geological, geophysical and tracer investigations to constrain the inversion to offer only plausible solutions [START_REF] Lochbühler | Structure-coupled joint inversion of geophysical and hydrological data[END_REF][START_REF] Zhao | On the importance of geological data for hydraulic tomography analysis: laboratory sandbox study[END_REF][START_REF] Soueid | Specific storage and hydraulic conductivity tomography through the joint inversion of hydraulic heads and self-potential data[END_REF]. In general, the use of geostatistical constraints is still widely applied in hydraulic tomography by using statistical properties such as mean and covariance to obtain an imagery of hydraulic properties with smooth features [START_REF] Kitanidis | Introduction to geostatistics: applications in hydrogeology[END_REF][START_REF] Yeh | Hydraulic tomography: Development of a new aquifer test method[END_REF].

Once the parameters of the objective function are formulated, its optimization is processed iteratively using one of the three categories of optimization algorithms: deterministic [START_REF] Li | A Geostatistical inverse modeling of transient pumping tests using temporal moments of drawdown[END_REF], stochastic [START_REF] Jimenez | Smart pilot points using reversiblejump Markov-chain Monte Carlo[END_REF] or global [START_REF] Castagna | A Bayesian approach for inversion of hydraulic tomographic data[END_REF]. These algorithms involve a repetitive numerical solving of the forward problem by applying finite elements or finite difference methods to assess the ability of the proposed model to reach the convergence state. Therefore, the computation time is dependent on the numerical tools, the dimensionality of the unknown parameters and the nature of the optimization algorithm employed. In general, deterministic algorithms are the most used in high-dimensional inverse problems due to their advantage of reaching convergence in a short time compared to stochastic or global algorithms by exploiting the geometric properties of the objective function (gradient and Hessian) [START_REF] Tarantola | Generalized nonlinear inverse problems solved using the least squares criterion[END_REF]. Nevertheless, their efficiency depends strongly on the choice of the initial model as they only provide local minima. On the other hand, stochastic methods (such as Markov chain Monte Carlo McMC) are based on the concept of sampling where the solution is obtained by exploring the optimization performance of many randomly generated models and selecting the best of them in terms of minimization of the objective function. Most of the time, algorithms belonging to this category are easy to implement and insensitive to the initial model.

However, their applicability to inverse problems with a large number of unknown parameters remains limited due to the high repetitive and heavy computation of the forward problem during the sampling process [START_REF] Oliver | Markov chain Monte Carlo methods for conditioning a permeability field to pressure data[END_REF]Fu and Gómez-Hernández, 2009b;[START_REF] Wang | A hybrid inverse method for hydraulic tomography in fractured and karstic media[END_REF][START_REF] Jardani | Geostatistical inverse modeling of the transmissivity field of a heterogeneous alluvial aquifer under tidal influence[END_REF][START_REF] Elsheikh | Bayesian Reservoir History Matching Considering Model and Parameter Uncertainties[END_REF]. Other global optimization algorithms such as particle swarm optimization, genetic, simulated annealing can also be used to identify the best solution by iteratively evaluating the objective function until the state of convergence. However, convergence can be slow, especially in the case of high-dimensional inverse problems (Scales et al., 1990;[START_REF] Fernández-Martínez | Uncertainty assessment for inverse problems in high dimensional spaces using particle swarm optimization and model reduction techniques[END_REF].

In this paper, we test a new generation of deep learning algorithm that can be considered as a global optimizer in which prediction does not depend on an initial model and the computation of the sensitivity matrix as it is the case for gradient-based methods [START_REF] Zio | Approaching the inverse problem of parameter estimation in groundwater models by means of artificial neural networks[END_REF]. The concept of this approach is different from previous optimization methods as the process focuses on approximating the inverse function by finding a universal relationship linking the input and output data (in the case of hydraulic tomography between hydraulic pressure and hydraulic transmissivity). Indeed, the deep learning is a powerful generalizer of non-linear and complex functions by identifying a set of parameters such as weights and biases allocated to the hidden layer neurons that process the input data and link them to their corresponding output. The first generation of deep learning algorithms uses a set of hidden layers with a large number of neurons fully connected via high-dimensional weights, which complicates computation and requires a lot of time, especially for handling the imaging tasks [START_REF] Shen | A transdisciplinary review of deep learning research and its relevance for water resources scientists[END_REF]. However, the emergence of new deep learning architectures with Convolutional Neural Networks (CNN) concept that have been proven in high-resolution image processing have opened the way for broadening their applicability in the hydraulic characterization [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF][START_REF] Sun | Discovering State-Parameter Mappings in Subsurface Models Using Generative Adversarial Networks[END_REF][START_REF] Zhu | Bayesian deep convolutional encoder-decoder networks for surrogate modeling and uncertainty quantification[END_REF]. Indeed, the CNN concept has succeeded in reducing drastically the image processing time thanks to the convolution calculation, which allows to retrieve the local features of the image through small filters [START_REF] Indolia | Conceptual understanding of convolutional neural network -a deep learning approach[END_REF]. These filters cover the whole input image but their activation is done zone by zone with local connections between pixels.

This local convolution reduces the number of parameters to optimize during the learning process [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF]. Among the first attempts at adopting the CNN concept in the realm of tomography by inversion in geosciences, we cite: Sun (2018) trained the Generative Adversarial networks for linking the hydraulic head map and the spatial distribution of hydraulic conductivity. [START_REF] Zhu | Bayesian deep convolutional encoder-decoder networks for surrogate modeling and uncertainty quantification[END_REF] proposed a deep convolutional encoderdecoder network to reconstruct the image of hydraulic conductivity from hydraulic head map. In this article, we discuss the use of the convolution neural networks formed according to the SegNet architecture to process hydraulic tomography. This technique was initially designed for the semantic segmentation of objects learned on the images [START_REF] Badrinarayanan | Segnet: A deep convolutional encoderdecoder architecture for image segmentation[END_REF]). To illustrate the inversion principle with the CNN-HT structure, we organize this manuscript with the following outline: the first section will be devoted to the generation of data used in training, validation and test, and the second to the introduction of concept of the CNN-SegNet structure.

In the application sections, we discuss the relevance of the approach on theoretical cases under different conditions on training data with high and low resolution and contaminated or not by the noise. We also devote a section to the comparison of the results obtained with CNN-SegNet code with those determined with the Gauss-Newton algorithm.

Theoretical concept:

2.1.Preparation of Training data:

In deep learning algorithm, the first step concerns the construction of training data that will condition the accuracy of the result obtained from the predicted inversion operator (see figure 1). Given the impossibility of having a real dataset connecting the spatial distribution of hydraulic properties of aquifers and the piezometric responses, we use synthetic models in which the hydraulic parameter fields are generated geostatistically and their corresponding piezometric responses are obtained by solving the numerical groundwater equation in steady state (forward problem).

(

Γ  ∇ ∇ =   =   p 0 p T h q - h h x x m δ , ( ) ( ) . , , d F( ) :
) 2
where m is the negative of logarithm of transmissivity To predict an inversion operator using convolutional neural networks with the encoder-decoder architecture, we need to reshape the piezometric data into a field with the same dimension than the hydraulic transmissivity field in order to link tow images with the same size. To do this, we can use one of two methods:

(i)

The first consists in interpolating the piezometric data recorded on the observation wells for each pumping test in order to form images with the same resolution in the target transmissivity field. In this case, the input data are composed of a number of maps equal to the number of pumping tests. In our example, we have 5 pumping tests and we will then obtain 5 maps that will be treated jointly as 5 channels in convolutional neural networks.

(ii) This second is more compact than the first in which all piezometric data set recorded for different pumping tests are introduced by a single map using a projection operation derived from the Gauss-Newton formulation. This formulation results from the minimization of the objective function presented in the introduction section (

( ) 0 ∂ = ∂ m m Ψ
) [START_REF] Tarantola | Generalized nonlinear inverse problems solved using the least squares criterion[END_REF]:

( ) ( ) 1 T 1 1 T 1 1 obs i i 0 d m d m i 1 i i i i F( ) - + + - - - - + = - - +         m m m m m h J C J C J C C , (3) where i
J is Jacobian matrix with a size of N×M, (N and M are number of the data and the unknown parameters respectively). At the first iteration, usually we start with the prior model then:

( )

1 T 1 1 T 1 obs 0 d m d 1 0 0 0 0 F( ) - - - - + = - +     m m m h J C J C J C , (4)
and this equation can be reformulated under this following form:

obs 1 Λ = + m h β , (5)
where

1 T 1 1 T 1 0 d m d 0 0 0 0 F( ) β - + - - - = -    m m J C J C J C , (6) and 1
T 1 1 T 1 d m d 0 0 0 Λ - + - - - =     J C J C J C , ( 7
)
1 m is the hydraulic transmissivity model resulting from the first iteration; obviously, it does not fit the observed data, but it can be applied as a projection term for all piezometric data in order to convert them into matrices with same size as the transmissivity fields. Jacobian can be derived analytically using the Thiem's equation for confined aquifers to express the radial flow field and its adjoint state operator (for details, see Appendix).

To further simplify this projection term so that it is less dependent on prior information, we remove the bias term β , replace the covariance matrices m C and d C by simple identity matrices, now the new projection term takes this form:

p obs obs Λ = h h With 1 T T m 0 0 0 Λ - + =     J J I J . ( 8
)
Thanks to this simple multiplication, the hydraulic data obs h will be transformed into a vector The comparison of the results of the two methods (linear interpolation and projection) in terms of inverse operator prediction is presented in the last section. However, the first sections all applications are done with the projection method.

Inversion by Convolutional Neural Networks with encoder-decoder structure

In this article, we will explore for the first time the efficiency of SegNet architecture to process the inversion in hydraulic tomography. The approach uses the Conventional Neural Networks concept built with an encoder-decoder architecture (see figure 3). This approach was originally conceived to process semantic segmentation of images by delineating the shapes of learned objects [START_REF] Badrinarayanan | Segnet: A deep convolutional encoderdecoder architecture for image segmentation[END_REF]. from [START_REF] Badrinarayanan | Segnet: A deep convolutional encoderdecoder architecture for image segmentation[END_REF]. Here, the piezometric data can be resized to the same dimension of the output field using the projection method

In the encoder network, the input image is processed sequentially by a set of encoders. In our case, we use only 2 encoders chosen after an analysis of the inversion results obtained with networks having 1, 2 and 3 encoders. The network with 2 encoders gave the best result in terms of training and generalization with a low sensitivity to overfitting. Each encoder starts with a convolutional process in which multiple filters with the small size (3×3×64) are convoluted on the input layer zone by zone to identify their features [START_REF] Badrinarayanan | Segnet: A deep convolutional encoderdecoder architecture for image segmentation[END_REF]. This calculus permits to establish a local connectivity between pixels with these shared small filters, which remain simple to predict in the training process. The convolutional operations are followed by Batch normalization which is a linear operation that plays the role of a regularizer in scaling the convolution outputs in order facilitate and accelerate the training [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF]. This step is followed by the application of the ReLU function, which introduces the nonlinearity in the process. Then, the dimensionality of the result of the previous sequence will be considerably reduced in the max-pooling layer with stride of 2 and 2× 2 windows to keep only the main features. Therefore, the construction of feature maps in the encoder network leads to a decrease in the spatial resolution that will be restored in the decoder network. Indeed, the decoder network is configured in a symmetrical form with 2 decoders to prepare the output by recovering the resolution lost in the decoder with up-sampling operations from the max-pooling indices obtained in the previous max-pooling layers carried out in each encoder [START_REF] Badrinarayanan | Segnet: A deep convolutional encoderdecoder architecture for image segmentation[END_REF]. These results of up-sampling will be convoluted with other filters to build feature maps with more details. Then the Batch normalization and ReLU

1 i i 1 i ârg min F ( , ) ϒ ϒ - =   = -     ∑ Nt m h , ( 9
)
where i is the sample index of the training data set with Nt as size of the training data , p obs,

1 i F ( , ) ϒ -h
denotes the inverse function to identify by using the training data in the prediction of the networks parameters ϒ with ADAM optimization algorithm (Kingma and Ba, 2014).

This optimizer was run with a batch size of 70 and a learning rate starting at 0.01 and decreasing by 0.1 every 50 epochs. The computation was performed on a workstation with (Intel(R) Xeon(R) Silver 4110 CPU @ 2.10GHz) and 128 G of RAM with a single GPU. After training is complete, we use a sample of generated data not used during training phase to check the quality of the inversion results by the CNN-HT code.

Applications to synthetic cases

In this section, we apply the inversion concept with CNN encoder-decoder networks on theoretical cases where the log of hydraulic transmissivity fields are generated randomly with a Gaussian variogram and their hydraulic head responses by the groundwater equation as explained in the data preparation section. This distribution has -5 μ = 10 (m 2 /s) as mean and an isotropic Gaussian variogram with 0.5 and 10 m as variance and range respectively. With this variogram model, we performed 11000 realizations in which the transmissivity is ranging over [10 -8.3 , 10 -1.6 ] (m 2 /s); however 95% of these models have transmissivities between 10 -6.4 and 10 - 3.5 (m 2 /s) (see figure 4). Then, these models are used in a forward operator that takes 2s per simulation to get the hydraulic head responses. These realizations aim to evaluate the

Effect of training data size

In the first analysis, we examine the impact of training data size on the reliability of the inverse operator prediction and determine the amount of data required to obtain a satisfactory prediction. For this purpose, we performed three assessments with different training data sizes: T1=1000, T2=3000, T3=10000, the validation was done with 1000 models and the test with 500 models different from those used in the learning phase (training and validation). During the optimization process, we followed the quality of the minimization of loss functions with the training and validation data until they reached optimal values to stop the process manually or wait until the maximum number of epochs was met (see figure 5). In this way, we avoid the problem of overfitting which appears especially when we tried to build the network with a limited number of training data (T1=1000).

The results obtained with these different sizes of training data are presented in the Table 1a, in which we can easily identify through the analysis of the correlation coefficient between the real and predicted models of test sample that the accuracy of the prediction improved with increasing training data size. Indeed, the reliability of the predictive model derived from a deep learning algorithm is highly dependent on the amount of data used in the training, in which a large amount of data guarantees a good prediction. This implies a long computational time in the generation of training data, in any case the degree of accuracy desired remains linked to the nature of the problem studied. In our case, training with T1=1000 data allowed us to predict 500 test models with a correlation coefficient R between real and predicted models ranging from 0.74 to 0.93 and a mean of 0.86 (see Table 1a). From these models, we choose 3 models: the first one (m1) with a relative poor prediction (R=0.75), a second model (m2) (R=0.83) close to the mean representing quality of the majority of predicted models, and the last one (m3) with a good reconstruction (R=0.91) (see figure 6). These models will be used to analyze the impact of training data size on the quality of their estimations (see figure 6). The first model m1 is characterized by a high degree of variability with the small heterogeneities at the boundary of the domain and a logarithm of the hydraulic transmissivity that varies between [-6 , -2.5]. Over this range, there are outliers in the Gaussian distribution of training data in particular which are around -2.5. Thus, in this case only a few models have been seen in the training with this feature, which explains the poor reconstruction (see figure 6). The same problem occurs in the prediction of models with low outliers < -6.5. However, the models having variabilities well covered by training models such as m2 and m3 with transmissivities values between [-6.5, -3.5], the reconstructions are more accurate. However, if the number of training models increases, the models with outliers will be more frequent and will be reconsidered in the generalization as shown in the table, where the reconstructions of these models have been well improved.

Effect of amount of observation wells

In the second analysis, we study the influence of the number of hydraulic data collected during the pumping tests on the efficiency of the network to approximate the inverse operator.

In this analysis, we first reduce the number of observation wells to half (22 observation wells)

and in the second step to one quarter (11 observations wells). However, we maintain the number of pumping tests and the size of the spatial discretization of the transmissivity fields ( 16 1b). Because the number of piezometers becomes insufficient to cover all the heterogeneities and the complexities of the medium, and this also occurs with classical inversion methods. This is further confirmed by analyzing the spatial aspect of the reconstructions of models (m1, m2 and m3), in which we observe a loss of resolution and an increase in smoothness degrees with the reduction of hydraulic data (figure 8). However, the main heterogeneities can still be identified with less data.

Effect of observation uncertainties

In this section, we analyze the effect of uncertainties that may be associated with hydraulic pressure measurements collected during pumping tests. This involves contaminating the numerical piezometric data with a random Gaussian noise of 0.5 as the standard deviation, and verifying its impact on the quality of hydraulic field reconstructions of test models.

As a first analysis, we use an inverse operator established on the basis of uncontaminated training data to interpret hydraulic head data from a test set that are contaminated by noise.

Here, the training operation did not take into account the noise. The result of this test is presented by the benchmark models in figure 9 and table 1c of correlation coefficients, which reveal a sharp degradation in the quality of the predictions compared to the noise-free test models, with a decrease in the correlation coefficient (from 0.92 to 0.79 mean for all test models). In the second strategy, noise is taken into account in the estimation of the inversion function by adding noise to the hydraulic data used in the learning process. This method has led to a significant improvement in the prediction quality of the test models, as shown by the three models (m1,m2 & m3) and the correlation coefficient noted in table 1c with the average R of the set of models increasing from 0.79 to 0.84. However, this prediction is less accurate when compared to the noise-free training and test data, as shown by the R-values. We can, therefore, conclude that uncertainties in the measured data may have an impact on the results of inversion, as is the case with classical optimization techniques, but to reduce their impacts, it is necessary to incorporate these uncertainties in the construction of the inversion model, i.e. in the learning process.

Comparison of CNN-HT and Gauss-Newton algorithms

In this section, we compare the inversion results obtained with the CNN-HT algorithm with those determined by the traditional Gauss-Newton (GN) method. To this end, we apply both algorithms to the reconstruction of a transmissivity field with a mean of 10 -5 m 2 s -1 and generated with a spherical covariance having 0.5 and 10 m as variance and range respectively.

However, this time, the models are highly discretized (32×32) compared to those used in the previous sections. The Gauss-Newton algorithm, as described in the introductory section, is a deterministic approach that relies on the minimization of an objective function (Eq. 1) by iteratively computing the Jacobian matrix until obtaining a local minimum. The calculation of the Jacobian matrix remains the most demanding part of this process, especially when the number of unknowns is important. This calculation is usually done either by the finite difference or the adjoint state technique. The finite difference method is simple to implement, but its computation is time consuming because it involves to solve the forward problem for each unknown parameters [START_REF] Cardiff | Efficient solution of nonlinear, underdetermined inverse problems with a generalized PDE model[END_REF][START_REF] Fischer | Application of large-scale inversion algorithms to hydraulic tomography in an alluvial aquifer[END_REF]. In our case study, where the model is discretized into 32×32 cells, the construction of the Jacobian will therefore involve solving the forward problem 1024 times. On the other hand, the adjoint state technique requires less computation, by iteratively solving an adjoint operator with the form of the forward problem on the number of observation wells, but its implementation is quite complex (for more details see [START_REF] Sykes | Sensitivity analysis for steady state groundwater flow using adjoint operators[END_REF][START_REF] Cardiff | Efficient solution of nonlinear, underdetermined inverse problems with a generalized PDE model[END_REF]. In this section, we test the implementation of Gauss-Newton algorithm to invert the synthetic hydraulic data determined with a hydraulic transmissivity field shown in figure 10, (called True model) and using the piezometric coverage shown in figure 2. These data were then contaminated with noise with a variance of 5cm 2 . The inversion of these datasets was constrained, as explained in the introduction, by geostatistical parameters such as covariance and mean. The results of the inversions with the GN algorithm using the finite difference and adjoint state methods are presented in Figure 10 where we can conclude that the GN algorithm succeeded in both cases in identifying the main heterogeneities of the real field with certain smoothness. In the case where the Jacobian is computed with the finite difference method, convergence was achieved after 7 iterations in 4h10min and solving the forward problem 7×(1024+1) =7175 times. In contrast, the adjoint state method took only 1min24s.

Regarding the application of the CNN-HT algorithm, we geostatistically generate 11000 transmissivity models with the same variogram as the real model used for the comparison. We keep 10000 models for the training, and the rest of the models are equally split for validation and testing. In this case study, we test two different ways of resizing the piezometric data, the first method called projection method uses the Hessian matrix to reshape the input data into the same size as the transmissivity field (output) as described in the section 2. This method requires the calculation of a Jacobian matrix with a homogeneous field for example the mean field value.

The second method is based on a linear interpolation of the hydraulic data recorded during each pumping test to form 5 maps of size 32×32 (i.e., one map for each pumping test), so the input data has the form of a matrix of 5 channels (32×32×5).

Both approaches were used to build the SegNet networks and then tested on the test sample. The analysis of the predictions obtained on this sample shows that both methods provide almost the same prediction qualities with a small advantage for the projection method according to the correlation coefficients (see table 2). However, when we reduced the training size to T=500 and T=2500 models, the network using the projection method was able to provide accurate generalization while the network with interpolated data failed. We believe this is due to the fact that the interpolation of piezometric data does not preserve all the information masking the degree of variability in hydraulic responses, which makes learning difficult with few training data. In addition, using the interpolation method increases the size of the input data, which requires more memory and makes the training process more cumbersome compared to the projection technique that brings the data into a compact form.

Regarding the application of these two networks on the hydraulic data of the real model, we find that the network built with compact data provides a better reconstruction of the spatial variability of transmissivity field with a correlation coefficient R=0.80 (see figure 10).

Moreover, this prediction has almost the same accuracy as the result determined by the Gauss-Newton algorithm (R=0.83). We also mention that the time used in the GN with the finite difference method is perhaps comparable to the time spent in the DL method mainly devoted to the construction of the data set. In particular, the DL method becomes faster when we reduce the amount of training data to T=500 and T=2500 (20 min and 1h41 min), because even with these small samples, the transmissivity field can be well reconstructed well in a reasonable amount of time, as demonstrated by the comparison of real and inverted fields (figure10).

At the end of this comparison, we conclude that both the CNN-HT and Gauss Newton codes are based on the forward problem for which the computation time depends on the numerical method and the degree of mesh refinement used to solve the groundwater flow equation. The machine learning technique uses this forward problem to generate a training database to build the CNN network. The time spent in this construction can be equivalent or less to the time spent in finding a local minimum with the Gauss Newton using the finite difference method to calculate the Jacobian matrix over a large number of unknown parameters.

However, it is difficult to have an idea in advance of the time that will be spent in building this database because the size needed to get accurate results is related to the degree of complexity of the input-output relationship and the type of network. The optimal size as well as the rest of the network parameters (epochs, number of encoders) are only established by a trial-and-error analysis, which can take a considerable amount of time. Therefore, the data size used in this paper cannot be generalized to all cases treated by this network or others.

Both methods also share their dependence on prior information. In the Gauss-Newton algorithm, this prior information is introduced as a constraint to guide the optimization towards realistic models and reduce the uncertainties associated with the non-uniqueness of the solution.

On the other hand, the CNN-HT algorithm uses it to generate the training models. The CNN-HT offers the possibility of incorporating multiple prior models into the construction of the training dataset, which is difficult to perform in a deterministic algorithm where the mathematical forms of these constraints must be differentiable. In addition, both methods are sensitive to the uncertainties associated with the measured data. The CNN-HT allows simultaneous learning with data affected by several degrees of uncertainty by adding noise to the training data, then combining this set with the original to double the size of the training data and take into account the impact of noise. Regarding the quantification of the prediction uncertainties, in the Gauss Newton method, we use the a posteriori covariance formulation to get an idea of the uncertainty interval concerning the local minimum [START_REF] Tarantola | Generalized nonlinear inverse problems solved using the least squares criterion[END_REF]. However, in the CNN-HT algorithm it is difficult to construct an uncertainty map of the predicted model. We can only build an uncertainty map on the model used in the test.

Table 2: We report the metric evaluations (correlation coefficient R) and computational time of inversions obtained with the CNN-HT and Gauss-Newton algorithms. The CNN-HT network was trained with hydraulic data introduced either by the linear interpolation method or by the projection method. We tested 3 different sizes of training data (500, 2500, 10000), but the network with the interpolation method failed in the training process when the data size became small (500 and 2500). We also present the evaluation of the Gauss-Newton metric using the finite difference method and adjoint state on a test model shown in figure 10.

Conclusion and Summary

In this work, we explore the relevance of a deep learning tool called the SegNet network in the mapping of spatial variability of hydraulic parameters by interpreting the piezometric data recorded during pumping tests. This approach relies on the approximation of the non-linear inverse operator connecting the hydraulic head data to the hydraulic transmissivity field. This approach is the outcome of an adaptation of the SegNet network that has been successfully applied in the realm of image segmentation by identifying the objects studied on the pixels of images. The SegNet is designed to establish a relationship between two images sharing the same size and certain similarities, which is not the case for hydraulic head data, which are less numerous than the hydraulic transmissivity values to be estimated. To meet this point, we test the projection and interpolation methods. Projection method is based on the classical iterative formulation of a least-square method at the initial iteration to build a projection operator with the Hessian and Jacobian matrices derived from the mean value of transmissivity fields used in the training data to resize the hydraulic data into the hydraulic transmissivity field size.

Interpolation method involves reshaping the data size by interpolating the hydraulic head data to form a map with the same size as the output (transmissivity field) for each pumping test. As a result, the interpolation significantly increases the size of the input data, which can slow down the numerical calculation without improving the quality of the prediction compared to the projection method, which better preserves the hydraulic information in compact form.

In the projection method, the operator is applied to the totality of the hydraulic pressure data used in training which are determined numerically as responses to pumping tests carried out on geostatistically generated hydraulic transmissivity fields with certain statistical properties. Once the training data has been gathered, the inverse operator is established by identifying the different filters used in the SegNet, which is composed of encoders and decoders networks. The encoder network allows, through multiple application of filters and sequential calculations of convolution, Batch normalization, ReLU and max-pooling operations to extract the main features of the input data. This encoder process induces a drop in spatial resolution, which will be recovered in the decoder network, which is designed to prepare the output with the applications of up-sampling, convolution, Batch normalization and ReLU operations. This network ends with a regression layer to evaluate the performance of the inverse operator approximation.

The quality of this approximation is highly dependent on the amount of data used in the training with an accurate approximation when the dataset is large. This makes the assembly of the training data the most expensive step in time, requiring several hours compared to the learning process that in this case does not exceed 30 minutes. In practice, it is preferable to launch the learning with few data and check whether or not the desired accuracy is achieved before generating large amounts of data unnecessarily. In terms of comparison, the time spent in forming a training database can be equivalent to the time spent to reach a local minimum with Gauss Newton using the finite difference method to compute the Jacobian matrix

In addition to the amount of training data, the nature of this data also controls the prediction results of the network. Indeed, the trained network is only intended to handle certain types of models with similar characteristics to the models used in training. Thus, it is by no means a universal inversion operator, which makes the choice of the nature of the training models a crucial step that must be based on realistic prior information. In real application case, it should also be verified that the transmissivity fields used in the training allow to generate hydraulic data with the magnitude of the piezometric data observed in the field. It is also necessary to ensure that the degree of heterogeneity of the generated models can be easily mapped to the piezometric data during training process. Because if the piezometric coverage does not allow to cover transmissivity models with high heterogeneity, the learning process will fail to establish this relation. The result of the inversion with CNN-HT depends on the prior information and this particularity is also found with classical inversion techniques such as Gauss-Newton.

Noise can alter the information carried by the piezometric data and is a frequent problem in inversion algorithms. The deep learning algorithms are no exception, but its impact can be reduced by integrating it into the learning step by contaminating the input of training data. The number of wells and their spatial configurations used in the study also had a determining effect on the effectiveness of the CNN-SegNet inversion tool, with more reliable reconstructions when these piezometers are well distributed over the study area, which was also observed in the use of conventional deterministic and stochastic methods.

At the end of this conclusion, we enumerate and discuss the main advantages and limitations of the proposed method. Regarding the advantages, we point out that:

The CNN-HT method uses the forward problem as a black box to generate the training data with user-selected transmissivity fields, thus avoiding numerical instability and outliers in hydraulic head simulations related to its solving. These numerical problems can be encountered when using conventional inversion methods (stochastic and deterministic) if they are illconstrained.

Once the inverse operator is trained, the network permits to interpret the piezometric data in a few seconds without the user having to add new parameters.

The CNN-HT code provides predictions that require any initial hydraulic transmissivity models or computation of sensitivity matrix, as is the case with deterministic methods. This advantage makes CNN-HT an effective inversion tool for dealing with highly nonlinear inverse problems with large-scale parameterization that are difficult to handle with deterministic codes, especially when the computation of sensitivity is complex.

The CNN-HT code performs reconstructions of the transmissivity field that depend on the a priori model used in the generation of transmissivity field for the training data. This dependence to the a priori model is a property shared with other inversion methods. However, in this CNN-HT tool, the mathematical form of this a priori model is not subject to any condition, which is not the case with deterministic methods where the model must have a quadratic and differentiable form.

The CNN-HT algorithm also provides great flexibility for the user to incorporate multiple prior models, especially when in doubt about which a priori information to use.

The proposed method is based on a simple network that can be generalized and applied to handle hydraulic tomography in 3D or in time domain in order to identify the storativity coefficient as well, by only modifying the training dataset. It is also easily adaptable to deal with hydraulic tomographies with complex heterogeneities having for example multiple discontinuous hydrofacies with contrasting transmissivities, or parameterized in discrete fracture network mode. These types of problems can be treated as a classification task by replacing the regression layer of the network with the Softmax layer. In this way, the network can identify fracture geometries or hydrofacies shapes.

The code allows the incorporation of uncertainties on the boundary conditions imposed in the numerical solving of the forward problem by randomly choosing their values within a confidence interval when building the training data set.

This code, like other deep learning algorithms that use the convolutional neural network technique are able to better reduce the effects of noise in the data. This is thanks to the convolution concept that processes the input data with some local connectivity and not in a pointwise manner as in other conventional inversion methods [START_REF] Vu | Convolutional neural networks with SegNet architecture applied to three-dimensional tomography of subsurface electrical resistivity: CNN-3D-ERT[END_REF]. The impact of noise could be minimized and accounted for in the inversion results by contaminating the hydraulic data used in the training data with a noise signal.

Regarding the disadvantages of this method, we mention that It is based on a repetitive and time-consuming numerical resolution of the forward problem, which is involved in the construction of the training database. The size of this database conditions the quality of the predictions, and its optimal size can only be obtained by a long trial-and-error procedure.

The method requires important computing resources to carry out the learning.

The code also uses a large number of parameters that intervene in the construction of the networks (such as: size of the decoder-encoder, number of filters and their size), and in the optimization algorithm used in the training phase (number of epochs, regularization parameters, initialization of filters, Batchsize, Dropout). The determination of all these parameters implies a tedious sensitivity analysis with several training operations to choose the best parameters.

The method does not provide an uncertainty assessment of the predicted transmissivity fields. This is a major weakness of this method compared to classical approaches. There are some attempts in the literature to estimate the uncertainty of the prediction determined with DL algorithms by performing the learning in a probabilistic way [START_REF] Kendall | Modelling uncertainty in deep learning for camera relocalization[END_REF][START_REF] Gal | Dropout as a Bayesian approximation: Representing model uncertainty in deep learning[END_REF]. However, this type of approach requires a lot of computational time and it does not take into account all the uncertainties that can come from all the parameters used in the learning.

The method is highly dependent on prior information. This dependence can be an obstacle to its application, especially when the network fails in the learning phase to relate the transmissivity fields to their piezometric data. This can occur when the transmissivity fields have complex heterogeneities that are not captured by low piezometric resolution. This dependence can also lead to a misinterpretation, even with a well-trained inversion operator, when the magnitudes of the piezometric data used in the training are so different from those we desire to predict. To avoid both of these problems, it is necessary to ensure that generated transmissivity fields have a degree of heterogeneity that can be mapped with the piezometer coverage and that they produce hydraulic head data that include the magnitude of the hydraulic observations to be interpreted.

In conventional inversion methods, the predicted transmissivity distribution is the result of a compromise between the a priori model and the piezometric data. However, here, the result is strongly linked to the a priori model and the piezometric data to be interpreted come only at the end of the learning process.

. x y with a rate p q . φ is the adjoint operator with a source term

  [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF] used the Generative Adversarial networks to generate randomly in 2D and 3D binary hydraulic conductivity fields and MCMC to perform the inversion process. In geophysics realm, Wu and Lin et al. (2019) applied CNN with encoder-decoder networks for mapping the subsurface velocity from seismogram data. Puzyrev and Swidinsky. (2019) used the CNN network to determine the vertical subsurface heterogeneity of electrical conductivity by training electromagnetic data.

  studied domain. h is the hydraulic head responses due to water extraction represented by a point source term well locations p x ( δ is Dirac delta function). The pumping tests will be done sequentially on several wells. The forward problem is solved by finite element technique using the Comsol software with a constant hydraulic level 0 h imposed as Dirichlet boundary at all limit of the domain d Γ that are positioned far from the studied area for reducing their impact.

Figure 1 :

 1 Figure 1: A sketch describing the concept of hydraulic tomography with the CNN-SegNet. This concept is based on the prediction of inverse operator

Figure 2 :

 2 Figure 2: Position of the wells used in the generation of the hydraulic head data. 5 pumping tests were carried out successively in the wells marked with a triangle and the responses of the hydraulic head were collected in 44 wells with the dot symbol. This configuration provides 220 hydraulic data that are then used in the reconstruction of the hydraulic transmissivity.

m

  representing the mean of the fields. In this case, the

h

  of the same size as the output hydraulic transmissivity field and used in matrix form as the input field in the network. This formulation allows to link linearly obs h

Figure 3 :

 3 Figure 3: SegNet Convolutional Neural Networks with encoder-decoder structure. Modified

  are applied successively on the convolutional outputs as in the encoder network. At the end of the encoder network, we add a convolution layer having the same size of the output and a regression layer to assess the performance of the training. This last layer replaces the Softmax layer used in the original version of SegNet network for the segmentation task. The set of weights and biases ϒ used in the various encoder-decoder operations are chosen in way to match the training data p

  effectiveness of the CNN-SegNet approach in the reconstruction of hydraulic transmissivity field and to analyze the sensitivity of this reconstruction to the size of the training dataset and the amount and uncertainties of the hydraulic head data collected during the pumping tests.

Figure 4 :

 4 Figure 4: Gaussian distribution of the log of transmissivity fields used in the training of CNN-SegNet network.

Figure 5 :

 5 Figure 5: Both curves represent the quality of predictions of the network on validation and training data (10000 models for training, 1000 models for validation) during the optimization process.(RMSE: Root Mean Square Error).

Figure 6 :

 6 Figure 6: Three models selected from test models to illustrate the prediction quality of the CNN-HT tool by comparing their true and predicted transmissivities. In this figure, the

  ×16), as used in previous test. These new well configurations are illustrated in the figure 7. In practice, only the hydraulic responses of abandoned wells are removed from previous training data. Thus, these data acquired with the two new observation well configurations are used as training input to determine the corresponding inverse operators, and the results obtained were compared to the previous configuration with 44 wells. The training operations were performed with 10,000 models, and the tests were also done on the previous unseen samples. The reconstruction results of the test models with different series of piezometric data expressed in terms of correlation coefficients reveal a degradation of the accuracy of prediction with the decrease in the number of hydraulic data (see Table

Figure 7 :

 7 Figure 7: Two new observation well configurations were retained to study the effect of the reduction in hydraulic data on the quality of the prediction. In the first configuration, the number of observation wells was reduced by half (22 observation wells) and in the second, by a quarter (11 observation wells) from the initial configuration shown in Figure 2.

Figure 8 :

 8 Figure 8: Predictions of the three models are provided with a set of training data (T3=10000) in which piezometric responses are collected with 3 different well configurations (44, 22 and 11 wells). These predictions showed a degradation in the reconstructions of hydraulic transmissivity with a decrease in the number of wells used in the survey.

Figure 9 :

 9 Figure 9: The purpose of these reconstructions is to analyze the effect of noisy data on the efficiency of the inversion with CNN-HT by comparing 3 types of predictions: 1) with the training and test data without noise; 2) with the training data without noise but the test models are noisy; 3) both the training and the test models are noisy. These comparisons show that the uncertainties in the measured data are affecting the quality of the reconstructions; however, their impacts could be attenuated by taking into account this noise in the training process.

Figure 10 :

 10 Figure 10: The results obtained when applying CNN-HT and Gauss-Newton algorithm in the reconstruction of the "True model". The results obtained with GN using the finite difference and adjoint state methods recovered the main heterogeneities of the field. The same reconstruction quality was obtained with CNN-HT, especially using the projection method for the preparation of the input data. Even with a small amount of data ( T=500 and T=25000), this CNN-HT network provided satisfactory reconstructions.

  T is transmissivity hydraulic with m=-log(T).When the transmissivity is constant, h and φ can be determined analytically with Thiem'values imposed as the boundary of the domain located at a distance R, away from the study area. r is radial distance. Then the sensitivity of hydraulic head at point O to a variation of the transmissivity at cell Ω j can be approximated with this integral:

Table 1

 1 Summarizing metric evaluations of the predictions quality (in term of correlation coefficient) of 500 test models with correlation coefficient. a. Evaluation of the impact of the size of training data (T) on predictions. b. Evaluation of the impact of the number of observation wells on predictions. c. Evaluation of the effect of uncertainties in the piezometric data on predictions with incorporation or not of these uncertainties in the learning process.

	a. Size of training data	Correlation coefficient
			R
		min	mean	max
	T1=1000	0.74	0.86	0.93
	T2=3000	0.77	0.90	0.95
	T3=10000	0.84	0.92	0.96
	b. Number of observation data	Correlation coefficient
			R
		min	mean	max
	T3=10000, N° data=220	0.84	0.92	0.96
	T3=10000 N° data=110	0.73	0.85	0.93
	T3=10000, N° data=55	0.57	0.77	0.91
	c. Degree of noise on the data	Correlation coefficient
			R
		min	mean	max
	T3=10000, N° data=220,	0.84	0.92	0.96
	std=0 free noise
	T3=10000, N° data=220,
	std=0.5 on Training and Test	0.55	0.84	0.93
	data
	T3=10000, N° data=220,	0.49	0.79	0.92
	std=0.5 only on Test data

Acknowledgements :

We thank Normandy region for its financial support to our consortium on the hydraulic characterization of aquifers.

Appendix

Here we present the mathematical formulation to analytically calculate the Jacobian matrix of a homogeneous transmissivity field. This matrix is used in the projection method to transform the hydraulic data to the same size as the transmissivity field to form the Segnet network.

We recall that in the joint state method, the sensitivity matrix can be expressed as an integral [START_REF] Sykes | Sensitivity analysis for steady state groundwater flow using adjoint operators[END_REF]: