# Multi-task neural network in hydrological tomography to map the transmissivity and storativity simultaneously: HT-XNET

M. T. Vu, Abderrahim Jardani

► To cite this version:

M. T. Vu, Abderrahim Jardani. Multi-task neural network in hydrological tomography to map the transmissivity and storativity simultaneously: HT-XNET. Journal of Hydrology, 2022, 612, 10.1016/j.jhydrol.2022.128167. insu-03779401

HAL Id: insu-03779401

https://insu.hal.science/insu-03779401v1

Submitted on 25 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

# Multi-task neural network in hydrological tomography to map the transmissivity and storativity simultaneously: HT-XNET

2

3

M. T. Vu ¹,*, A. Jardani ¹

4

¹ Université de Rouen, M2C, UMR 6143, CNRS, Morphodynamique Continentale et Côtière,

5

GeoDeepLearning Consortuim, Mont Saint Aignan, France

6

* Corresponding author: minh-tan.vu@univ-rouen.fr

7

Emails: abderrahim.jardani@univ-rouen.fr; minh-tan.vu@univ-rouen.fr;
**Abstract:**

In this paper, we present a novel method to simultaneously map transmissivity and storativity in a heterogenous aquifer using a multitask convolutional neural network. This coupled inversion algorithm translates transient hydraulic head data from pumping tests into two independent tomography of transmissivity and storativity in a two-dimensional problem. Based on the SegNet architecture, the multitask neural network provides an effective solution with minimal weights and biases, works as an end-to-end operator that directly approximates the inverse function. Multiple sharing mechanisms enable the multitasking approach to outperform single-task models by 10% accuracy. Application to synthetic experiments shows that the quality of inverted maps relies on the complexity level of heterogeneity in the hydraulic fields. Mapping accuracy also depends on the data coverage, but being resistant to data noise due to the featuring mechanism in the convolutional network adopted. Comparison with conventional inversion method reveals that the deep learning approach provides better generalization in the inversion while performing each inference on the order of milliseconds without further calibration from users.

Keywords: Neural Network, Multi-task, Groundwater, Inversion, Transmissivity, Storativity.
23

**Introduction**

24

25

26

27

28

29

30

31

Understanding flow and contaminant transport mechanisms in water reservoirs requires reliable assessments of the spatial patterns of hydrogeologic properties, including hydraulic transmissivity and specific storage coefficient (Fitts, 2012; *Şen, 2015*). Methods for evaluating these parameters include hydraulic tomography, which involves the acquisition and analysis of multiple piezometric responses associated with multiple pumping tests conducted sequentially. Hydraulic tomography analysis often involves a joint interpretation of all hydraulic head measurements in which the inversion approach is used to find the best transmissivity and storativity fields in terms of numerical fitting of the hydraulic measurements (e.g., Berg & Illman, 2015; Castagna, et al., 2011; Franssen, et al., 1999).

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

The inverse problem is by nature ill-posed and requires additional so-called *a priori* information about the properties of the unknown hydraulic parameters to narrow the field of possibilities and arrive at a solution that satisfies both the measured data and *a priori* data criteria. This ill-posed problem is mainly due to the fact that the measured piezometric data are not sufficient to identify all heterogeneities of the environment (Tarantola, 2005; Yeh, 1986). Among the approaches used to constrain the inverse problem is the geostatistical approach. In which, the statistical properties (mean and covariance) of the target field are included in the optimization process as known parameters to regulate probability distributions of transmissivity and storativity (Tarantola, 2005). In optimization, deterministic algorithms such as Gauss-Newton or Conjugate Gradient, are usually applied to access a local minimum in handling the hydraulic tomography inversion. The algorithms can incorporate a large dimension by computing the sensitivity of all measured values with respect to a small parameter variation at each unknown transmissivity and storativity parameter assigned to the discretized cell of the study area. This computation can be effectively performed using an adjoint state formulation (Sun & Yeh, 1990), where an adjoint partial differential equation system further needs to solve for each observation point. However, in the transient mode, the inversion of hydraulic tomography becomes more complex, as the optimization aims to determine the spatial variabilities of both the transmissivity and storativity variables from a time-recorded hydraulic head data set. Under these conditions, the number of
parameters increases, making the problem more ill-posed and underdetermined in practice (Tarantola, 2005; McLaughlin & Townley, 1996).

In this study, we develop a deep learning algorithm based on the concept of convolutional neural networks to directly approximate the inverse function. A function that performs an end-to-end operator to interpret the input hydraulic data as an image into multiple output hydraulic tomographies: the transmissivity and the storativity. A neural network uses multiple linear and nonlinear operators executed sequentially on hidden layers, involving multiple multidimensional coefficients. This high-dimensional parameter system allows us to efficiently reconstruct the highly nonlinear inverse operator in our problem (Hornik, et al., 1989). The coefficients are derived via an optimization process known as the learning task. This learning task consists of determining the coefficients that process the input maps to match the corresponding output maps, both of which are stored in the database. Recent advances in learning algorithms have helped facilitate this optimization, while most of the effort has been devoted to building a generalized database (Schmidt et al, 2019; Wani et al, 2019). Since the hydraulic tomography database can only be built based on a theoretical approach, in this approach we generate thousands of transmissivity and storativity fields that are used to numerically generate the piezometric responses related to the pumping tests. The entire dataset is thus used to build a network linking the piezometric responses and the hydraulic properties fields.

Such deep-learning algorithms are becoming an essential modeling tool in many fields, including the geosciences, thanks to the development of new networks such as convolutional neural networks. The introduction of convolutional neural networks (CNNs) has offered an efficient alternative to solve complex image processing problems such as tomography, which are common in geophysical, hydrogeological or remote sensing applications. Their effectiveness lies in convolutional computation, which consists in extracting the main feature of the input images by scanning them with a small, shared filter to reconstruct the output image. Processing only the main feature with small filters makes the computation more flexible and facilitates learning in the training phase. CNNs are typically formed as a multi-layer architecture with interconnection schemes to make the network more efficient in handling
complex tasks (Albawi, et al., 2017). Among them, the encoder-decoder networks such as Segnet (Badrinarayanan, et al., 2017) and Unet (Ronneberger, et al., 2015) have been effectively applied in various tomography tasks in geosciences. In seismic studies, CNN networks have been used to reconstruct the spatial velocity of geological structures by inputting seismogram images into the networks (Zhao, 2018; Zhang & Lin, 2020; Yang & Ma, 2019; Moseley, et al., 2020). They were also adopted to perform 2D or 3D imaging of the subsurface by processing the apparent resistivity data recorded during the electrical resistivity tomography investigations (Liu, et al., 2020; Vu & Jardani, 2021). In hydrology, the CNN with encoder-decoder has been effectively adapted to hydraulic tomography to process hydraulic data associated with pumping tests in reconstructing the transmissivity field in porous aquifers and in delineating the fracture structure in fractured aquifers (Jardani, et al., 2022; Vu & Jardani, 2022). Within the framework of this topic, other CNN architectures were also tested such as the generative adversarial network (GAN) to image hydraulic parameters in 2D or 3D binary models from the inversion of hydraulic data (Laloy et al., 2017 & 2018; Sun, 2018).

Proposed by (Badrinarayanan, et al., 2017) for segmentation tasks, SegNet architecture has been quickly adopted in several applications in geosciences such as the classification of seismic facies, hydraulic tomography in the reconstruction of hydraulic transmissivity, detection of geochannels and mapping of 3D electrical resistivity (Pham, et al., 2018; Mukhopadhyay & Mallick, 2019; Vu & Jardani, 2021). While these works focus in single-task networks, the architecture proposed in this study attempts to deal with multiple tasks in a joint manner. When dealing with a multivariate problem consisting of multiple output variables, such as hydraulic tomography where both transmissivity and storativity fields are determined, the interpretation of each field can be made independently using a single-task neural network. In this case, each network is trained separately to derive the corresponding single field, which can be cumbersome. However, the issue can be addressed jointly in a single neural network, a multi-task algorithm. Multi-task learning (MTL) involves designing a network that can perform multiple tasks simultaneously, in which each task benefits from the others through mutual learning. The network interconnection can be established at hierarchical levels, through which the extracted features are mutually transferred, acting as data augmentation during the training process. By exploiting the
similarity and contrast between mutual features, MTL-based strategy promotes the generalization, efficiency, and accuracy of each task (Zhang & Yang, 2017; Ruder, 2017).

The effectiveness of MTL networks has been recently studied especially in image processing and text translation (Zhang & Yang, 2017; Liu et al., 2016). However, their application in geosciences is still limited and at an early stage of development. In the field of seismic exploration, Siahsar et al. 2017 developed a multi-task algorithm to reduce the effects of random noise on seismic data. Wu et al. 2019 suggested a Unet-based algorithm to detect faults, seismic images, and seismic normal vectors from training images of synthetic 3D seismic map. Li et al. 2020 used a multitasking neural network to improve the resolution of seismic tomography by pairing an encoder-decoder with a dense-skip connection structure. In this work, we develop a multi-task architecture to perform hydraulic transmissivity and storativity field reconstructions by analyzing hydraulic head data acquired in pumping experiments, which has rarely been discussed and used for inversion problems in hydrogeology. The proposed neural network consists of two encoder-decoder structures, both of which map transmissivity or storativity fields from temporal measurements of piezometric level related to pumping tests. Both architectures are based on SegNet, which is known for its ability to identify objects learned from images through segmentations with low computational cost compared to other autoencoders (Badrinarayanan, et al., 2017). In order to enhance the performance of multi-task modeling, each task element is previously optimized and then interconnected at multiple hierarchical levels.

This study, which focuses on the ability of developing an MTL network to perform inversion tasks mapping hydraulic properties of aquifers from hydraulic data, will be organized into 4 sections. Section 2 is dedicated to describing the theoretical background, including groundwater modeling, solving the inversion problem with the multitask network and its application to a synthetic case. Section 3 focuses on evaluating the performance of the trained network and the impact of various probable conditions including the quality and quantity of the trained data. Section 4 provides a comparison with a conventional Least-square inversion method, outlining the differences between the two methods and
their impact on prediction accuracy. Finally, we conclude the paper with the main results, discussions, and an outlook on future work.

## 2. Methodology

In this section, we begin by summarizing the theoretical concept of groundwater flow in a heterogeneous aquifer, often referred to as the forward problem, to numerically simulate hydraulic responses to pumping tests. In this approach, the forward problem is applied to form the training data that will be used to train the multi-task networks that will be described later.

### 2.1 Forward problem: Groundwater modelling

The horizontal groundwater flow in a 2D confined aquifer is governed by (Bear, 1972)

$$S \frac{\partial h}{\partial t} - \nabla(T\nabla h) = q \delta(x), \qquad (1a)$$

subjected to the following boundary conditions

$$h = h_D \quad at \quad \Gamma_D, \qquad (1b)$$

where $h$ (m) is the hydraulic head, $T$ (m²/s) is the transmissivity, $S$ (.) is the storage coefficient, $q$ (m²/s) is a sink/source at the well location $x$, and a constant hydraulic head $h_D$ is imposed as a boundary condition on the far-away boundary $\Gamma_D$.

Solving the forward problem with heterogeneous distributions of hydraulic parameters can be done only by a numerical method such as the finite element method. The solution involves in determining the spatio-temporal distribution of the hydraulic load related to a water withdrawal imposed on the borehole by introducing the hydraulic properties of the aquifer and the boundary conditions. Solving the flow equations in this work is carried out using COMSOL software (COMSOL Multiphysics) to create a database for training the neural network. In fact, the database consists of thousands of T and S fields randomly generated with characteristics regularized by a priori information. Such an a priori model, usually constrains the inverse problem to avoid non-uniqueness in optimization, is recalled here in the machine learning algorithm to form the generation of training database. As geostatistical
constraint is one of the most common methods in the inversion of hydraulic tomography, it is employed
here to generate the fields of hydraulic properties while the corresponding hydrological responses are
obtained from computing the forward problem.

## 2.2 Inversion problem: Hydraulic tomography from GW level measurements

The forward problem from the previous section can be reformulated as follows

$$h = f(S, T), \qquad (2)$$

where $h$ represents the hydraulic head data, $f$ denotes the forward operator, and $S, T$ are the storativity and transmissivity fields. In practice, however, the goal of a hydrogeological investigation is usually to map the parameter field from the observed data

$$(S, T) = f^{-1}(h), \qquad (3)$$

To estimate the parameter fields $(S, T)$ from a set of measurements, a conventional inversion solver attempts to fit the given $h$ by iteratively adjusting $(S, T)$ through an optimization process

$$(S, T) = \operatorname{argmin}(\|h - f(\bar{S}, \bar{T})\|) + g(\bar{S}, \bar{T}), \qquad (4)$$

where $(\bar{S}, \bar{T})$ denotes the estimated variable field $(S, T)$ in each step, while $g$ represents the regularization terms.

In this work, however, we design a neural network to directly approximate the operator inverse problem that linking the hydraulic head data and aquifer properties with the use of a training database. Then the approximation of this inverted function can be formulated as follows

$$f^{-1} = \Psi(\Theta), \qquad (5)$$

where $\Psi$ represents the neural network consisting of successive linear/nonlinear operators (the layers) with a set of parameters $\Theta$. Thus, $\Psi$ acts as an end-to-end function that directly reconstructs the storativity and transmissivity fields from a given measurement of $h$

$$(S, T) = f^{-1}(h, \Theta) = \Psi(h, \Theta), \qquad (6)$$
Determination of network parameters $\Theta$ implied in various linear and nonlinear operations in the networks is performed by a training process, in which these parameters are chosen so that the network reproduces the models of the training data through optimizing the following loss function

$$ \Theta = \operatorname{argmin}(\| (S, T)^i - \Psi(h^i, \Theta) \|), \quad (7) $$

where $(S, T)^i$ denotes a model and its corresponding hydraulic head measurement $h^i$ obtained from solving the forward problem.

Hydraulic parameter field reconstructions rely primarily on the training database, building this database is hence the most critical task in this deep learning approach. Since a large database often promotes the accuracy of the network outcome, the database contains a large number of models, which requires a computational effort. As a direct solution, but also as a black box, it is not involved in further computations of a sensitivity matrix, as is the case with deterministic algorithms. Once trained, however, the neural network inverts the hydraulic parameters immediately and straightforwardly.

## 2.3 Multi-task Architecture

In this work, we develop a convolutional neural network based on the SegNet architecture that can interpret the measurement of hydraulic head at the input into the mapping of storativity and transmissivity fields at the output. SegNet is a autoencoder neural network that has been proven to effectively capture the features embedded in live video under real-world and complex circumstances. Originally, SegNet was designed for identifying trafict object, but here the architecture is redimensioned and adapted into a network that maps two fields of hydraulic properties from temporal images of hydraulic heads captured during pumping tests. The proposed multi-task architecture then consists of two SegNet-based networks individually optimized for each task (Figure 1). Each network consists of three encoder-decoder pairs, entailed with a regression layer at the output to map the corresponding field.
<figure><img src="image_2.png" /><figcaption>Figure 1: Multi-task architecture in X-shape composing two encoder-decoder network: one network (upper) detects the transmissivity heterogeneity connected to another network (lower) to map the tomography of heterogeneous storativity in matrix. The convolutional architecture consists of 95 layers with 0.6 million learnable parameters with multiple linking and bridging mechanisms.</figcaption></figure>

In each pair, the encoder extracts and condenses the key features found in the input images, and then transfers them to the next encoder while linking the indices to its decoder. This decoder, in turn, extracts the features from the previous layer and compensates for the loss of resolution due to the condensation process in the encoders by using the indices to prepare the images to have the same resolution as the input images. With the adoption of the gateway scheme, the accuracy of the resulting reconstruction maps can be significantly improved, which is an asset of this type of architectures. This SegNet network shares some topological similarities in its structure with other common encoder-decoder networks such as VGG16, UNet. However, by removing the fully connected layers from VGG16, SegNet encoders are significantly smaller and easier to train in practice. As the transfer of max-pooling indices instead of
values between encoders and decoders incurs only a fraction of the cost compared to UNet (Siam, et al., 2018). Thus, SegNet consists of a much smaller set of trainable parameters than competing architectures and can be determined using standard stochastic gradient descent. In fact, the total number of trainable parameters in this SegNet-based multi-task architecture is 0.6M, which is significantly smaller than other recent architectures proposed in the literature (Gao, et al., 2019; Misra, et al., 2016). Therefore, forming this small network requires less time and memory, making it an ideal candidate for multitasking.

Before building a multitasking network, it is necessary to closely examine the structure of the single-tasking network to understand how it works and to verify its limitations and advantages. In SegNet, each encoder and its associated decoder are arranged symmetrically to keep the input and output maps the same size. The encoding task is based on three operations such as convolution, Batch normalization and Relu performed sequentially on the encoders. The convolution calculation consists of scanning the input image with small filters to extract the features formed by the neighboring pixels. The resulting maps are then passed and filtered through a batch normalization (BN) and a rectified linear unit (ReLU), before being reduced by a Max-pooling operation. The encoders retain only the most important input information while reducing the dimensionality of the problem, thereby lowering computational costs.

The condensation process at the encoder stage builds low-level features and reduces the size of the parameters to be processed to better prevent the overfitting problem. On the other hand, the decoders unfold the condensed feature maps using an up-sampling layer formed with transfer of pooling indices from the encoders. Then the resulting maps are processed by CN, BN, and ReLU to successively translate the features and resize the output to the same resolution as the input. Understanding the functionality of each part of the network helps in designing an efficient multitasking architecture by choosing the optimal position to ensure information transfer between tasks in the connected network. The interest of multi-task learning is the establishment of links between different networks performing different tasks during simultaneous learning. Each task can benefit from the others by sharing information. Mutual learning serves as a complement of data for the others, reducing the risk of overfitting and achieving better generalization.
Moreover, both tasks are formalized with a single loss function that summarizes the objective function of each task. In this case each task sees the other as a regularization term (Ruder, 2017). The regularization term acts as a numerical stabilizer in the optimization process that helps to better determine the parameters of the network. However, one of the main challenges in designing an efficient multi-task architecture is determining how to share the features between tasks. Sharing at different layers often leads to a different feature distribution that can alter the efficiency of one or even all tasks (Misra, et al., 2016). Multi-task architectures usually have a shared trunk between tasks, where all tasks have shared layers and task-specific layers. The first one retrieves the generic features in the shared part of the tasks, while the latter handles the task-specific features. For optimal performance, it is recommended to find a compromise between these two parts in an architecture. In this aspect, a multi-task encoder-decoder structure with a shared encoder and a task-specific decoder seems promising. While the encoder focuses on low-level, cross-task feature learning, the decoders perform high-level feature interpretation to better learn task-specific details.

In practice, multi-task architectures with an encoder-decoder structure have been used in image processing (Kuga, et al., 2017) where the output variables share certain characteristics. However, this is not the case in our application where the two fields (transmissivity and storativity) are generated independently. For this reason, information transfer was only established at the decoder level by sharing low-level features and the decoder processes very distinct features. On the other hand, if we assume that the two fields in the output have common characteristics, the transfer can be performed both at the encoder and decoder (Kokkinos, 2017).

To connect the trunks, Misra et al. (2016) proposed the use of trainable scalars to perform the weighted sum of the features of the two tasks at multiple feature levels by constructing the cross-point units. These units provide shared representations between the two tasks, where their placements are decided in functioning of the dependencies between the tasks. Instead of trainable scalars, Gao et al. (2019) suggested a fusible style scheme in which the feature maps of multiple tasks are concatenated, and their dimensionalities are reduced by 1×1 convolutional layers. Following their work, we use 1×1
convolutional layers to link hierarchical features, in which the result of one task is added to another. Such additional layers are formed with learnable weights, that are more flexible and require less memory than concatenation. We retain the pooling index sharing mechanism in SegNet, which is computationally efficient and improves the accuracy of task-specific reconstructions. Thus, the new structure combines the linking and bridging mechanisms, reduces the number of learned parameters to 0.6 million that can be trained efficiently from a standard gradient descent scheme.

The training phase aims at estimating the network parameters by optimizing the following objective function, which is formulated as a sum of two objective functions for the two tasks.

$$L = L_1 + L_2, \qquad (10)$$

$$L_1 = \sum_{j=0}^{H-1} \sum_{k=0}^{W-1} \|m_T(j, k) - \hat{m}_T(j, k)\|, \qquad (8)$$

and

$$L_2 = \sum_{j=0}^{H-1} \sum_{k=0}^{W-1} \|m_S(j, k) - \hat{m}_S(j, k)\|, \qquad (9)$$

where $m_T = \log_{10}T$ and $\hat{m}_T = \log_{10}\hat{T}$, with T and $\hat{T}$ denote for the true and predicted transmissivity; H and W are the height and width of the output map. $m_S = \log_{10}S$ and $\hat{m}_S = \log_{10}\hat{S}$, with S and $\hat{S}$ denote for the true and predicted storativity.

## 2.4 Hydrogeological model

To construct a training database, we randomly generate the transmissivity and storativity fields from geostatistical variograms. Using SGEMS in MATLAB (Remy et al, 2009), two generations consist of 43,000 aquifer models, where the distributions of $\log_{10}T$ and $\log_{10}S$ are independent and regulated by Gaussian variograms whose mean, range, and sill parameters are randomly chosen with fixed intervals. In detail, both fields obey a variogram function $\gamma(u) = s[1 - \exp(-u/r)]$, but each field is generated independently, with the sill $s$ ranging from 0.1 to 1 and the range $r$ varying from 10 to 30 m. For the transmissivity field, its mean value is randomly selected at each generation from $10^{-3}$ to $10^{-5}$ m$^2$/s, with
transmissivity magnitudes ranging from $10^{-1}$ to $10^{-7}$ m$^2$/s. Similarly, the storativity mean ranges from $10^{-4}$ to $10^{-6}$ while storativity amplitudes range from $10^{-2}$ to $10^{-8}$. Regarding to the degree of the heterogeneity, both fields vary from 1 to 5 orders of magnitude. With this method, a database is created with different degrees of heterogeneity in a favorable generalization. In this study, no prior relationship is taken between transmissivity and storativity in the aquifers.

The parameterized fields of 32×32 cells are assigned to a 50 m × 50 m aquifer, on which 49 observation wells are regularly installed to measure variations in hydraulic head associated with water injection and extraction over time (see Figure 2 and Table 1). In this study, four injection/pumping pairs are successively implemented with a flow rate of 10 l/s. Two wells on opposite sides of the domain boundaries help to capture better the effects of aquifer heterogeneity on hydraulic head data.

<figure><img src="image_3.png" /><figcaption>Figure 2: Study area with monitoring setup including 49 observation wells and 4 pairs of pumping/injection wells at the border.</figcaption></figure>

In this theoretical study, the hydraulic head data are numerically solved from the forward problem as described in Section 2. The training database thus consists of transmissivity (32×32×1) and storativity (32×32×1) fields as the outputs and their hydraulic head responses collected with 49 wells at 10 timesteps in 2 minutes. To construct an input of the same size as the output, we use the nearest-neighbor spatial interpolation algorithm to form piezometric variation maps with 32×32 cells for each timestep.
Ultimately, the input data for each model is organized as a multi-channel matrix with (32×32×40) as the dimension, where the depth 40 is derived from the number of injection/pumping tests multiplied by the number of observations over time.

<table><caption>Table 1: Parameters used in the synthetic model</caption><thead><tr><td>Description</td><td>Value</td></tr></thead><tbody><tr><td>Study area</td><td>50 m ×50 m</td></tr><tr><td>Number of pumping/injection tests</td><td>4</td></tr><tr><td>Number of observation wells</td><td>49</td></tr><tr><td>Test time</td><td>2 minutes per test</td></tr><tr><td>Transmissivity</td><td>10<sup>-1</sup> to 10<sup>-7</sup> m<sup>2</sup>/s</td></tr><tr><td>Storativity</td><td>10<sup>-2</sup> to 10<sup>-8</sup></td></tr></tbody></table>

## 3. Results and discussion

This section focuses on analyzing the results and some interferences that control the reconstruction quality with a multi-task network, such as the size of the training data, the level of noise in the hydraulic head data, and the characteristics of hydraulic property fields when generating the data for learning. We also devote this section to include comparisons of the reconstruction accuracy from the proposed multi-task network versus those from the simple networks and the classical Gauss-Newton method.

### 3.1 Learning process and results

To train the network, we randomly select 30,000 models with the corresponding hydraulic data and another 3000 models for validation (Case C0). The training process is performed in MATLAB using the ADAM optimization algorithm with a constant learning rate of 0.001 and a batch size of 128. The training process is manually interrupted at epoch 227 to avoid overfitting and lasts 21 hours on an
Nvidia Quadro K2200 GPU in a Dell Precision Tower 5810. The evolution of the objective function during the training is detailed in Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3: Evolution of the loss functions during the training process.</figcaption></figure>

To quantitatively validate the resulting maps with the trained network, we use the coefficient of determination and the root mean square deviation, which are defined as follows

$$R_T^2 = 1 - \frac{SSR_T}{SST_T}, \qquad (11a)$$

$$R_S^2 = 1 - \frac{SSR_S}{SST_S}, \qquad (11b)$$

$$RMSE_T = \sqrt{\frac{\sum_{1}^{n} (m_T^{true} - m_T^{pred})^2}{n}}, \qquad (12a)$$

$$RMSE_S = \sqrt{\frac{\sum_{1}^{n} (m_S^{true} - m_S^{pred})^2}{n}}, \qquad (12b)$$

where $SSR_T = \sum_{1}^{n} (m_T^{true} - m_T^{pred})^2$ and $SST_T = \sum_{1}^{n} (m_T^{true} - \bar{m}_T^{true})$, with $m_T^{true}$ and $m_T^{pred}$ are the true and predicted log₁₀T, while $\bar{m}_T^{true}$ is the mean of the true log₁₀T field; n denotes the number of pixel (32×32), similar definitions for the metrics related to the storativity field S.
To evaluate the effectiveness of the trained network, we perform the test for 10,000 models that were not used in the learning phase. In this test, the network takes only 0.006 seconds to perform an interpretation of the simultaneous transmissivity and storativity tomography. Reconstructions of the entire test dataset yield average RMSE and $R^2$ values for the transmissivity models of 0.22 and 0.85, respectively, while the corresponding values for the storativity models are 0.33 and 0.63. The network appears to reconstruct transmissivity more accurately than the storativity field. This difference is due to the fact that storativity has less influence on the flow field in the aquifer, especially at steady state during the pumping test. To summarize, the reconstruction accuracy for the test data set is presented in Table 2 while its histogram is collected in Figure 7.

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Representative examples demonstrating the effectiveness of the trained architecture in reconstructing tomographies from simple to complex features (Case C0). From top to bottom: the true predicted maps and their comparisons; from left to right: three examples with their transmissivity and storativity fields. The proposed network provides accurate reconstructions for both fields, with accuracy depending on the complexity of target fields.</figcaption></figure>
To examine the quality of reconstructed maps in detail, three representative examples with different degrees of heterogeneity are selected and shown in Figure 4. In which, the trained model maps transmissivity and storativity from the measured hydraulic heads, which are then compared with the real models to evaluate the reconstruction accuracy. Generally, the proposed architecture successfully reconstructs the main heterogeneities of transmissivity and storativity fields. The quality of the reconstruction also exposes a dependence on the complexity of the target maps. In the case of Ex1, where both fields have smooth distributions, the network effectively detects their features including structure and magnitude. In the case of Ex2, although the transmissivity is well reconstructed, the strong heterogeneity in the storativity field is only mapped with some smoothness, since only the overall pattern is detected. Similarly, Ex3 is characterized by the presence of small heterogeneities that make their reconstruction a challenge for the process. These high heterogeneities are thus captured with smoothness. Indeed, the issue is a common challenge for all inversion methods where the reconstruction of complex distributions depends on the density of the data, demanding a high resolution to enable capturing the fine details (Jardani, et al., 2022; Ren & Kalscheuer, 2020; Liu, et al., 2018).

## 3.2 Effect of multi-task learning

To further evaluate the effectiveness of multi-task learning, we compare the results of the trained multi-task architecture with single-task networks. While a single-task network translates hydraulic head data into transmissivity and storativity maps, but by two separate networks (Case C1). As described in Section 2.3 for a multi-tasking design, the multi-task architecture is achieved by adding the interconnections between two single-task networks. Thus, the difference between the single-task and multi-task networks is only in the connection layers that are absent in the single-task network. We train the single-task network using the same hyperparameters, and data set as for the multi-task architecture. The trained networks are then tested on the same test dataset of 10,000 unseen models, whose overall metrics are summarized in Table 2 and Figure 9.
<figure><img src="image_6.png" /><figcaption>Figure 5: Reconstruction of transmissivity and storativity tomographies from two single-task networks (Case C1) for the same examples from Figure 4. The single-task networks are smaller with fewer learnable parameters but providing lower reconstruction accuracy.</figcaption></figure>

Validation of the single-task network on all test models yields an average RMSE and $R^2$ for the transmissivity reconstructions of 0.25 and 0.81, while the corresponding storativity reconstructions are 0.37 and 0.56. The results of single-task networks are indeed encouraging and competitive, as they are smaller, thus trained more quickly with fewer computational resources. However, the multitask networks offers more accuracy, with a 10% higher $R^2$ and 10% lower RMSE for both reconstructed tomograms. This leverage is confirmed by analyzing of the resulting reconstructions of representative examples performed by the single-task networks in Figure 5, in which the reconstructions for the same representative samples are made and compared with the truth models. The reconstructions for both fields with the single-task networks are relatively less accurate than those obtained with the multi-task network. Indeed, the reconstruction of Ex1 with the single-tasks is satisfactory, but the advantage always comes with outcome maps from the multi-task architecture.

The advantage is more highlighted when applying the multi-task network to complex heterogeneities that appear in Ex2 and Ex3. The trained multi-task excels in capturing these heterogeneities compared
to the single-task networks, where the reconstructions were found to be less accurate. As explained in Section 2.3, in multi-task learning, each task augments the input and regularizes the unified loss function, making the network more effective to capture details embedded in a complex environment. A similar conclusion has been drawn in the literature with applications of multitask networks in image processing (Horwath et al., 2020; Koziarski & Cyganek, 2018; Shorten & Khoshgoftaar, 2019).

### 3.3 Effect of number of observation points

Regardless of the optimization algorithms used, the reconstruction accuracy in an inverse problem is highly dependent on the quantity and quality of the measurements employed in this inversion. In this section, we examine how the quantity of hydraulic head data affects the reconstruction accuracy with the multi-task network. To this end, we compare the reconstructions from hydraulic head acquisition in three different configurations of observation wells. The first configuration represents a dense coverage with 49 wells, from which the testing reconstructions are presented beforehand in Figure 4, the second with only half of the original configuration (25 wells in Figure 6) and the last with 16 wells. The hydraulic head data retrieved with the new arrangements of observation wells are employed in the interpolation process to create input images for training. For each configuration, the multi-task network is then trained with corresponding altered databases to optimize the network performance and the resulting reconstructions for the test models are presented in the Table 2 and Figure 9.
<figure><img src="image_7.png" /><figcaption>Figure 6: Layout of a coarser monitoring scheme consisting of only 25 observation wells of unequal density, maintaining the highest coverage in the lower right corner.</figcaption></figure>

Overall, the number of observation points determines the reconstruction accuracy. In fact, the reduction of piezometric data leads to a degradation of the information described in the hydraulic level maps, which in turn leads to a significant reduction of the network effectiveness in predicting transmissivity and storativity fields. This effect is even clearly visible, especially when the aquifer properties are highly heterogeneous, as shown in Figure 7 for the representative samples where the reconstructions are performed using the coarse dataset and compared with the true ground models. For coarse features in Ex1, the neural network can reliably grasp the key features in both maps with some smoothing in the range of values. However, the trained network fails to reconstruct the finer details in Ex2 and Ex3 where the degree of heterogeneity in both fields is substantial. Differences in piezometric coverage over the study domain also affect the accuracy of reconstructions in Figure 7ghi, where the area better covered by more wells is better mapped. This impact of data coverage is also confirmed by the results for 16 observation wells, as shown in Table 2, with an apparent dependence of reconstruction quality on the input data density.
<figure><img src="image_8.png" /><figcaption>Figure 7: Reconstructions of tomographies with a monitoring scheme of 25 observation wells as a half of the original design (Case C2). The results refer to the same representative models in Figure 4. Low coverage of the monitoring data leads to a low accuracy of the predicted maps when only major features are determined.</figcaption></figure>

As mentioned before, it is well known that the effectiveness of an inversion algorithm in imaging heterogeneities in the subsurface depends on the ability of the data to cover these heterogeneities. More often than not, data sparsity makes reconstruction difficult and results in solution uncertainty (Ren & Kalscheuer, 2020; Liu, et al., 2017; Tarantola, 2005). In this test, we attempt to determine two tomograms with 2024 coordinates (2×32×32) from only 100 measurement points (4×25), the problem is thus highly underdetermined with a relatively low input/output data rate of about 5%. Although multi-interval temporal measurement can improve the data coverage, which are sometimes redundant and cannot completely solve the space constraint issue.
439

## 3.4 Effect of observation uncertainty

440

441

442

443

444

Noise is an inherent component of any field measurement and assessing the noise effect on the inversion result is a crucial element in determining model uncertainty. For this analysis, we contaminate the hydraulic head data with variable order Gaussian random noise with standard deviations of 5%, 10%, and 25% of the hydraulic head amplitude. The trained neural network then evaluates 10,000 sets of contaminated data, whose average metrics are shown in Table 2 and Figure 9.

445

446

447

448

449

450

451

452

453

454

455

456

457

458

The results from input data with different levels of noise show a clear impact of noise on the output reconstructions of the neural network, with the mapping quality deteriorating as the noise level increases. The impact is initially confirmed by the representative models with a noise level of 5% in Figure 8, which remains the likely noise amplitude for hydraulic measurements. Here the reconstructions are evaluated with contaminated measurements and compared with the truth model of both fields. The impact extent then appears limited when the features of the transmissivity and storativity fields are correctly reconstructed despite the noise. Mapping bias due to noise is then not apparent when comparing the reconstructed maps with the uncontaminated data in Figure 4 and the contaminated data in Figure 8. This noise resistance is related to the working mechanism of convolutional layers in the proposed architecture, where the processing of hydraulic data in these layers is operated by filters that perform region-wise sampling over the entire input maps. Such information extraction on the basis of neighboring points minimizes the potential abrupt noise at certain points and thus reduces the uncertainty of interpretation, which has been confirmed in other work (Jardani, et al., 2022; Vu & Jardani, 2021; Wu & Lin, 2018).
<figure><img src="image_9.png" /><figcaption>Figure 8: Tomography reconstructed from data contaminated with 5% noise (Case C3). The results refer to the same representative models in Figure 4. The reconstructions appear to be resistant to noise, with unclear degradation compared to the uncontaminated case in Figure 4. The advantage lies in the operation of the convolutional neural network used, where the translation is based on featuring of the maps.</figcaption></figure>

However, the mechanism does not seem to work effectively in the presence of high noise, since the output accuracy drops significantly when the noise level exceeds 10% and 25%. On the one hand, the noise in these situations is sufficiently high to deform the input maps completely, which is beyond the ability of convolutional layers to limit the damage (see Table 2 and Figure 9). On the other hand, this random noise also degrades the time evolution of the hydraulic signal and eventually accumulates errors in the input maps. In short, the network exhibits some resilience to a certain range of noise, although excessive noise can lead to significant degradation of the resulting maps.

To summarize, the accuracy metrics for the test dataset from the previous discussions are grouped in Figure 9 for histograms and in Table 2 for average ranges.
<figure><img src="image_10.png" /><figcaption>Figure 9: Histogram of reconstruction accuracy for 10,000 testing models in the discussion above. From top to bottom by row: Cases from C0 to C3; from left to right by column: Context of cases and corresponding histograms for transmissivity and storativity. Impacts of the interferences differ from case to case: while the outcome maps are most sensitive to the data scarcity, they are resistant to the measurement noise.</figcaption></figure>

<table><caption>Table 2: Summary of prediction accuracy tabulating the mean values of comparison between the truth and predictions for 10,000 testing models, performed by the neural networks from the discussions above.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Data</th><th rowspan="2">Case</th><th colspan="2">Transmissivity</th><th colspan="2">Storativity</th></tr><tr><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>Multi-task</td><td>Uncontaminated data from 49 observation wells</td><td>C0</td><td>0.22</td><td>0.85</td><td>0.33</td><td>0.63</td></tr><tr><td>Single-task results</td><td>Uncontaminated data from 49 observation wells</td><td>C1</td><td>0.24</td><td>0.81</td><td>0.37</td><td>0.56</td></tr></tbody></table>
<table><tr><td rowspan="2">Lower resolution dataset</td><td>Data from 25 observation wells</td><td rowspan="2">C2</td><td>0.28</td><td>0.75</td><td>0.39</td><td>0.52</td></tr><tr><td>Data from 16 observation wells</td><td>0.30</td><td>0.72</td><td>0.41</td><td>0.47</td></tr><tr><td rowspan="3">Contaminated data with a Gaussian noise</td><td>5% noise</td><td rowspan="3">C3</td><td>0.23</td><td>0.83</td><td>0.44</td><td>0.22</td></tr><tr><td>10% noise</td><td>0.27</td><td>0.76</td><td>0.54</td><td>-0.22</td></tr><tr><td>25% noise</td><td>0.48</td><td>0.17</td><td>0.71</td><td>-1.04</td></tr></table>

483

484

## 3.5 Influence of features in training models

485

486

487

488

489

490

491

492

493

494

The algorithm used in this paper is primarily based on a learning process that employs a database synthetically constructed. As a priori information, some choices on the nature of features in the hydraulic property fields have been made to reflect the known conditions of the water system under study, e.g., following Gaussian distribution in this study. Indeed, we reintroduce an essential element of an inverse problem - the a priori model, which is generally used to constrain the inversion solving, is critical to find a realistic solution and to avoid the problem of non-unity. Therefore, when analyzing the performance of a deep learning algorithm, evaluating significance of the a priori model is of interest to identify its applicability limit and generalization degree. To get a better idea about the limit and effectiveness of the proposed approach, in this part the trained network in Case C0 is recalled interpreting three test cases, of which the features differ from those predefined in the training database.

495

496

497

498

499

500

501

In the first case D1, transmissivity and storativity have a binary distribution (see Figure 10). The trained network is able to provide an accurate representation of the transmissivity field. However, in the case of the storativity map, some minor artifacts appear in the reconstruction. As the discrepancies between the regions of the two tomograms are somewhat less clear, this reflects the impact of the nature of the trained models, which are Gaussian with continuous distributions that disagree with the target maps. In the second case, the model also exhibits binary heterogeneity with a geochannel crossing the domain diagonally. The prediction effectively reconstructs the location and shape of the target channel,
but the *a priori* information of the training model also impacts the reconstruction by the presence of some artifacts, which are then minor.

<figure><img src="image_11.png" /><figcaption>Figure 10: Use the trained network to interpret models with features that differ from those of the trained dataset. From top to bottom by rows: the true, predicted maps and their comparison; from left to right by columns: three models from D1 to D3. The trained network succeeds in reconstructing the test models, with minor deformations caused by the constraints imposed by the given training dataset.</figcaption></figure>

In the last test D3, the network aims to detect a circular anomaly in both transmissivity and storativity. Two small anomalies are located in different places in the transmissivity and storativity models. The trained network succeeds in detecting the anomalies with some tolerance in terms of shape and values. Finally, although this trained network learns with generated patterns by geostatistical models, the proposed architecture has proven its ability to reconstruct various patterns that deviate from the learned features. The extensive generalization of the proposed network roots from an abundance of heterogeneity in the learning models which are generated from variograms with widely varying parameters. This property is indeed a key feature of any data-driven model that also works for a single-task architecture, as described in the literature, e.g., Vu & Jardani, 2021 & 2022. However, in any case,
the selection of learning models must be based on a priori information from other hydro-geo-physical investigations on the studied area. A favorable generalization of this algorithm yet offers a possibility to train from a database with a wide range of features, which should be helpful to handle uncertain in the a priori information (Jardani, et al., 2022).

## 4. Comparison with conventional Gauss-Newton algorithm

In order to test the effectiveness of the new method against a conventional inversion method, in this section we perform a simultaneous inversion of the transmissivity and storativity fields using the Gauss-Newton algorithm. Both fields are assumed to be independent and to follow the predefined Gaussian distributions with a variogram used as an a priori constraint of the database. The algorithm is formulated to minimize the following objective function:

$$L = [f(m_T, m_S) - h_{obs}]^T [f(m_T, m_S) - h_{obs}] + \alpha_1 (m_T - \bar{m}_T)^T G^{-1} (m_T - \bar{m}_T) + \alpha_2 (m_S - \bar{m}_S)^T M^{-1} (m_S - \bar{m}_S), \quad (13)$$

where $m_T$ and $m_S$ are parameter vectors of logarithm transmissivity and storativity respectively, while their corresponding prior values $\bar{m}_T$ and $\bar{m}_S$ are the mean of fields; $f$ is the forward problem and $h_{obs}$ are observed hydraulic head data; $G$ and $M$ are the parameter covariance matrices of transmissivity and storativity respectively; $\alpha_1$ and $\alpha_2$ are the regularization coefficients, in this test, get a value of 0.05.

Thus, the parameter fields $m_T$ and $m_S$ are determined on the basis of an iterative minimization of the function (13) in which both fields are proposed, tested and modified at each step until the convergence state. For the transmissivity

$$m_T^{k+1} = m_T^k + \beta_1 \Delta m_T, \quad (14a)$$

$$\text{with: } \Delta m_T = [J_T^T J_T + \alpha_1 G^{-1}]^{-1} J_T^T [h_{obs} - f(m_T, m_S)], \quad (14b)$$

$$J_T = \frac{\partial f}{\partial m_T}, \quad (14c)$$

and for the storativity field
$$ \mathbf{m}_S^{k+1} = \mathbf{m}_S^k + \beta_2 \Delta \mathbf{m}_S, \qquad (15a) $$

$$ \text{with: } \Delta \mathbf{m}_S = [\mathbf{J}_S^T \mathbf{J}_S + \alpha_2 \mathbf{M}^{-1}]^{-1} \mathbf{J}_S^T [\mathbf{h}_{\text{obs}} - f(\mathbf{m}_T, \mathbf{m}_S)], \qquad (15b) $$

$$ \mathbf{J}_S = \frac{\partial f}{\partial \mathbf{m}_S}, \qquad (15c) $$

where $\beta_1, \beta_2$ are the scalars determined by a line search at each step to speed up the minimization process; $\mathbf{J}_T$ and $\mathbf{J}_S$ are Jacobian matrices for the transmissivity and storativity fields respectively, which are computed with classical finite difference method, however, can also be derived using the adjoint state method to reduce computation time.

This inversion algorithm is applied to solve three representative examples in Figure 4, where the prior information for both fields recall the mean and covariance matrix used in generating the true distributions. The measurements used for the inversion are the same as those for the neural network described above, from the monitoring system with four pumping/injection tests and 49 observation wells. The solver employs the same computer as for training the proposed neural network, with a multicore CPU running on 12 processors. Each inference takes about 2 days and the resulting inversion for three cases are illustrated in Figure 11. Overall, the conventional inversion method yields comparable results to the neural network, although a detailed analysis reveals some differences between the two. The reconstruction quality relies on the selection of variogram parameters in the regularization conditions, where changing in the selection can reduce the mapping accuracy to a half. However, stepwise processing in a conventional approach provides an opportunity to intervene during optimization to adjust or improve the final solution with various considerations, which differs from an end-to-end deep learning operator. Further differences have been also described extensively in Jardani et al. (2022) between the Gauss-Newton and the deep learning algorithm with a CNN concept.

In terms of computation time, computing Jacobi matrices using the difference method typically requires a massive and tedious repetition of the previous problem solution. For example, in our scenario, each Jacobi computation requires the use of 32×32 direct solvers at each solution step. In addition, each inversion process often involves dozens of iterations, which is equivalent to generating the entire training data set. In contrast, Jacobian matrices can be solved within minutes using the adjoint state method, but it is a much more complex procedure, especially in the time domain (Sun & Yeh, 1990).
With such approach, operating a Gauss-Newton algorithm becomes faster. Other equivalent versions of Gauss-Newton, such as the Principal Component Geostatistical Approach (PCGA) in which the Jacobian matrix is introduced as a multiplication by a vector, significantly reduce the optimization time (Kitanidis & Lee, 2014). Thus, inference in a neural network can be superior in terms of time and computational cost only when compared to the classical Gauss-Newton method using a finite difference approach on a highly discretized domain. In general, the Gauss-Newton algorithm can be executed in a reasonable amount of time compared to a statistical learning algorithm where the construction of a training database remains a cumbersome task.

<figure><img src="image_12.png" /><figcaption>Figure 11: Reconstructions of tomographies from a conventional inversion method using the Gauss-Newton algorithm for the same representative models in Figure 4. The inverted reconstructions are of comparable accuracy, but require intensive computation and time for each inference, and their accuracy relies on the calibration of preconditions and the choice of initial models.</figcaption></figure>

Gauss-Newton is a deterministic algorithm in which the resulting solution is a local minimum that depends on the initial model with which the optimization begins. In general, the prior model is used as
the initial model; but if this model is uncertain, the computation may take a long time or even not lead to convergence. Deep learning algorithms are otherwise considered as global optimizers and their solutions hence do not depend on a given initial model (Hornik, et al., 1989). However, both algorithms share a dependence on the *a priori* model. While the conventional model is formalized in an objective function to be optimized with a deterministic algorithm, the deep learning algorithm relies on featuring the database. Therefore, deterministic algorithms must have a differential form, which is not the case with deep learning algorithms as they offer more possibilities for *a priori* models, even with a complex form (Jardani, et al., 2022). In this case, for example, a wide range of variograms is used to generate a database with a high degree of pattern diversity in order to achieve generalization of the trained network. Such a strategy allows to handle the uncertainty of *a priori* model in the DL algorithms. However, this uncertainty is complex to include in a deterministic algorithm which is often tied with a fixed *a priori* model, where the uncertainty in turn affects the quality of reconstructions. As shown in Figure 11 for the result obtained with a Gauss-Newton algorithm, the reconstruction quality varies from one case to another, depending on how close the *a priori* conditions are to the truth.

## Conclusion

In this paper, we introduce a deep learning approach to map transmissivity and storativity fields from hydraulic heads measured during pumping/injection tests. This mapping directly approximates the inversion operator using a neural network with a multi-task structure. To perform multiple tasks, the proposed network interconnects two single-task convolutional neural networks with encoder-decoder structure based on SegNet architecture, which maps two tomograms at the outcome. The key technique for building the multi-task network is to share hierarchical features between task processings that enhances the resulting transmissivity and storativity reconstructions. Even the two fields, transmissivity and storativity, may be independent, their information is both embedded in the measured hydraulic heads, with a complex interaction between the two in driving the hydraulic dynamics. However, by exploiting the similarity and contrast between the shared information of the tasks, the multitask learning filters out the interaction, thus improves the reconstruction quality of both fields. The network
interconnection can also act as a data augmentation leading to an improvement of the reconstruction
generalization in each task. At the same time, this strategy preserves the specific information of each
field by separating the process in the decoder part of both networks. With minimal learnable parameters,
the proposed network is efficient in terms of memory and inference processing time. Since this learning
process adopts a single unified loss function for both reconstructions, yet it can be effectively trained
with a standard gradient descent technique.

This new network is trained with a large synthetic database where aquifer properties are randomly generated using geostatistical variograms to provide models with different degrees of heterogeneity that ensures generalization of the network. These transmissivity and storativity models are used in solving the groundwater flow equation to obtain numerical piezometric responses, which are the input for training the network with two outcome fields. The yielding reconstructions of the multi-task network on the test models proved the effectiveness of multitask learning, which outperforms the single-task neural networks with 10% higher accuracy for both tomograms. The performance of the trained architecture is then analyzed by considering different situations that may arise in real applications, such as the scarcity of piezometric data, the alteration by noises, and uncertainties related to the *a priori* model used to generate the training models.

In this analysis, the accuracy of neural network reconstructions shows a clear dependence on data
coverage. Insufficient coverage often leads to a smooth interpretation of the fields, as the small volume
of data does not capture the full complexity of the aquifer. Indeed, the dependence of inversion quality
on data resolution is a common problem in inversion which is independent of the optimization tool
used. The proposed architecture also exhibits some resilience to data noise since the process is based
on feature extraction from the input image, an important advantage of the convolutional neural network.
The effectiveness of deep learning inversion is then tested for cases where the binary spatial distribution
of the T and S fields deviates from the Gaussian distribution used to generate the database. As the
network correctly identifies the main features of both fields, the proposed approach proves its broad
generalization in detecting even some unlearned features. However, to achieve this generalization, a
large training database is required, which can be time consuming to build.

Concerning the comparison of our network with the classical deterministic Gauss-Newton method,
four main points can be drawn: The first is that the deep learning algorithm can be considered as a
global optimizer, of which the result is independent of the initial model that is not the case for Gauss-
Newton algorithms. Second, both methods share a strong dependence on the *a priori* model where a
complementary information is required with the piezometric data to assist the algorithm in finding
realistic models of the subsurface. As in most cases, measurement data alone are not sufficient to reveal
the complexity of environments under study. However, the way introducing this *a priori* model in a
deep learning algorithm is totally different from a Gauss-Newton method. In a Gauss-Newton
algorithm, the *a priori* information is introduced into the objective function in the form of derivable and
quadratic terms, whereas a deep learning approach includes *a priori* constraints when generating
database for the training process.

Without requiring an *a priori* model in derivable form, a deep learning algorithm can otherwise offer
an advantage of including uncertainties of *a priori* conditions into the inversion. When the construction
of the database is built with a wide range of features to cover a variety of *a priori* conditions, the trained
network thus can gain a better generalization. The third point is sensitivity to noise. A deep learning
algorithm minimizes the impact of data noise on the result maps thanks to the convolutional operation
which processes the data image on a region-wise basis, as opposed to deterministic algorithms based
on every data point. The last point concerns the computation time. Building a deep learning algorithm
is more computationally intensive than a Gauss-Newton method when the generation of database often
requires computing the forward problem for thousands of times. However, this becomes very interesting
when the Gauss-Newton code adopts a finite difference method to compute the Jacobian matrix for a
large number of unknowns, which would also be extremely heavy. These pros and cons of both methods
should be further validated in light of a real dataset, which should be the objective of future work.
Albawi, S., Mohammed, T. A. & Al-Zawi, S., 2017. *Understanding of a convolutional neural network*. 2017 International Conference on Engineering and Technology (ICET), 1-6.

Badrinarayanan, V., Kendall, A. & Cipolla, R., 2017. SegNet: A Deep Convolutional Encoder-Decoder Architecture for Image Segmentation. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 39(12), 2481-2495.

Bear, J., 1972. *Dynamics of Fluids in Porous Media*. New York: American Elsevier Publishing Company.

Berg, S. & Illman, W. A., 2015. Comparison of Hydraulic Tomography with Traditional Methods at a Highly Heterogeneous Site. *Groundwater*, 53(1), 71-89.

Castagna, M., Becker, M. W. & Bellin, A., 2011. Joint estimation of transmissivity and storativity in a bedrock fracture. Water Resources Research, 47(9).

Castagna, M. & Bellin, A., 2009. A Bayesian approach for inversion of hydraulic tomographic data. *Water Resources Research*, 45(4).

COMSOL Multiphysics® v. 5.6. www.comsol.com. COMSOL AB, Stockholm, Sweden.

Fitts, C. R., 2012. *Groundwater Science*. Second Edition, Academic Press.

Franssen, H.-J. H., Gómez-Hernández, J., Capilla, J. E. & Sahuquillo, A., 1999. Joint simulation of transmissivity and storativity fields conditional to steady-state and transient hydraulic head data. *Advances in Water Resources*, 23(1), 1-13.

Gao, Y. et al., 2019. NDDR-CNN: Layerwise Feature Fusing in Multi-Task CNNs by Neural Discriminative Dimensionality Reduction. *IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, 3200-3209.

Hornik, K., Stinchcombe, M. & White, H., 1989. Multilayer feedforward networks are universal approximators. *Neural Networks*, 2(5), 359-366.
Horwath, J. P., Zakharov, D. N., Mégret, R. & Stach, E. A., 2020. Understanding important features of deep learning models for segmentation of high-resolution transmission electron microscopy images. *npj Computational Materials*, 6(1), 108.

Jardani, A., Vu, T. M. & Fischer, P., 2022. Use of convolutional neural networks with encoder-decoder structure for predicting the inverse operator in hydraulic tomography. *Journal of Hydrology*, 604, 127233.

Kitanidis, P. K. & Lee, J., 2014. Principal Component Geostatistical Approach for large-dimensional inverse problems. *Water Resources Research*, 50(7), 5428-5443.

Kokkinos, I., 2017. UberNet: Training a Universal Convolutional Neural Network for Low-, Mid-, and High-Level Vision Using Diverse Datasets and Limited Memory. *30th IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*.

Koziarski, M. & Cyganek, B., 2018. Impact of Low Resolution on Image Recognition with Deep Neural Networks: An Experimental Study. *International Journal of Applied Mathematics and Computer Science*, 28, 735-744.

Kuga, R. et al., 2017. Multi-task Learning Using Multi-modal Encoder-Decoder Networks with Shared Skip Connections. *2017 IEEE International Conference on Computer Vision Workshop (ICCVW)*.

Laloy, E., Hérault, R., Jacques, D. & Linde, N., 2018. Training-Image Based Geostatistical Inversion Using a Spatial Generative Adversarial Neural Network. *Water Resources Research*, 54(1), 381-406.

Laloy, E. et al., 2017. Inversion using a new low-dimensional representation of complex binary geological media based on a deep neural network. *Advances in Water Resources*, 110, 387-405.

Liu, B. et al., 2020. Deep Learning Inversion of Electrical Resistivity Data. IEEE Transactions on Geoscience and Remote Sensing, 58(8), 5715-5728.

Liu, P., Qiu, X. & Huang, X., 2016. Recurrent Neural Network for Text Classification with Multi-Task Learning. *Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence*, 2873-2879.

Liu, Y., Farquharson, C., Yin, C. & Baranwal, V., 2018. Wavelet-based 3-D inversion for frequency-domain airborne EM data. *Geophysical Journal International*, 213, 1-15.
Liu, Y., Farquharson, C., Yin, C. & Baranwal, V. C., 2017. Wavelet-based 3-D inversion for frequency-domain airborne EM data. *Geophysical Journal International*, 213(1), 1-15.

Li, Y. et al., 2020. Multitask Learning for Super-Resolution of Seismic Velocity Model. *IEEE Transactions on Geoscience and Remote Sensing*, 1-12.

McLaughlin, D. & Townley, L. R., 1996. A Reassessment of the Groundwater Inverse Problem. *Water Resources Research*, 32(5), 1131-1161.

Misra, I., Shrivastava, A., Mulam, H. & Hebert, M., 2016. Cross-stitch Networks for Multi-task Learning. *IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, 3994-4003.

Moseley, B., Nissen-Meyer, T. & Markham, A., 2020. Deep learning for fast simulation of seismic waves in complex media. *Solid Earth*, 11(4), 1527-1549.

Mukhopadhyay, P. & Mallick, S., 2019. Bayesian deep learning for seismic facies classification and its uncertainty estimation. *SEG Technical Program Expanded Abstracts 2019*.

Pham, N., Fomel, S. & Dunlap, D., 2018. Automatic channel detection using deep learning. *SEG International Exposition and Annual Meeting*.

Remy, N., Boucher, A., Wu, J., 2009. *Applied Geostatistics with SGeMS: A User's Guide*. Cambridge University Press.

Ren, Z. & Kalscheuer, T., 2020. Uncertainty and Resolution Analysis of 2D and 3D Inversion Models Computed from Geophysical Electromagnetic Data. *Surveys in Geophysics*, 41(1), 47-112.

Ronneberger, O., Fischer, P. & Brox, T., 2015. U-Net: Convolutional Networks for Biomedical Image Segmentation. arXiv:1505.04597

Ruder, S., 2017. An Overview of Multi-Task Learning in Deep Neural Networks. arXiv:1706.05098v1

Schmidt, J., Marques, M.R.G., Botti, S. et al., 2019. Recent advances and applications of machine learning in solid-state materials science. *npj Computational Materials*, 5, 83.

Şen, Z., 2015. *Applications in Hydrogeology for Geoscientists*. Elsevier.
739

Shorten, C. & Khoshgoftaar, T. M., 2019. A survey on Image Data Augmentation for Deep Learning.
*Journal of Big Data*, 6(1), 60.

740

741

Siahsar, M. A. N. et al., 2017. Data-driven multitask sparse dictionary learning for noise attenuation of 3D seismic data. *Geophysics*, 82(6), 385-396.

742

743

Siam, M. et al., 2018. A Comparative Study of Real-Time Semantic Segmentation for Autonomous Driving. *2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops* (CVPRW), 700-710.

745

746

Sun, A. Y., 2018. Discovering State-Parameter Mappings in Subsurface Models Using Generative Adversarial Networks. *Geophysical Research Letters*, 45(20), 11137-11146.

747

748

Sun, N.-Z. & Yeh, W. W.-G., 1990. Coupled inverse problems in groundwater modeling: 1. Sensitivity analysis and parameter identification. *Water Resources Research*, 26(10), 2507-2525.

749

750

Tarantola, A., 2005. *Inverse Problem Theory and Methods for Model Parameter Estimation*. Society for Industrial and Applied Mathematics.

751

752

Vu, M. T. & Jardani, A., 2021. Convolutional neural networks with SegNet architecture applied to three-dimensional tomography of subsurface electrical resistivity: CNN-3D-ERT. *Geophysical Journal International*, 225(2), 1319-1331.

753

754

755

Vu, T. M. & Jardani, A., 2022. Mapping discrete fracture networks using inversion of hydraulic tomography data with Convolutional Neural Network: SegNet-Fracture. *Journal of Hydrology*, 609, 1277-1292.

756

757

758

Wani, M., Bhat, F., Afzal, S. & Khan, A., 2019. *Advances in Deep Learning*, Springer Singapore.

759

Wu, X. et al., 2019. Multitask learning for local seismic image processing: fault detection, structure-oriented smoothing with edge-preserving, and seismic normal estimation by using a single convolutional neural network. *Geophysical Journal International*, 219(3), 2097-2109.

760

761

762

Wu, Y. & Lin, Y., 2018. InversionNet: A Real-Time and Accurate Full Waveform Inversion with CNNs and continuous CRFs. *The Journal of the Acoustical Society of America*, 144(3), 1683.

763

764

Yang, F. & Ma, J., 2019. Deep-learning inversion: A next-generation seismic velocity model building method. *Geophysics*, 84(4), 583-599.

765
766

Yeh, W. W.-G., 1986. Review of Parameter Identification Procedures in Groundwater Hydrology: The Inverse Problem. *Water Resources Research*, 22(2), 95-108.

767

768

Zhang, Y. & Yang, Q., 2017. A Survey on Multi-Task Learning. *IEEE Transactions on Knowledge and Data Engineering*.

769

770

Zhang, Z. & Lin, Y., 2020. Data-Driven Seismic Waveform Inversion: A Study on the Robustness and Generalization. *IEEE Transactions on Geoscience and Remote Sensing*, 58(10), 6900-6913.

771

772

Zhao, T., 2018. Seismic facies classification using different deep convolutional neural networks. *SEG Technical Program Expanded Abstracts* 2018, 2046-2050.

773