Parameter Identiﬁcation by High-Resolution Inverse Numerical Model Based on LBM/CMA-ES: Application to Chalk Aquifer (North of France)

: The present paper proposes the numerical solution of an inverse problem in groundwater ﬂow (Darcy’s equation). This solution was achieved by combining a high-resolution new code HYSFLO-LBM (Hydrodynamic of Subsurface Flow by Lattice Boltzmann Method), based on LBM, to solve the direct problem, and the metaheuristic optimization algorithm CMA-ES ES (Covariance Matrix Adaptation-Evolution Strategy) to solve the optimization step. The integrated optimization algorithm which resulted from this combination, HYSFLO-LBM/CMA-ES, was applied to the hydrogeological experimental site of Beauvais (Northern France), instrumented by a set of sensors distributed over 20 hydrogeological wells. Hydrogeological parameters measured by the sensors are necessary to understand the aquifer functioning and to serve as input data for the identiﬁcation of the transmissivity ﬁeld by the HYSFLO-LBM/CMA-ES code. Results demonstrated an excellent concordance between the integrated optimization algorithm and hydrogeological applied methods (pumping test and magnetic resonance sounding). The spatial distribution of the transmissivity and hydraulic conductivity are related to the heterogeneous distribution of aquifer formations. The LBM and CMA-ES were chosen for their proven excellent performance and lesser cost, in terms of both money and time, unlike the geophysical survey and pumping test. The model can be used and developed as a decision support tool for integrated water resources management in the region.


Introduction
Chalk formations form the most important and considerable aquifer in the northern part of France as groundwater represents 97% of the water supply in the "Hauts-de-France" region-especially for drinking water, public consumption, agricultural and industrial activities, and supporting river flows. Hydraulic relations between groundwater chalk and rivers allow a favorable exploitation for industrial sectors by pumping the water surface resources.
Groundwater flow is governed by the estimation of hydrodynamic properties, namely the transmissivity, hydraulic conductivity, and water content values. The mapping of these parameters is very complex because: (i) the cost of experimental tests in hydrogeological wells is very high in terms of time and budget, especially in the hydrogeological context related to deeper wells; and (ii) the application of geophysical surveys, especially magnetic resonance soundings (MRS), is not obvious as the quality of results are influenced and perturbed by noise measurements in the field, which depend on urban areas and can makes measurements very difficult to acquire. On the other hand, identified transmissivity values (the important parameter in the pumping test interpretation) are rare in the North of Paris basin (Hauts-de-France region) and, if they exist, are in private scientific reports and not available in the public domain. Indeed, methods and techniques to estimate this hydrogeological parameter, such as geotechnical approaches which are based on the grain size, the intrinsic permeability, or on the analyses or pumping tests, show several limitations in terms of economics, spatial scale, time, and the heterogeneity of the hydrogeological formations. On the other hand, the pumping test, which generally requires a longer duration of two, three, or more days to produce good-quality results, presents several practical and technical problems, especially in urban sectors (e.g., disturbances to the local population, flooding of the water supply networks, and drawdown in the wells of the city). Therefore, the numerical approach is considered to be the better solution of the transmissivity estimation. Hence, it is necessary to identify other approaches to deriving suitable values in order to understand and to define the principal parameters which govern the heterogeneity of the groundwater flow and, especially, to identify the transmissivity by adopting indirect processes.
Numerous works have focused their theses on inverse modeling, which is an important tool in the management and planning of water resources as it simulates and characterizes hydraulic parameters to erratic properties in hydrogeological models [1][2][3][4][5]. In hydrogeology, inverse modeling is widely used for prospective pumping tests [6,7], to identify the permeability coefficient of rock mass [8], to manage irrigation networks [9], to manage water in the urban context [10], and to analyze the groundwater quality monitoring network [11]. In addition, several numerical models in hydrogeology have used genetic algorithms in order to find the optimal solution of inverse problems in many fields of hydrogeology and in particular to identify the field of transmissivity of aquifers [12].
Solutions of inverse problems are determined in two stages. The first is the resolution of the direct model to estimate the calculated simulated values of state variables necessary for the construction of the objective function (the error between the calculated and observed state values). The second step is the minimization of the objective function whose minimum is the desired identified value. The resolution of the direct problem is generally accomplished by classical discretization methods such as the finite difference method, the finite element method, and the finite volume method. Moreover, for the minimization phase, most of the classical optimization algorithms can be implemented to solve the inv erse problem, but have the disadvantage to convergence towards the local minimum. To avoid this drawback, several works have focused on the use of metaheuristic optimization algorithms such as genetic algorithms (GA) [13][14][15] or evolutionary strategy algorithms (ESA) [16][17][18].
For two decades, the Lattice Boltzmann method (LBM) has been considered a serious alternative to classical computational fluid dynamic (CFD) methods for solving Navier-Stokes equations in complex flow configurations. Since its introduction, LBM has been successfully used in all disciplines related to flows in general, including porous medium, multiphase flows, wave propagation, multiphysics flows, and many others. The foundations of LBM are rooted in statistical physics and, in particular, kinetic theory. Indeed, LBM reproduces the movement of particles virtually placed on a predefined grid (called a lattice) which collide with each other and then propagate on this grid. In addition, the bases of the kinetic theory, on which LMB is based, consist of solving the Boltzmann equation, which describes the spatio-temporal evolution of the distribution function with a source term which represents the collision operation. By its design, LBM is naturally parallelizable and particularly suited to implementations on GPU (Graphics Processing Unit) architecture. It is this advantage which results in a considerable saving of computing time, which continues to make LBM more and more attractive. There is an abundance of literature on all aspects (theoretical and practical) of LBM, but we will only cite a few books that present it in detail [19][20][21][22].
In this work, we present the solution of an inverse problem allowing the identification of the heterogeneous transmissivity field of an experimental site instrumented in 20 wells. To do this, we developed a new integrated optimization algorithm called HYSFLO-LBM/CMA-ES, in which the direct model is solved by the code HYSFLO-LBM and the metaheuristic optimization algorithm CMA-ES due to [23]. Both the LBM and CMA-ES are presented in Sections 2.2 and 2.3 of this paper. Finally, we underline that the main novelty of the integrated optimization algorithm (HYSFLO-LBM/CMA-ES) lies in the fact that this digital tool allows identification of the transmissivity (or diffusivity) tensor and its spatial heterogeneity. With this numerical tool, we do not use the concept of zonation 6(the computational domain is represented by a limited number of areas with constant transmissivity), which is not desirable for study areas with high spatial variability. This paper is organized as follows. Section 1 presents the introduction and the motivation behind the choices of LBM and the CMAE-ES algorithm as the basis of the HYSFLO-LBM/CMA-ES code. The second section is devoted to the mathematical description of the different phases implemented to solve the inverse problem. Application of the new integrated optimization algorithm HYSFLO-LBM/CMA-ES to a realistic case is described in Section 3. Finally, Section 4 is dedicated to the discussion.

Governing Equations
The flow of a fluid in a porous medium is formulated by Darcy's law. For a porous medium of negligible deformations crossed by a slow flow fluid, this law expresses that the infiltration velocity → q is proportional to the gradient of the hydraulic head ∇ϕ. It should be noted that Darcy [24] established this law for a flow of water in a vertical cylindrical column homogeneously filled with sand. Moreover, Darcy's law was also found to be applicable to heterogeneous porous media and non-uniform flow. If we denote by ϕ the hydraulic head and by = T the transmissivity tensor, Darcy's law is then expressed by:

T∇ϕ
(1) In d-dimension space, the transmissivity tensor is represented by a matrix t ij 1 ≤ i ≤ d Note that the system of Equation (1) is not closed (2 unknowns ( → q and ϕ) and 1 equation). To complete this system, we will use the continuity equation which assumes that the fluid is incompressible, or in other words, the divergence of the velocity is zero (∇. → q = 0). Finally, the association of the continuity equation and Darcy's law will give an elliptical partial differential equation (PDE) with unknown ϕ. For a given transmissivity field t ij 1 ≤ i ≤ d 1 ≤ j ≤ d , this PDE makes it possible to calculate the hydraulic head ϕ and the velocity of the flow to be deduced in post-processing using Darcy's law. Moreover, if the flow occupies the domain Ω Ω ⊂ R d of the boundary ∂Ω such as ∂Ω = Γ D ∪ Γ N and Γ D ∩ Γ N = ∅, the mathematical model describing the stationary flow in a porous medium is given by: where ϕ D is a known function for imposing the Dirichlet condition Γ D . ϕ N is also a known function for imposing the Neuman condition on the domain boundary Γ N . f is the source/sink term in the domain Ω. The vector → n denotes the vector normal to Γ N oriented towards the outside of the domain.
For a given transmissivity field, the mathematical model (2) describing the stationary flow in a porous medium is also called the direct problem, which we denote hereafter as (DP). Alternatively, if the hydraulic head ϕ is known at a few points in the Ω (of the measured values for example), the model (2) allows us, in this case, to determine the tensor of transmissivity T. In this case, the solution of the inverse problem is denoted by (IP).
Under certain regularity conditions, Ciarlet [25] theoretically showed that the mathematical model (2) admits a unique solution. For more details on the mathematical analysis of the problem (DP), the reader may also consult the book of Dautray and Lions [26].
To solve the direct problem (DP), we developed a numerical model based on the Lattice Boltzmann method. The foundations of this method are presented in the following subsection.

Lattice Boltzmann Method
The Lattice Boltzmann method (LBM) is a relatively new numerical method when compared with the classical approaches used in numerical simulation. It appeared in the 1990's and was initially deduced from the methods of gases on the network and cellular automata [27]. Unlike classical approaches based on the discretization of the Navier-Stokes equations (NSE), the Lattice Boltzmann method is based on the formalism of statistical physics, which consists of the numerical resolution of the Boltzmann equation. This equation is concerned not only with macroscopic quantities (speed, pressure, density), but directly with the distribution of the various particles constituting a fluid. This is called a mesoscopic representation. The simplicity and locality of the LBM method algorithm allows its easy and efficient parallelization. Thus, very quickly, this method was used for unsteady and incompressible CFD calculations [28]. Therefore, as long as the Mach number of the flow remains sufficiently low, LBM allows us, by its nature, to simulate the behavior of a fluid governed by unsteady, weakly compressible Navier-Stokes equations.
The Boltzmann equation (BE) plays a main role in the kinetic theory of gases. In the absence of external force, it expresses the convective transport equation of the velocity This distribution is none other than the probability of a particle to be found at the position where Ω denotes the rate of change of f due to the collision of particles (also called the collision operator). Ω can be approached by several simple models [29], but the most popular in LBM is the Bhatnagar-Gross-Krook (BGK) model (1954). BGK assumes that Ω is a function of the f distribution, the equilibrium distribution f (eq) , and the relaxation time τ (the time necessary to bring the system back to the state of equilibrium), as: The Equations (3), (4), and (6) constitute BE which must be solved numerically.
From the second principle of thermodynamics (or the application of the H-theorem), we can deduce that f (eq) is of Maxwellian type of the form: where d is the space dimension (d = 1,2,3), ρ and → v are, respectively, the macroscopic fluid density and velocity, R = k B /m with k B is the Boltzmann constant, m is the molecular mass, and T is the temperature of the system. Note that the artificial sound speed c s is defined as c 2 s = RT. Numerical evaluation on a computer is very costly in terms of CPU computing time. We therefore often look for an approximate expression, based on polynomials for example. If we assume that → v c s then the Taylor expansion of the exponential function allows us to approximate (5) by: To solve the BE, we start by defining a computational grid (called a lattice) which will be used for both spatial discretization To solve the BE, we start by defining a computational grid (called a lattice) which will be used for both spatial discretization ⃗ ⃗ and speed discretization ⃗ . To recover the NSE by the numerical resolution of BE, it is imperative that the discrete speeds must be chosen so that the following quadratic equation of the moment ℳ ( ) of order (k) is exact: where and ⃗ are the coefficients and points of the quadrature Equation (7). In 2D computations, LBM uses the Gauss-Hermite quadrature equation with at least five points for the numerical integration. The quadrature equation being chosen allows us to determine the coefficient and the direction ⃗ . With these two parameters, we can then discretize BE by setting (⃗ ⃗ , ) = (⃗ ⃗ , ⃗ , ), which must verify the BE in turn as: In 2D, there are two types of lattice: square and hexagonal. The most commonly used lattice, which leads to more accurate results, is the square lattice D2Q9 (calculation in 2D with a discretization of the velocity variable ⃗ in 9 directions, see Figure 1). If the D2Q9 lattice is adopted ω i , → c i , and f (eq) i are given as: where c is lattice velocity c = ∆x/∆t, with ∆x being the lattice size and ∆t being the time step.
where c s = c/ √ 3. Note that, at this stage, we only give the discretization according to the velocity space. For the spatio-temporal discretization, LBM uses the characteristic method but without interpolation, since the lattice should be chosen in such a way that ∆ → x = → c ∆t. Thus, the final discretization of the BE is given by: 8 (12) Once this equation is solved, the macroscopic density and the velocity of the fluid can be recovered by: In this paragraph, we wanted to give a general description of LBM, and it is for this reason that we treated it as a fluid flow governed by the macroscopic model of NSE. As the aim of this paper is to solve the Darcy equation by LBM, we have to give only the new expression of f (eq) (as mentioned above) which allows us to cover Darcy's macroscopic equation. The Darcy equation is a special case of the NSE for slow flows (low Reynolds number) in which the inertia term is considered negligible. Therefore, the stationary Darcy's equation is elliptical (flow dominated by diffusion processes) and resembles a diffusion equation. In a D2Q9 lattice, the equilibrium distribution function for a diffusion equation is given by:

Solution of Darcy Equation by LBM
The solution of Darcy's equation by LBM is done in two steps: collision and streaming.
If we consider that the transmissivity tensor streaming: where λ is the dimensionless relaxation time (λ = ∆t/τ). Developments in LBM show that the viscous effects of the Navier-Stokes equations are related to the relaxation time. It is the same for the transmissivity for the Darcy equation. Thus, the macroscopic transmissivity coefficient T(x, y) can be related to the relaxation factor λ as: Finally, the expression (15) is used to compute the equilibrium distribution function with the values of the weighting coefficients ω i given by expression (14). The hydraulic head ϕ can be calculated as: For the boundary conditions, LBM requires conditions on the different directions of the distribution function f i while the boundary values are given on the primitive variable of the macroscopic model. It is, then, necessary to transform the boundary values from the macroscopic model to the mesoscopic model. To answer this question, several transformations have been proposed in the literature [30,31]. These transformations allow to impose various types of boundary conditions (Dirichlet, Neuman, Robin, etc.) on the pressure as well as on the velocity of the flow. In this paper, we use the Dirichlet-type boundary conditions on the four boundaries of the studied domain by applying the flux conservation principle which is expressed by: where opp(i) denotes the direction opposite to the direction i. For instance, suppose we know f 7 and ϕ on the boundary denoted ϕ wall , then expression (19) allows us to transform the value of ϕ wall into f 5 by: Finally, we give the Algorithm 1 implemented in the computer code HYSFLO-LBM to solve the Darcy equation, which was used for the direct problem.  (14)) loop: t + ∆t compute λ ij from (Equation (17)) compute f (eq) k from (Equation (14)) collision from (Equation (15)) streaming from (Equation (16)) set boundary condition from (Equation (19)) compute ϕ from (Equation (18)) test ϕ − ϕ 0 < tol exit set ϕ 0 = ϕ go to the next time step

CMAES Algorithm
CMAES belongs to the category of evolution strategy algorithms. It has become the primary algorithm for free-gradient optimization. It is well suited to solving continuous optimization problems where the objective function is not known explicitly and is not necessarily convex. Although it is stochastic in nature, the mutation stage of the CMA-ES algorithm is considered to be correlated and deterministic. This property makes this algorithm easy to implement and less computationally intensive. Moreover, like metaheuristic optimization algorithms, CMA-ES has the advantage of converging towards a global optimum. CMA-ES is gaining popularity and is becoming the benchmark algorithm in metaheuristic optimization. It has been successfully applied to several engineering disciplines including: environmental engineering [32], acoustics [33], electronics [34], hydrogeology [35], medicine [36], thermal and fluid flow [37], structural mechanics and failure [38], and many others. CMA-ES is particularly efficient for non-convex, poorly conditioned, multimodal optimization problems and with noisy evaluations of the objective function.
The CMA-ES algorithm is from the family of strategy evolution algorithms and, like genetic algorithms, is inspired by the Darwinian theory of evolution. CMA-ES is performed in four steps: initialization, selection, recombination, and mutation. These steps operate on a set of µ parents to produce λ children, which we now denote by CMA-ES (λ, µ). In order to explore the search space, CMA-ES uses candidate populations according to a multivariate random distribution. The mutation step is considered to be the main operation in the CMA-ES algorithm. It allows us to produce a new population by adding a multivariate random vector to the parent population. The mutation also acts as the guide of CMA-ES for the different transformations (rotation and homothetic) of the adapted covariance matrix from the generated population. It should be noted that the mutation process is based on a set of parameters (called strategy parameters) which are updated automatically without user calibration, but by exploiting the various information from previous generations [16].
Like any iterative algorithm, CMA-ES begins with an initialization step for all of the algorithm's control parameters. This step also consists to initialize a population of λ individuals randomly according to the multivariate normal distribution. Hansen [39] suggests considering a population of size λ = 4 + 3ln(d), where d is the size of the optimization variable. It also indicates that this value is reasonable for large optimization problems. In the rest of this paper, we denote by (X (g) i=1,λ ∈ R d ) the individual's population of generation g, and by F the objective function.
After the initialization stage, CMA-ES enters the loop of generations until convergence towards the optimal solution. This loop begins by evaluating all of the individuals in the population by the objective function to obtain the fitness value F X CMA-ES then selects the "best" µ (µ ≤ λ) individual among the λ individuals based on their fitness value ("best" here designates the rank of the X (g) i according to its fitness value: "smaller" if this is a minimization problem or "bigger" if it is a maximization problem). Then, we evaluate the weighted average vector on these µ individuals noted m (g) and given by: where w k are the weighting coefficients given by Hanssen [39]. It is the expression (21) which will then contribute to the construction of the next generation of parents X (g+1) i=1,λ . The CMA-ES algorithm mutation step consists of adding to the population mean vector of the generation (g), a noisy component according to a multivariate normal distribution as: where σ is a positive parameter and N (0, C) denotes a vector of independent normal random numbers with zero mean and covariance matrix C (symmetric and positive definite). One can observe that for a good mutation which leads to a rapid convergence towards the optimum, it is necessary to choose judiciously the parameters σ and the matrix C. However, in practice this choice turns out to be difficult. It is exactly at this point where the power of the CMA-ES algorithm appears, which automatically adapts these two parameters during successive generations and without user intervention. Thus, the CMA-ES algorithm can be summarized in the following steps named Algorithm 2: Algorithm 2: CMA-ES Steps set numerical parameters (d: size of the optimization variable) initialization σ, m, C generation: g + 1 selection → X, m, µ mutation → X = m + σN (0, C) evaluation → F(X) test F(X) < tol exit adaptation → σ, C go to the next generation

The Integrated Optimization Algorithm (HYSFLO-LBM/CMA-ES)
In hydrogeology, the transmissivity tensor = T is one of the key parameters in the modeling of groundwater flows. Its estimation with precision is capital to deduce (by the Darcy equation) an accurate velocity field. A transmissivity field identified with good precision also guarantees precision of the diffusion coefficient of the studied porous medium, since the diffusion is deduced from the velocity field (and therefore from the transmissivity). Therefore, the identification with high accuracy of the tensor = T is also necessary for the prediction of pollutant propagation in a porous medium. It is therefore important to deploy the most sophisticated numerical tools to identify the transmissivity tensor = T. This is explained in the following paragraphs. If we have a series of observed values of the hydraulic potential (ϕ obs ), we can formulate the identification of = T by a mathematical optimization problem with constraints.
In other words, the problem is equivalent to the search for the tensor = T which minimizes the error between the observed values (ϕ obs ) and those computed by the direct problem ϕ comp . Recently, this problem was solved by a few researchers using manual regional adjustment of = T until an error was obtained which they considered reasonable! With the spectacular advances in the field of mathematical optimization, and especially in metaheuristic algorithms, the identification of the tensor = T can be treated as the resolution of an inverse problem (IP). In other words, identification of = T becomes automatic without any manual adjustment. In this case, this problem can be formulated as the minimization of an objective function F (which is the error between ϕ obs and ϕ comp ) with direct model (Equation (2)) itself as a constraint. Thus, the inverse problem to solve is: where F is the objective function defined by F Note that the problem (23) can be confronted with the problem of computational stability (sensitivity to measurement noise), but also with the problem of the uniqueness of the solution = T opt . This problem is called, the "ill-posed problem". In fact, in the groundwater flow, different hydrogeological conditions can provide identical observations of the hydraulic potential or solute concentration. It is then impossible to determine uniquely the transmissivity (or diffusivity) tensor only from observations. Consequently, problem (IP) most often requires a well-chosen regularization strategy to guarantee a certain stability of the result.
Several solutions have been proposed in the literature to remedy the ill-posed problem in groundwater flow modeling. We cite for instance: (i) minimization of non-linearity and non-convexity during the execution of the optimization algorithm by transforming the Because it was formulated within a rigorous mathematical framework, the strategy of regularization made it possible to prove both the existence, the uniqueness, and the stability of the solutions of the inverse problems under certain regulatory assumptions. Several regularization strategies have been proposed in the literature. In his book, ref. [40] analyzes all of these strategies by indicating their conditions of applicability. In this paper, we use the Tikhonov regularization method which consists of modifying the objective function by introducing a regularization term based on the a priori information. In this case, the new objective function which will be considered in problem (IP) is: where = T * are some known values of transmissivity in the study domain (in absence of = T * , this value must be zero) and alpha is the regularization parameter. The choice of α should be consistent with the inaccuracy of the input data. Ref. [41] proposed various practical choices which also prove to be useful for the efficiency (in the sense of convergence) of the algorithm implemented for solving the inverse problem (23).
) < tol then convergence to the = T opt and exit adaptation → σ (g+1) , C (g+1) go to the next generation

Realistic Case: Application to the Experimental Hydrogeological Site of Beauvais (Unconfined Aquifer)
This case is presented to demonstrate the successful application of the integrated LBM and CMAE-ES algorithm as the basis of the HYSFLO-LBM/CMA-ES code to solve the inverse problem in handling transport through fractured media, especially the chalk aquifer in the north of the Paris basin (Figure 2a).   [43] and designed to measure, especially, water pressure), which contains electrical conductivity and temperature sensors, memory for storing measurements, and a battery. Recorded data are afterwards stored in the DIVER's internal memory. Indeed, the DIVER consists of a pressure sensor, designed to measure water pressure, and a temperature sensor. The DIVER contains an autonomous datalogger and the duration of the measurements can be chosen by the user, ranging from seconds to hours. The various recorded measurements can be recovered either in the laboratory or at the measurement site using an optical communication system linking DIVER to a laptop computer.
Modeling of the Hydrogeological Experimental Site of Beauvais system requires expertise in terms of hydrodynamic characteristics, water table levels, the computed hydrological balance, transmissivity, and hydraulic conductivity. Generally, these parameters result from geostatistical concepts, laboratory experiments, pumping tests [44,45], and magnetic resonance sounding (MRS) for their hydrogeologic estimation [46]. Indeed, MRS investigations summarize the definition of the hydrodynamic parameters, especially in Beauvais city (in the North Paris Basin) and near the Hydrogeological Experimental Site, and will serve to ameliorate the database which should be used in the numerical processes of this paper and contribute to the knowledge of the heterogeneous permeable bodies [47] in the chalk aquifers.

Modeling Results
The integrated optimization algorithm HYSFLO-LBM/CMA-ES was applied to the instrumented zone described in the previous paragraph. It is a rectangular area of 123 m wide and 133 m long. For high computational resolution, the size of the square lattice ∆x was chosen to be 1 m. The study area has 20 measuring points of the hydraulic head ϕ. An interpolation procedure based on the Kriging techniques was performed to estimate the value of the hydraulic head on the entire computational grid (123 × 133 nodes) from the 20 measured points. These Kriged values were subsequently considered to be the field of observed values and denoted by ϕ obs . The integrated optimization algorithm aimed to identify the field of the aquifer transmissivity during the year 2016. Table 1 presents the precipitation values of this year from which we have deduced the values of the source term f (see Equation (2)) necessary for the direct model HYSFLO-LBM. Finally, the direct model also required the boundary conditions for its forcing. The values of the Dirichlet boundary condition ϕ D (see Equation (2)) at the four edges limiting the domain were calculated by evaluating ϕ obs on these four limits as ϕ D = ϕ obs (Γ D ). Moreover, the values of = T * necessary for the regularization of the inverse problem (see Equation (25)) were determined by using some values found in the literature from previous studies carried out on neighboring areas with the same hydrogeological characteristics.
The aquifer profits from principal refill, which is materialized by flow collected by the HESB. The lithological nature of the aquifer system is composed, especially, by the chalk deposits. It is very important to note that this formation is marked by dual porosity, with both interstices and cracks. The primary type of hydraulic conductivity linked to the interstitial porosity of the aquifer remains very low and generally does not exceed 10 −0.5 m/s [48]. The chalk fracturation at the surface permits higher hydraulic conductivities, between 10 −3 and 10 −2 m/s, which influence the runoff of the groundwater and create turbulent flows [49]. The comprehension of the infiltration and the identification of recharge and discharge periods is determined by the analysis and the interpretation of hydraulic head and rainfall evolution ( Table 1). The hydraulic gradient is about 0.0047 between the wells F1 and F4 and 0.0066 between F2 and F4. Based on hydro-geological knowledge, physical arguments, and piezometric analysis and interpretation, the hydrodynamic model is built in order to simulate the groundwater flow and to determine the distribution of the transmissivity values of the chalk aquifer in the HESB. The numerical process of achieving the objectives in this modeling expertise involves a series of procedures. They can be classified according to three principal steps: data collection; computer simulations, calibration, and analysis; and numerical results (Figure 4a-c). At the domain boundary and after hydrodynamic conditions, the piezometric values are imposed. Simulations are conducted in steady state with the reference values from 2016. The process of georeferenced nodes allows us to accelerate the modeling calibration (Figure 4b). This numerical method highlights an overlay between measured and computed hydraulic head, where the maximum recorded was 72.6 m and 71.92 m, respectively.   In order to stimulate the quality of the calibration, the absolute and the relative hydraulic head errors and are estimated using the following expression: Prior to interpreting the results, we would like to point out that we distinguish two areas: a diagonal band around the 20 measurement points and an area on either side of this diagonal. We recall that in the second zone it was necessary to complete the hydraulic head values by interpolation in order to have a reference values field to formulate the objective function. It is quite obvious then that the results obtained in this zone will be more precise than that of the measurement zone, and occasionally with errors much lower than the model precision. Consequently, the interpretation and the discussion of the obtained results that we will present in the rest of the paper should concern only the diagonal band around the points of measurements.
The result, in the form of an errors map, displays map results displays relative and absolute errors in order to locate the points where there is a very large difference between measurement and simulates hydraulic heads. Generally, we note three principal zones which showed the existence of relatively small margins of absolute and relative error (Figure 5a,b): the first margin of absolute error from: 8.0 × 10 −3 m to 2.0 × 10 −2 m which corresponds to margin of the relative error between 1.0 × 10 −2 % and 2.5 × 10 −2 %; -the second margin of absolute error from: 2.0 × 10 −2 m to 2.8 × 10 −2 m which corresponds to margin of the relative error between 2.5 × 10 −2 % and 3.5 × 10 −2 %; In order to stimulate the quality of the calibration, the absolute and the relative hydraulic head errors E a and E r are estimated using the following expression: Prior to interpreting the results, we would like to point out that we distinguish two areas: a diagonal band around the 20 measurement points and an area on either side of this diagonal. We recall that in the second zone it was necessary to complete the hydraulic head values by interpolation in order to have a reference values field ϕ obs to formulate the objective function. It is quite obvious then that the results obtained in this zone will be more precise than that of the measurement zone, and occasionally with errors much lower than the model precision. Consequently, the interpretation and the discussion of the obtained results that we will present in the rest of the paper should concern only the diagonal band around the points of measurements.
The result, in the form of an errors map, displays map results displays relative and absolute errors in order to locate the points where there is a very large difference between measurement and simulates hydraulic heads. Generally, we note three principal zones which showed the existence of relatively small margins of absolute and relative error (Figure 5a the third margin of absolute error from: 2.8 × 10 −2 m to 3.8 × 10 −2 m which corresponds to margin of the relative error between 3.5 × 10 −2 % and 5.0 × 10 −2 %. From these two figures we should retain that the absolute and relative errors do not exceed 4 cm and 5.0 × 10 −2 %, respectively. These two errors are located exactly at the hydrogeological well "F3" where the hydraulic mechanism is characterized by the groundwater dividing axis and a probable network of fractures, allowing control of the groundwater flow. Finally, in view of the relative errors (Figure 5a), and of the calibration figure (Figure 4b), the integrated optimization algorithm HYSFLO-LBM/CMA-ES allows us to obtain very good results both for the identification of the transmissivity and the simulation of the flow in the experimental site of Beauvais.
In the same way and for the same reasons cited above for the error fields (Figure 5a,b), we subsequently interpret the transmissivity and hydraulic conductivity fields only for the diagonal band, including the 20 measurement points.
The HYSFLO-LBM/CMA-ES code identifies the distribution of the hydraulic conductivity (or transmissivity) and its possible heterogeneity (Figure 6a,b). It is for this reason that the proposed integrated optimization algorithm punctually identifies (point by point) the tensor = T. This step was carried out by solving an optimization problem of the variable = T as a real vector of 123 × 133 components which corresponds to the total number of the computation grid nodes. It should be noted that the CMA-ES algorithm is a stochastic search algorithm affected by certain errors induced by the generation of individuals randomly. In order to reduce these errors, we carried out about ten simulations to consider, as the final value = T opt , only the average of these simulations. Consequently, Figure 6a,b shows, outside the diagonal band of the measurements, the heterogeneity of the hydraulic conductivity and of the tensor = T, but not of the random noises that could probably be interpreted by the reader. Figure 6b shows lower, middle, and higher transmissivity values are about 0.0025, 0.0045, and 0.007 m 2 /s, respectively, but globally, the magnitude order is about 10 −3 m 2 /s. This variation is very heterogeneous in the hydrogeological experimental site and is represented by the spatial distribution of transmissivity values (Figure 6a).
The same spatial distribution concerns the hydraulic conductivity. These values were obtained using simulated transmissivity values and the potential hydraulic head of the experimental site, and by considering the hydraulic character of the unconfined chalk aquifer. The hydraulic conductivity heterogeneity ranged from 5.0 × 10 −4 to 10 −4 m/s.

Discussion
The use of the Lattice Boltzmann method (LBM) in the numerical modeling of the groundwater allows to identify hydraulic characteristics of the chalk aquifer of the Hydrogeological Experimental Site of Beauvais, which revealed numerous transmissivity values ranging from 0.0025 to 0.007 m 2 /s. The spatial repartition of hydrogeological characteristics is accompanied by the anisotropic character of the chalk, which is characterized by the horizontal conductivity (10 −14 to 10 −12 m/s) and the vertical conductivity (10 −15 to 10 −13 m/s) [50]. The transmissivity and the hydraulic conductivity constitute principal characteristics allowing us to understand the groundwater flow. The variation of the transmissivity (0.0045 and 0.007 m 2 /s) and the hydraulic conductivity (5.0 × 10 −5 to 10 −4 m/s) is accompanied by the spatial distribution of the water content in the chalk aquifer, with maximum values of about 30-35% [46], and by the heterogeneity of lithological formations which are deduced from the analysis and the interpretation of drilling cuttings (dating from 2014 and with GPS coordinates: 49 • 27 36.88" N 2 • 04 17.38" E) belonging to the experimental hydrogeological site of Beauvais. It is mainly lithological deposits which are composed of clay with flint (Figure 7a), chalk with or without clay, and flint. This variation could be explained by primary and secondary hydraulic conductivities as defined by [50] and which are from 10 −8 to 10 −5 m/s and 10 −3 to 10 −1 m/s, respectively. The analysis and the interpretation of transmissivity values in different media were deduced from: the resonance magnetic sounding (Figure 7b,c); and the pumping test, especially in the 47 wells (chalk aquifer) and concerned plateau, dry valley, and humid valley with orders of magnitude about 10 −3 m 2 /s in plateau and an average transmissivity of about 6.7 × 10 −3 m 2 /s [51]. From the geological structures point of view, value ranges which are deduced from the HYSFLO-LBM/CMA-ES code could be explained by the presence of fracture networks. Indeed in 2011, we realized near the hydrogeological experimental site 424 of fractures measurements. Three orientations have been identified: N-S, WNW-ESE, and NNE-SSW. In fact, numerous works highlighted a positive correlation between fracture size and transmissivity [52] on the one hand and, on the other hand, the influence of the fracture network and the faulting mechanism on the spatial distribution of hydrogeological characteristics (transmissivity and hydraulic conductivity) and on the groundwater flow [45].
It is difficult to highlight in the field the relationship between fractures and the repartition of transmissivity values. Note the concordance transmissivity values resulting from hydrogeological experimentation (pumping tests and MRS) and the numerical modeling simulation by using the LBM. Hence, the accuracy of the results indicates that the use of the HYSFLO-LBM/CMA-ES code is recommended to simulate the flows governed by partial differential equations (PDEs) given by the model presented in this paper.
The integrated optimization algorithm HYSFLO-LBM/CMA-ES could be used in the industrial field (agricultural and food activities) as well as in the sector related to the supply of drinking water. The HYSFLO-LBM/CMA-ES code, as a numerical tool, constitutes a scientific opportunity to guide decision-makers towards favorable sectors for hydrogeological investigations and more transmissive areas in terms of water resources.   Informed Consent Statement: Not applicable.

Data Availability Statement:
The data used to support the conclusions of this study are available on request and by agreement from the corresponding author.

Acknowledgments:
The identification of transmissivity values will help us to know the hydrodynamic function of the chalk aquifer in the Oise region, especially in the Hydrogeological Experimental Site of Beauvais (HESB, northern part of the Paris basin). The construction of this site has benefited from the support of the European Regional Development Fund (FEDER), the Picardy region (Haut-de-France), and the Ministry of Higher Education and Research. This site and its hydrogeological platform are attached to the AGYLE Research team "Agroecology, Hydrogeochemistry, Environments, and Resources" of the Institut Polytechnique UniLasalle.

Conflicts of Interest:
The authors declare no conflict of interest. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data used to support the conclusions of this study are available on request and by agreement from the corresponding author.