Frederick Delay
email: fdelay@unistra.fr

Hamid Badri

Philippe Ackerer

Heterogeneous hydraulic conductivity and porosity fields reconstruction through steady-state flow and transient solute transport data using the continuous adjoint state

Keywords: Parameter estimation, adjoint state, adaptive parametrisation, groundwater, hydraulic conductivity, porosity

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.

Introduction

Parameter estimation through inverse methods for flow and solute transport simulations is still challenging despite the variety of concepts and methods developed since the pioneer works in inverse methods of [START_REF] Vemuri | Sensitivity analysis method of system identification and its potential in hydrologic research[END_REF], [START_REF] Emsellem | An automatic solution of the inverse problem[END_REF], [START_REF] Yeh | Quasi-linearization and the identification of aquifer parameters[END_REF], among others. Recent reviews of conceptual frameworks and methods can be found in Hendricks [START_REF] Hendricks Franssen | A comparison of seven methods for the inverse modelling of groundwater flow. Application to the characterization of well catchments[END_REF], [START_REF] Zhou | Inverse methods in hydrogeology: Evolution and recent trends[END_REF], and [START_REF] Linde | Geological realism in hydrogeological and geophysical inverse modeling: A review[END_REF]. Almost all the developed methods target the estimation of hydraulic conductivities (or transmissivities) constrained (conditioned) by hydraulic heads and, through the last 40 years, by additional data such as solute concentrations, geological information or geophysical properties.

In this work, we focus, as a first goal, on estimating spatially-heterogeneous hydraulic conductivity and porosity fields conditioned by both hydraulic heads and solute concentrations.

Many authors have already addressed this challenge, and Table 1 presents an overview of various approaches dedicated to this topic. The table is not an exhaustive presentation of the existing literature but a selection of works (17 contributions ranked as they appear in time in the literature) that we consider representative of methods and trends.

Solving the inverse problem for flow and transport in ground water systems is tightly associated with methodological choices. Among the possible choices, we exemplify, via continuous adjoint state calculations, how gradient-based inversion techniques can cope with highly parameterized problems. Mixing data of various types and weighting them in the objective function of an inverse problem is also a key feature. This is why this study emphasizes methodology as a second goal.

Notably, when modeling natural systems (or synthetic test cases supposedly close to natural systems), it is always complicated to discuss on identifiability, uniqueness of inverse solutions, quality of conditioning, etc., simply because highly parameterized inverse problems encompassing multiple elementary processes result in fuzzy theoretical frameworks. That being said, it cannot be overlooked that complex systems might result in multiple inverse solutions.

Repeatability, in the sense of estimating how similar various solutions are when initiated at different locations in the parameter space, appears as a simple way to determine success or failure of inversions. Incidentally, when dealing with synthetic test cases employing fullyknown reference problems, the comparison between inverse solutions and references is also informative.

In this study, observations associated with flow and transport parameter estimations are scarce piezometric head data from steady-state flow conditions and solute concentrations under transient transport conditions. These data are used to build an objective function assembling two terms: the sum of the weighted squared differences between measured and computed piezometric heads h F , and the sum of the weighted squared differences between measured and computed solute concentrations c F . Notably, the weights may sometimes serve to narrow down the range spanned by the squared differences when data of different units are used. This type of objective function is usually employed within: (i) the non-linear Least Square framework, which also often adds a regularization term on model parameters, or (ii) the Maximum Likelihood framework earlier proposed by [START_REF] Carrera | Estimation of Aquifer Parameters Under Transient and Steady State Conditions: 1. Maximum Likelihood Method Incorporating Prior Information[END_REF], which allows for handling prior information on model parameters. Summing two quantities that have different units in an objective function requires an adapted strategy and/or the use of a weighting coefficient applied to one of these terms (usually c F ). In fact, the weighting should be directly given by the probability density functions of errors between model outputs and observations. The point is that these densities are generally conjectured, which results in weights given to h F and c F becoming parameters of the inversion procedure. Pioneer works from [START_REF] Umari | Identification of aquifer dispersivities in twodimensional transient groundwater contaminant transport: An optimization approach[END_REF], [START_REF] Strecker | Parameter identification of a ground water contaminant transport model[END_REF], and [START_REF] Keidser | A comparison of four inverse approaches to groundwater flow and transport parameter identification[END_REF] proposed a two-stage procedure to minimize successively h F and c F . Simultaneous inversion of both heads and solute concentrations (and simultaneously decreasing both h F and c F ) is nowadays the common way to estimate the parameters. Nevertheless, the optimal weighting of the two terms h F and c F has been under debate [START_REF] Carrera | Estimation of Aquifer Parameters Under Transient and Steady State Conditions: 1. Maximum Likelihood Method Incorporating Prior Information[END_REF][START_REF] Doherty | Ground Water model calibration using pilot points and regularization[END_REF][START_REF] Medina | Geostatistical inversion of coupled problems: dealing with computational burden and different types of data[END_REF], and is still not fully resolved.

Even though abundant literature and many studies are available for other inversion techniques, gradient-based methods are the most common algorithms applied to minimize the objective function. The minimization employs either sensitivity coefficients with a Gauss-Newton method [START_REF] Levenberg | A Method for the solution of certain non-linear problems in least Squares[END_REF][START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF] or adjoint state variables coupled with a Quasi-Newton method [START_REF] Byrd | A limited memory algorithm for bound constraint optimization[END_REF]. Over the last decade, Kalman filtering [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] and its extensions such as Ensemble Kalman filtering (e.g., [START_REF] Iglesias | Ensemble Kalman methods for inverse problems[END_REF][START_REF] Deng | Identification of hydrological model parameter variations using ensemble Kalman filter[END_REF] have been applied and improved to efficiently solve the inverse problem. The method has become a valuable alternative approach to gradient-based methods for estimating flow and transport parameters.

Parameterization techniques have also shown an interesting evolution over the last 50 years. Zonation was the most popular method during the seventies period, probably due to lack of powerful computers that hampered the deployment of stochastic (Monte Carlo) simulations.

Stochastic parameterization now seems to be the standard approach for seeking hydraulic conductivity, but recently, [START_REF] Pool | A comparison of deterministic and stochastic approaches for regional scale inverse modeling on the Mar del Plata aquifer[END_REF] recalled that zonation could be a relevant option to represent geological patterns for large-scale aquifers (i.e., systems of several tens of km extension with flow patterns modeled at the 100-1000 m scale). In most applications, transport parameters (effective porosity and dispersivities) are usually considered as uniform over the modeled domain, and transport heterogeneity is mainly associated with heterogeneous hydraulic conductivity fields and subsequent variations of the mean fluid velocity in the system.

In this work, we explore the capacity of the adaptive parametrization technique detailed in [START_REF] Ackerer | Inversion of a set of well-test interferences in a fractured limestone aquifer by using an automatic downscaling parameterization technique[END_REF] and Hassane and Ackerer (2017) coupled with an independently computed adjoint equation [START_REF] Delay | A comparison of discrete versus continuous adjoint states to invert groundwater flow in heterogeneous dual porosity systems[END_REF] to estimate hydraulic conductivity and porosity spatial distributions conditioned by head and concentration measurements.

In the following section, the mathematical models describing the flow and transport processes are presented. Section 3 provides a synopsis of the continuous adjoint state method used to minimize the objective function under the constraints resulting from solving the flow and transport equations.

The forward problem

Steady-state flow and transient transport conditions in a 2D confined aquifer are assumed. These conditions are quite common for practical cases where steady-state flow may represent average flow conditions. These conditions also work for shallow unconfined aquifers that show small variations of water levels (which keep the flow equations linear regarding hydraulic heads). We also assume constant aquifer thickness for simplicity, which allows for inverting hydraulic conductivities instead of transmissivities, but without loss of generality.

Under these conditions, groundwater flow and solute transport are modeled over the domain Ω and, for transport, over the time interval ] ] 0,T by (e.g., [START_REF] Bear | Dynamics of Fluids in Porous Media[END_REF]:

. over over . over

w D D H w N N H q K h h h K h q Γ ∇ = Ω   = -∇   = ∂ Ω   -∇ = ∂ Ω  q q n (1) ( ) ] ] ] ] ] ] 0 .
; over and 0, over at 0 over and 0, . over and 0,

c D D C c N N C c ω c ω c q T t c c t c c T c q T ω ω Γ ∂  + ∇ - ∇ = = Ω  ∂   = Ω =   = ∂ Ω  -∇ = ∂ Ω   u D u q D n (2)
where h [L] is the hydraulic head; K [LT -1 ] is the hydraulic conductivity (considered here as a scalar); w q [T -1 ] is the flow source/sink term; c [ML -3 ] is the concentration; ω [-] is the effective porosity; q [LT -1 ] is the Darcy flux; u [LT -1 ] is the average water velocity; c q [ML - ( )

2 2 y x xx L T 2 2 y x yy T L x y yx xy L T u u D = + u u D = + u u D = D = - α α α α α α u u u u u (3)
where ( ) Boundary conditions and source/sink terms are assumed to be known and prescribed at constant values over time. It is noticeable that for the above equations, the porosity should not vary over time, otherwise flow would no longer be steady-state. Boundary conditions for flow and transport have also been written as independent from one another when it might not be the case.

For example, Neumann boundary conditions for flow with fluxes exiting the domain are often associated with Neumann boundary conditions for concentrations, but for entering water fluxes, concentrations are often prescribed via Dirichlet-type boundary conditions.

Parameter identification methodology

The objective function

The parameter identification procedure is based on the minimization of the quadratic differences between measured and computed variables. These differences are gathered in a socalled objective function defined as:

( ) ( ) ( ) , ( ), ( ) , ( ) , ( )
h c F F F κ = + p h p c p p h p p c p (4) with ( ) ( ) ( ) ( )
( ) ( ) 1 ˆ, ( ) ˆ, ( ) t T h h N T j j j j j c c j F F =  = - -     = - -   ∑ p h p h h W h h p c p c c W c c (5)
where () T denotes the transpose operator; h and j c are the vectors of computed heads and concentrations with a size corresponding to the number of observations; j is the index of observation times for transport; t N is the number of times at which concentrations have been measured; p is the vector of parameters (here, hydraulic conductivity and porosity) of size p N and ĥ and ˆj c are the vectors of observed values. ˆ- h h and ˆj j c c are usually named as "measurement" errors in the jargon of inverse problems, as they simply are the difference between model outputs and observations in a framework considering that the model is "exact" (for a given set of parameters) and data are "flawed". In essence, it must be understood that measurement errors encapsulate actual errors on measures, differences between model outputs and data, and the consequences of conceptual errors associated with a model or its settings that are always approximations of reality. Assuming that the measurement errors are not correlated in space and invariant over time, the W matrices are diagonal, and [START_REF] Medina | Geostatistical inversion of coupled problems: dealing with computational burden and different types of data[END_REF], and should therefore be modified during the minimization procedure.

Different minimization techniques of the objective function ( )

, ( ), ( ) F p h p c p exist.

Because the number of estimated parameters might be large, we rely upon the adjoint state method rendering an estimate of the gradient components of the objective function at a cost independent of the type and the number of parameters (see Appendix A). Even though the convergence toward a minimum of the Quasi-Newton methods is usually slower compared with the Gauss-Newton method [START_REF] Cooley | A comparison of several methods of solving nonlinear-regression groundwater-flow problems[END_REF], for highly parameterized problems, it is wise to rely upon the adjoint method seeking gradient components which are then introduced into a Quasi-Newton algorithm [START_REF] Chavent | Identification of distributed parameter system: about the output least square method, its implementation and identification[END_REF][START_REF] Neuman | Maximum-likelihood adjoint-state finite-element estimation of groundwater parameters under steady-and nonsteady-state conditions[END_REF][START_REF] Sun | Inverse Problems in Groundwater Modeling[END_REF][START_REF] Ackerer | Flow in double-porosity aquifers: Parameter estimation using an adaptive multiscale method[END_REF][START_REF] Sun | Model Calibration and Parameter Estimation[END_REF][START_REF] Delay | A comparison of discrete versus continuous adjoint states to invert groundwater flow in heterogeneous dual porosity systems[END_REF]. Variants of the adjoint state exist (e.g., [START_REF] Medina | Geostatistical inversion of coupled problems: dealing with computational burden and different types of data[END_REF] to calculate model sensitivity to parameters which then feed Gauss-Newton algorithms. The latter option is not employed here as the calculation costs strongly depend on the number of measurement values.

The continuous adjoint state

The adjoint state technique is associated with the class of problems ensuring the optimization (minimization) of an objective function under constraints (e.g., [START_REF] Bertsekas | Constrained optimization and Lagrange multiplier methods[END_REF]. These constraints can be, for example, the equations that are solved by a model for which we seek the optimal parameters rendering model outputs close to observations. Recently, [START_REF] Delay | A comparison of discrete versus continuous adjoint states to invert groundwater flow in heterogeneous dual porosity systems[END_REF] presented a thorough discussion on the adjoint state to invert the spatially distributed problem of flow in dual-porosity systems. Their discussion dealt with both the discrete and continuous forms of the adjoint state and with a comparison of their ability to assist inversions in providing rapid calculations of the gradients of the objective function. The present work inherits the method proposed by [START_REF] Delay | A comparison of discrete versus continuous adjoint states to invert groundwater flow in heterogeneous dual porosity systems[END_REF] and extends it to the inversion of solute transport treated via an advection dispersion equation. The adjoint state is employed in its continuous form, which can be differentiated, and then implemented, without prior knowledge about the structure of the discrete equations in the forward model.

Let us take a forward problem as a continuous operator ( )

, 0 c p = V handling a single state variable ( ) , c t
x and a single type of parameter ( )

, p t
x over the domain Ω and the period ] ] 0,T . References to space x and time t coordinates are dropped herein when feasible for better readability. The objective function of the inverse problem is usually discrete because it compares a finite set of observations of c with equivalent simulated values (see e.g., Eqs. ( 4) and ( 5)). Nevertheless, the objective function is rewritten in a continuous form

( ) ( ) 0 , , T F c p f c p d dt Ω = Ω ∫ ∫
, with f an elementary function, usually piecewise (or pointwise)

non-null at locations in time and space where observations are available. It must be noted that the simple integrals over time and space defining F are here employed for ease of notation, but without loss of generality. When dealing with discrete objective functions as those in Eq.

(5), the formulation of F becomes complicated with four integrals concealed in the continuous form (see, [START_REF] Delay | A comparison of discrete versus continuous adjoint states to invert groundwater flow in heterogeneous dual porosity systems[END_REF]. For example, with point observations of concentrations ˆb a c

denoting a value at location a x and time b t , a generalized least-square objective function on concentrations in a continuous form could be written as:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , , 0 0 ˆ, , ' ' ' T T i i j j m m n n m n i j m n i j c t c W c t c t t t t d d dtdt δ δ δ δ Ω Ω - - - - - - ∑∑∑∑ ∫ ∫ ∫ ∫
) ( ) ( ) 0 0 , , , , T T c p f c p d dt c p d dt µ µ Ω Ω = Ω + Ω ∫ ∫ ∫ ∫ L V (7)
As the second term in Eq. ( 7) is null, the Lagrangian L is equivalent to the objective function F , even though it encloses the addition of the constraintV multiplied by the variable ( )

V

, a variation of the Lagrangian L can be written as:

( ) ( ) ( ) 0 0 0 , , T T T f f c d dt p d dt c c p p c p d dt c p δ δ δ δ δ µ Ω Ω Ω ∂ ∂ = Ω + Ω + + + - Ω ∂ ∂ ∫ ∫ ∫ ∫ ∫ ∫ L V V (8)
where a first-order Taylor series expansion has been used to rearrange the terms in f .

Identifying the equations ruling the adjoint state µ relies upon the development of the variations in the constraint ( ) ( ) δ , which in turn modifies the variation of the Lagrangian that becomes:

( ) 0 0 , T T f p d dt c p d dt p δ δ µ δ Ω Ω ∂ = Ω + Ω ∂ ∫ ∫ ∫ ∫ L G (9)
with G a continuous operator combining the state variable c and the adjoint state µ , both previously calculated by the forward problem and the adjoint state equations, respectively. This operator G depends on the type of parameter (the perturbation p δ ) to which it is associated in the scalar product in Eq. ( 8). Usually, parameters in spatially distributed models are defined as uniform values over subdomains i Ω of Ω (e.g., zones, cells, etc.) and eventually over periods ,p t x to the subdomain i Ω and the period n t ∆ , the variation of the Lagrangian δL in Eq. ( 8) renders an approximation of the gradient of the objective function with respect to n i p ( )

, n i n i n n n i i i t t dF f d dt c d dt p dp p δ µ δ ∆ Ω ∆ Ω ∂ ≈ = Ω + Ω ∂ ∫ ∫ ∫ ∫ L G (10)
where integrals in Eq. ( 10) are restricted to the subdomain i Ω and the period n t ∆ . It is worth noting that the adjoint state µ is independent of the type of gradient component

Adjoint state equations and gradients of the objective function

A crucial point of the continuous adjoint state application is in developing the equations ruling the adjoint state and setting the operators G mentioned above for calculating the gradient components. Several attempts have appeared in the literature even though the technique has not been widely employed in Hydrology. The most important contribution is probably that proposed by Sun and co-workers [START_REF] Sun | Inverse Problems in Groundwater Modeling[END_REF] Regarding the inversion of the coupled problem of steady-state flow and transient transport, a variation of the continuous Lagrangian operator can be written as:

( ) ( ) ( ) ( ) ( ) ( ) ( ) 0 0 0 , , , , , , , , h h T
T T c c f f h c h d d h h h d h f f c d dt d dt c c c d dt c δ λ µ δ δ δ δ λ δ δ δ δ µ Ω Ω Ω Ω Ω Ω ∂ ∂ = Ω + Ω + + + - Ω ∂ ∂ ∂ ∂ + Ω + Ω + + + - Ω ∂ ∂ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ p p p p p p p p p p p L U U V V (11)
The elementary continuous objective function f has been dissociated in the two components h f and c f enclosing the state variable h (head) for flow and the state variable c (concentration) for transport, respectively;

( )

, 0 h p = U is the continuous steady-state flow equation; ( ) , 0 c p = V
is the continuous transient-transport equation; λ and µ are the continuous adjoint states associated with flow and transport, respectively; and p is the undifferentiated vector of factors ruling the flow and transport equations. It could include initial and boundary conditions, source/sink terms, and flow and transport parameters. Here, the aim is to retrieve hydraulic conductivity (and porosity) fields on the basis of head and concentration measurements for pre-identified dispersivity parameters, known source/sink terms, and known initial and boundary conditions. That being said, the hydraulic conductivity influences transport parameters, such as the mean water velocity u and the dispersion tensor D (see Eq. ( 2)). It makes sense to calculate via the adjoint state the gradient components of the objective function with respect to u and D and link them with the total variation of the objective function with respect to K.

The continuous adjoint state equations for steady-state flow ruled by Eq. ( 1) are derived as (see [START_REF] Delay | A comparison of discrete versus continuous adjoint states to invert groundwater flow in heterogeneous dual porosity systems[END_REF]:

( ) . 0 over 0 over . 0 over h D H N H f K h K λ λ λ Γ ∂  ∇ -∇ + = Ω  ∂  = ∂ Ω   ∇ = ∂ Ω   n (12)
Transport equations in Eq. ( 2) are flanked with continuous adjoint state equations in the form (see Appendix A):

( ) ] ] ] ] ] ]
c D C N C f ω T t c t T T T µ ω ω µ µ µ µ µ Γ ∂ ∂  - - ∇ -∇ ∇ + =  ∂ ∂   = Ω =   = ∂ Ω  ∇ = ∂ Ω   u D D n (13)
The minus signs appearing in Eq. ( 13) for terms in t ∂ ∂ and in .∇ u assume that Eq. ( 13) is solved backward over time from T to 0 and with a reversed water velocity field -u . Otherwise, Eq. ( 13) lacks physical meaning by letting a negative diffusion tensor appear.

As told earlier, continuous adjoint state equations look very similar to their equivalent forward problem in Eqs. ( 1) and (2). Eqs ( 12) and ( 13) can be solved numerically by any means (any discretization and numerical scheme) independent of the way the forward problem is solved. The forward model and the adjoint state calculations can be decoupled, which is not the case of the discrete adjoint state technique which is differentiated (from) and calculated with the same discrete equations as those of the forward model (e.g., [START_REF] Ackerer | Flow in double-porosity aquifers: Parameter estimation using an adaptive multiscale method[END_REF] .

h dF f d h d dK K λ Ω Ω ∂ = Ω+ ∇ ∇ Ω ∂ ∫ ∫ (14)
If there is no information (or no prior guess) on hydraulic conductivity values and their comparison with sought values, the term f K ∂ ∂ cancels out in Eq. ( 14) and the gradient component of the objective function is a simple scalar product between simulated heads and the associated adjoint state. For the part of the objective function handling concentrations, the gradient component with respect to K can be obtained as:

c c c dF F F dK K K ∂ ∂ = + ∂ ∂ u D (15) with ( ) 0 0
. .

T T T T c c c u F F F c d dt c h d dt K K ω µ µ Ω Ω ∂ ∂ ∂ ∂   = - ∇ Ω ⇒ = = ∇ ∇ Ω   ∂ ∂ ∂ ∂   ∫ ∫ ∫ ∫ u u u (16) ( ) ( ) 0 0 1 1 . . . T T T T c c y x x y x y F F c d dt c d dt K K u u h h K u K u K u x u x ω µ ω µ ω ω Ω Ω - - ∂ ∂ ∂ = ∇ ∇ Ω ⇒ = ∇ ∇ Ω ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = + = - - ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∫ ∫ ∫ ∫ D D D D D D D D (17) and ( ) ( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) 2 2 2 2 2 2 3 3 3 3 2 2 2 2 2 2 3 3 3 3 2 2 2 2 yy xx x x L x y T y T x y L y x x xy yx y L T x x y yy y xx T y x L x L y x T x y y xy yx x L T y y D D u u u u u u u u u u D D u u u u D u D u u u u u u u u D D u u u α α α α α α α α α α α α ∂ ∂ = + - = + - ∂ ∂ ∂ ∂ = = - ∂ ∂ ∂ ∂ = + - = + - ∂ ∂ ∂ ∂ = = - ∂ ∂ u u u u u u (18)
The gradient components of the objective function with respect to transport parameters, source/sink terms, and initial and boundary conditions are derived in Appendix A. As mentioned earlier, the non-intrusive character of the continuous adjoint state implies exchanges of information between the forward model and the adjoint state. This is also exemplified by the relations of Eq. ( 14) to Eq. ( 17), showing that gradient components of the objective function combine values of the state variables calculated by the forward model and values of the adjoint states.

Parameter estimation strategy

The objective function defined by Eq. ( 4) is a sum of two terms, depending on two very different quantities: piezometric heads expressed as a length and concentrations expressed as mass or moles per volume. In addition to differences in measurement units, the numerical values can also be significantly different. As an example, a difference of 1.0 m between computed and measured heads for heads varying between 100 m and 110 m is less significant than a difference of 1.0 mg.l -1 between computed and measured concentrations for values varying between 0.0 mg.l -1 and 2.0 mg.l -1 .

In the context of the Maximum Likelihood or Generalized Least Squares with regularization, the coefficient κ in Eq. ( 4) has a physical meaning if the priors on measurement errors are known. However, this information is not often available. [START_REF] Doherty | Ground Water model calibration using pilot points and regularization[END_REF] and [START_REF] Carrera | Estimation of Aquifer Parameters Under Transient and Steady State Conditions: 1. Maximum Likelihood Method Incorporating Prior Information[END_REF] suggested including this coefficient in the parameter estimation procedure. [START_REF] Medina | Geostatistical inversion of coupled problems: dealing with computational burden and different types of data[END_REF] analyzed the optimal weight for identifying hydraulic conductivities conditioned by piezometric heads and concentrations, and they tested different values. The expected value of the likelihood function provides robust weights to be assigned to concentration and hydraulic conductivity data, except for the steady-state flow case. We followed the suggestion of [START_REF] Medina | Geostatistical inversion of coupled problems: dealing with computational burden and different types of data[END_REF] and define the weighting coefficient as:

( ) ( ) 1 , ( ) , ( ) k k k h c F F θ κ + = + p h p p c p (19)
where θ is a user-defined variable, and k is the iteration index in the optimization procedure.

An alternative to evaluating hydraulic parameters using both types of data h and c is to consider ( ) -S1: hydraulic conductivities are estimated using head data only.

-S2: hydraulic conductivities are estimated using head and concentration data simultaneously.

-S3: the conditioning is alternated during the minimization in the sense that a first set of Mh iterations considers head data only, which is then followed by Mc

Numerical experiments

We rely upon synthetic test cases to evaluate the performance of the parameter estimation methodology. Measurement errors in data (as they could be defined here by adding random fluctuations to heads and concentrations extracted from the reference problems) are not considered in the inversion exercises as we mainly focus on comparisons between inversion strategies in their ability to identify reference parameter fields over diverse flow and transport conditions. In the context of inverse problems, the notion of measurement errors is often taken as the difference between outputs of a supposedly "exact" model and "flawed" data. Therefore, the absence of noisy measures does not go against the use of an objective function inheriting from the notion of measurement errors.

Settings of the test cases

The

2 log K σ =0.
12 . Hydraulic conductivity values vary between 3×10 -3 ms -1 and 10 -5 ms - 1 (Fig. 2). Usually, sequential Gaussian simulation techniques, when generating random fields over a rough mesh, do not fully match the prescribed covariance. In practice, the effective covariance (variogram) of the reference hydraulic conductivity field is slightly distorted and shows a correlation length of approximately 85 m, but the right prescribed variance. This feature is unimportant when the question is to retrieve a reference parameter field via inversion without any prior guess on its structure. In the case of inversion exercises seeking both hydraulic conductivity and porosity fields, the reference porosity field follows a normal distribution generated with a spherical covariance function that has a correlation length of 80 m and a variance 2 0.002

ω σ =
. Porosity values vary between 0.04 and 0.36 (Fig. 2), and, as for hydraulic conductivity, the effective covariance is slightly distorted with an effective correlation length of 100 m. Steady-state flow is considered. The test case C1 is designed without source/sink terms whereas five pumping wells and two injection wells (locations in Fig. 1) are added in the domain for the other test cases C2 and C3 (see Table 2). Pumping and injection rates are set to

=1

T α m, respectively. Notably, test case C4 re- handles the settings of test case C2, but adds an unknown heterogeneous porosity field to be retrieved in addition to the hydraulic conductivity field. For all test cases, the solute transport simulations are run under transient conditions with initial null concentrations within the system, and the simulations are performed over 120 days. The transport equations are solved using a combination of mixed and discontinuous finite elements methods, which limits the effects of numerical dispersion (Hassane et al., 2017).

Steady-state hydraulic heads and time-varying concentrations are observed at 20 wells (see Fig. 1) uniformly distributed in the domain. Concentrations are sampled every day. Overall, the observations consist of 20 head values and 20×120 = 2400 concentration values. It is worth noting that a large number of data points, as is the case here for concentration values, can mislead the inversion procedure when the corresponding computed variables are not sensitive to the parameter values, as, for example, close to the boundary conditions. The weighting parameter θ in the κ coefficient of the objective function (see Eq. ( 19)) is set to 10 -4 for each test case.

The initial parameter mesh of the adaptive parameterization (see, e.g., Ackerer and Delay, 2010) consists of 12 triangular elements and 11 nodes uniformly distributed over the domain. The parameter mesh is refined 3 times and the total number of estimated hydraulic conductivities, which may vary between runs according to the way the parameter mesh is refined, is between approximately 60 and 80. Each inversion is repeated 50 times, and the initial parameter grid is assigned with initial parameter values chosen randomly in a uniform distribution. The aim of this duplication is to evaluate the reproducibility of the results and compare them to the reference. The sought hydraulic conductivity values are bounded between -5

1.0×10 m.s -1 and -3

3.0×10 m.s -1 during inversion, and for test case C4, porosity is sought in the range [0.01, 0.40]. These upper and lower bounds of hydraulic conductivities and porosities are the only prior information on parameters brought to the inversion procedure.

Adjoint state variables associated with both flow and transport are computed with a different numerical code but with the same grid and the same numerical methods as those employed for the calculation of the state variables. Codes for the adjoint states were developed in the lab as tools independent of the forward codes with the aim of building an inversion toolbox working in parallel with any forward model, even on non-proprietary codes for which the calculation structures are unknown. Separate codes between the forward problem and the adjoint state calculations also allowed us to check how the exchanges of information between the grids of adjoint states and that of the forward flow and transport models could eventually hamper inversions (see e.g., [START_REF] Delay | A comparison of discrete versus continuous adjoint states to invert groundwater flow in heterogeneous dual porosity systems[END_REF]. This feature is not reported in this study as both the forward problems and adjoint state calculations share the same computation grids with exchanges of information free from any interpolation.

Inversion of the hydraulic conductivity with known uniform porosity

The settings of the test cases are motivated by three objectives that are: (i) to evaluate the different strategies for estimating the parameters (strategies S1 to S3), (ii) to assess the effects of flow conditions by comparing test cases C1 and C2, and (iii) to assess the effects of transport boundary conditions by comparing test cases C2 to C3.

Irrespective of the addressed test cases, the observed hydraulic heads are always matched within an error of 5-10 cm for maximal variations of heads over the flow fields of approximately 2 m. This very good matching (see also a comparison of heads between inverted and reference flow fields in Fig. 3) justifies that results on heads will not be addressed in the following discussion. Consequently, the comparisons are based on different criteria and only involve concentration and hydraulic conductivities. These criteria deal with the overall spatial distribution of concentrations at a given time, the analysis of breakthrough curves (BTC) at 3 locations (see Fig. 1), and the calibrated parameter fields compared with the reference. Fifty different solutions per test were analyzed. These solutions were obtained by starting the inversion procedure from different randomly picked locations in the parameter space. It is also worth noting that the multiscale parameterization randomizes the current inverse solution each time the parameter grid is refined. These solutions cannot be distinguished in view of the final values of the objective function, meaning that they are of similar quality in terms of quadratic differences between the reference and estimated heads and concentrations.

Regarding solute concentrations, the BTC stemming from 50 solutions were gathered into a confidence interval computed for each time step n and defined by was computed for each element E of the mesh with its associated standard deviation. An average error for each element E E ε was also used as an assessment criterion and defined by:

( ) ( ) log log E E E K K ε = - % (20)
where E K is the exact value (from the reference field) of the hydraulic conductivity at element E , and

E K %
is the estimated value. The variograms of ( ) log E K % were also compared with the reference as an overall indicator of quality of estimated parameter fields. In the same way as for breakthrough curves, the 50 variogram functions were gathered into a confidence interval for each lag.

The concentration distributions are provided in Fig. 3. The main features of the reference concentration distributions are reproduced by the inverse solutions for all test cases and all inversion strategies (regarding the objective function). However, the estimated concentrations show additional smearing compared with the reference even though the same mesh and same dispersivities were used for calculating reference and inverse solutions. As shown hereafter, the retrieved hydraulic conductivity fields are smoother than the reference.

This feature is consistent with the robustness of the flow equation with respect to hydraulic conductivity, in the sense that the head variable is not very sensitive to local contrasts in hydraulic conductivity values (e.g., [START_REF] Giudici | Spectral analysis of the balance equation of groundwater hydrology[END_REF].

The strategy overlooking concentration data in the objective function (S1) provides results as good as the other strategies accounting for concentrations in test cases C1 and C3; however, in test case C2, the effects of the injection wells are not completely reproduced. The quality of the match for test case C3 is still quite good, even though solute concentrations remain null for some observation wells that are not downstream of the local injection source. This means that including observations at locations that do not see any concentrations is helpful to retrieve the fluid velocity field and its subsequent relationship with hydraulic conductivities.

The different strategies can be distinguished when considering BTC (Figs. 4 to 6) as a quality criterion of the parameter estimation. For test case C1 (Fig. 4), the BTC are well reproduced when employing the two strategies S1 and S2. For strategy S3, which consists of incorporating concentration data in the objective function while calculating the first iterations of convergence by relying upon head data only, the method results in a few discrepancies regarding the arrival times of concentrations. This bias on arrival times also occurs with strategy S1 (concentration data are never used) and test case C2 (Fig. 5). These two attempts show that concentration data provide useful and valuable information on solute travel times, which in turn inform on hydraulic conductivity values.

Test case C3 is harder to decipher because solute concentrations remain equal to the initial conditions (null concentration values) at a significant number of monitored wells not located downstream of the local injection. Therefore, test case C3 is more sensitive to flow direction variations within the domain than test cases C1 and C2. A poor estimation of the flow direction in test case C3 leads to either a difference in mean travel time or a difference in the total solute "mass" passing through the wells (the mass being defined as the integral of the BTC over time). Discrepancies in mass can be observed at wells P7 and P9 (except for strategy S2).

Again, strategy S2, which includes weighted information on both head and concentration data in the objective function, is the most efficient strategy to retrieve the flow field. It is worth noting that for all test cases and all strategies, the early arrival times of concentrations are well estimated. The accuracy of the estimated parameters decreases with the travel distance (late arrival times) as shown by the increase in the range of the confidence interval from well P11 to well P9 for test cases C1 and C2.

The estimated hydraulic conductivities are analyzed by comparing the average value of ( ) log K , its related standard deviation, and the average error in Eq. ( 20) (Figs. 7 to 9). The three test cases provide very similar results, showing that:

-In general, the hydraulic conductivity field is fairly well reproduced irrespective of the inversion strategy employed.

-The largest differences between estimated and reference ( ) -The third strategy, S3, leads to the poorest parameter estimation compared with the other two.

log
-The ensemble of inverse solutions is wider for strategy S2 than for strategy S1

(see test cases C2 and C3 in Figs. 8 and9). These results may appear counterintuitive as it is usually expected that more information (of good quality) brought to the inverse problem should yield better results, or at least more constrained solutions (which is consistent with the robustness of these solutions). This notwithstanding, the additional information provided by concentrations does not have the same meaning as information on heads. In short, concentrations are sensitive to flow directions and associated travel times along these directions as heads are sensitive to conductivities in a neighborhood roughly centered on the head measurement location. Increasing the complexity of the inverse problem by multiplying the number of sensitive factors or phenomena is conducive to an increase of local minima. This feature spreads out the set of possible inverse solutions even though these solutions resulted in similar values of the objective function.

Finally, the results were analyzed via a global criterion, comparing the variogram of ( ) log K from the reference solution to the variograms computed from the estimated solutions (Fig. 10). All variograms from inverse solutions underestimate the spatial variability of hydraulic conductivities within lag distances less than one-half the size of the domain. This feature is consistent with the fact evoked above stating that hydraulic head fields are not very sensitive to local contrasts in hydraulic conductivity values. In addition, concentrations at a given time and location depend on the average travel time between the source and the observation well. Therefore, they are rather less sensitive to local hydraulic conductivity variations than to the mean of conductivities experienced along the flow path between injection and observation. This could explain the inference of smaller variances for variograms from inverted conductivity fields. In the same vein, hydraulic heads, as the variables of a diffusion process (Darcian flow), are known to be rather more sensitive to mean values of conductivities over large patches than to contrasts between local values.

Identification of hydraulic conductivity and effective porosity distributions

The additional test case C4 re-handles the settings of test case C2 but considers the inversion of both the hydraulic conductivity field and a heterogeneous effective porosity field.

A single strategy S2 was employed to estimate the parameters, with the minimization in a single phase of an objective function mixing both hydraulic head and concentration data (Eq. ( 4)).

This choice of a single strategy is motivated by the inversions of hydraulic conductivity fields that showed strategy S2 to be the most efficient (even though the three strategies employed render almost similar results). It was thought that strategy S2 would also be the best suited to the joint inversion of hydraulic conductivity and porosity. As the mean porosity of the heterogeneous field is approximately 0.12-0.15 in the western half of the domain compared with the uniform value of 0.05 in test case C2, the simulation time for transport is increased up to 280 days (120 days for test case C2) to let the solute plume widely invade the system.

The reliability of the simulated concentrations (Figs. 11 and12) and of the estimated parameters (Figs. 13 and14) has been significantly diminished compared with the test cases solely inverting hydraulic conductivities. This can tentatively be explained by several features, among which the most likely is that the main driver to solute concentration propagation is the mean water velocity, which depends on the ratio of hydraulic conductivity to effective porosity.

It is always challenging to identify both terms of a ratio when the latter is the sensitive parameter, even though both terms of the ratio in the present study are independent quantities with different correlation lengths, etc. It is worth noting that actual data may show relationships between porosity and conductivity values, following for example the well-known Archie's law in its various forms. However, these relationships mainly apply to the fine scale of core samples, but rarely apply to the elementary mesh of a regional aquifer model. In the latter case, small variations in porosity values do not generate enough variability in hydraulic conductivities that might vary over orders of magnitude. That being said, it seems that choosing independent porosity and conductivities values does not help/hamper inversion more than choosing loosely correlated values. Another feature hampering inversion is the increased number of sought parameters, which is approximately 70 hydraulic conductivity values in test case C2 inversions and becomes approximately 210 (in our tests with two independent parameter grids, one for conductivity and one for porosity) in test case C4 inversions.

The difficulties experienced by the inverse procedure in test case C4 are also witnessed by the important variability in the estimated concentrations resulting from flow and transport through solution parameter fields (Fig. 12). The spreading of the 50 different possible solutions, and the mean local errors between solutions and the reference are also enlarged (compare, for ( ) log K plots, Fig. 13 and test case C2-S2 in Fig. 8). The variograms of the solution parameter fields (Fig. 14) render the same appraisal by underestimating the variability at small scales (small lag-distances). They also highlight the difficulty in estimating with confidence the spatial distribution of porosity.

Conclusion

We have provided a rigorous and detailed presentation of the partial differential equations ruling the continuous adjoint state associated with solute transport. The adjoint state obeys an advection-dispersion equation that must be solved backward in time and with a reversed flow field compared to the forward problem of transport. It also obeys specific boundary conditions. It can be identified uniquely in a single calculation very similar to that of the forward problem; then, it can serve as calculations for all gradients of the objective function with respect to parameters, initial and boundary conditions, and source/sink terms. We have also extended the adaptive parameterization technique to solute transport and to coupled inversion of flow and transport data. The methodology has been applied to four synthetic test cases targeting the estimation of hydraulic conductivity and effective porosity fields, using both piezometric heads and solute concentration data. These numerical experiments showed:

1. Employing the continuous adjoint state solved via an independent numerical code, that is, in a non-intrusive way with regard to existing (proprietary) transport models, is a viable alternative to more common methods (Gauss-Newton and Quasi-Newton methods with embedded computations of sensitivities or adjoint variables) for joint parameter estimation.

2. Concentration data has to be taken into account in the parameter estimation and included in the objective function in a balanced way. If the objective function is not weighted, usually hydraulic head observations are rapidly fitted by the inverse procedure, and the conditioning by concentrations will not bring any additional information. In the weighting procedure, the weighting coefficient can be adapted during the minimization and is not necessarily a parameter to include in the minimization.

3. Concentration data render information on travel times and flow directions as a valuable means to improve the inference of unknown hydraulic conductivity fields. Null concentrations at some locations are also helpful. In the process of incorporating transport data to retrieve conductivity fields, it could be conjectured that a poor prior knowledge of the fluid velocity field could mislead the inversion relying upon transport data. The consequence would be to partly separate the inversion procedure in a first step, predefining the flow field (conductivities) on the basis of hydraulic head measurements; then, in a second step, launching calculations of transport. Even though successes of inversions are problem dependent, we have shown that the above conjecture is wrong. It is wise to invert both flow and transport scenarios in the same procedure, which means that the objective function of the inverse problem should mix information on heads and concentrations. The objective function would then be minimized in a single phase, handling both terms in heads and concentrations.

4. In the test cases reported by this study, it is worth noting that hydraulic head observations were always rapidly fitted, and when used alone, these data rendered valuable hydraulic conductivity fields. The addition of concentration data in the inversion procedure brought few cosmetic features to the conductivity fields, except for better assessment of early and mean travel times in the system, which applies to both sweeping-uniform or tortuous flow fields. A key question raised by this study can be formulated as: Is it worth a try to include concentration data for identifying hydraulic conductivities when it is known how cumbersome and costly concentration measurements can be?

( ) ( )( ) ( )( )( ) ( ) ( )( ) ( ) ( ) ] ]
D D C C N N C c c c t c c c T c c q T δ δ δ δ Γ  + = Ω =  + = ∂ Ω   -+ ∇ + = ∂ Ω  D D n (A5)
Developing to the first order (i.e., by neglecting terms in 2 δ ), the difference δ -V V renders:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ] ] , , 0 . . .
+ + - = ≡ ∂ ∂ + + ∇ + ∇ + ∇ ∂ ∂ -∇ ∇ -∇ ∇ -∇ ∇ = Ω p p p u u u D D D V V (A6)
with initial and boundary conditions:

] ] ] ] 0 over at 0 0 over and 0, . . . . 0 over and 0,

D C N C c t c T c c T δ δ δ δ Γ Γ  = Ω =  = ∂ Ω   -∇ - ∇ = ∂ Ω  D n D n (A7)
The basic idea prevailing for the identification of the continuous adjoint state equations is to isolate the terms in δ in the integrals

( ) 0 T d dt δ µ Ω - Ω ∫ ∫ V V
of the operator δ L in Eq.

(A3). This occurs via integration by parts for all terms in Eq. (A6), which are handled below as they appear ranked in Eq. (A6).

( )

0 0 0 0 T T T T c d dt c d dt c d t t c d dt t ωδ µ µ ωδ ωδ µ µ ωδ Ω Ω Ω Ω   ∂ ∂ Ω = - Ω + Ω   ∂ ∂   ∂ = - Ω ∂ ∫ ∫ ∫ ∫ ∫ ∫ ∫ (A8)
In Eq. (A8), simplifications come from the following properties: at 0 t = , 0 c δ = is due to boundary conditions in Eq. (A7); at t T = , µ is set to zero as it is the initial condition for calculating µ (see below), and µ can be defined up to the addition of a constant.

Ω Ω Ω Ω Ω   ∂ ∂ Ω = - Ω + Ω   ∂ ∂   ∂ = - Ω - Ω ∂ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ (A9)
where it has been accounted for the fact that 0 µ = at t T = . We note that:

[ ]

0 0 0 0 0 0 0 T T T T T a b b bdt adt a b adt a b if b t t t ∂ ∂ ∂ = - + = - - = ∂ ∂ ∂ ∫ ∫ ∫ (A10)
Using the property 0 µ = at t T = and applying Eq. (A10) in Eq. (A9) results in:

( )

0 0 T T c c d dt d dt t t δω µ δωµ Ω Ω ∂ ∂ Ω = Ω ∂ ∂ ∫ ∫ ∫ ∫ (
ω δ µ δ ω µ δ ω µ ωµ δ δ ω µ δ ω µ Γ Ω Ω ∂Ω Γ ∂Ω Γ Ω ∂ Ω -∇ ∇ Ω = - ∇ ∇ Ω + ∇ Γ - ∇ Γ = - ∇ ∇ Ω + ∇ Γ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ D D D n D n D D n (
ωδ µ ω µ δ ωµδ ω µ δ Γ Ω Ω ∂Ω Ω -∇ ∇ Ω = ∇ ∇ Ω - ∇ Γ = ∇ ∇ Ω ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ D D D n D (A16)
In Eq (A16) the contour integral cancels out because

Γ Ω Ω ∂Ω Ω -∇ ∇ Ω = ∇ ∇ Ω - ∇ Γ = ∇ ∇ Ω ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ D D D n D (A17)
In Eq (A17), the contour integral cancels out because Solving the adjoint state equations in Eq. (A18) changes the expression of the variation of the Lagrangian operator in Eq. (A3), which simplifies into:

D C N C f ω ω t c t T ω µ ω µ µ µ µ µ ω µ Γ ∂ ∂  - - ∇ -∇ ∇ + = Ω  ∂ ∂   = Ω =   = ∂ Ω  + ∇ = ∂ Ω   u D u D n (
( ) 0 0 , T T f d dt c d dt δ δ δ µ Ω Ω ∂ = Ω + Ω ∂ ∫ ∫ ∫ ∫ p p p L G (A19)
With the three types of parameters ω , u, and D in p, the operator δ L can be separated into three components:

From Eqs. (A11), (A14), and (A17):

( )

0 0 . . . T T T f c d dt c c d dt t ω δ δω δω µ µ µ ω Ω Ω ∂ ∂   = Ω + - ∇ + ∇ ∇ Ω   ∂ ∂   ∫ ∫ ∫ ∫ u D L (A20)
From Eq. (A13):

0 0 . T T f d dt c d dt δ δ ω δ µ Ω Ω ∂ = Ω - ∇ Ω ∂ ∫ ∫ ∫ ∫ u u u u L (A21)
From Eq. (A16):

( )

0 0 . . T T T f d dt c d dt δ δ ω µ δ Ω Ω ∂ = Ω + ∇ ∇ Ω ∂ ∫ ∫ ∫ ∫ D D D D L (A22)
The expressions in Eqs. (A20), (A21), and (A22) give the form of the gradient components of the objective function with respect to parameters ω , u, and D . For example, let us take ,

n i n i n n n x i x i x i t t dF f d dt c d dt u u u x δ µ ω δ δ ∆ Ω ∆ Ω ∂ ∂ ≈ = Ω - Ω ∂ ∂ ∫ ∫ ∫ ∫ L (A23)
The strength of the continuous adjoint state also allows us to retrieve source/sink terms, initial conditions, and boundary conditions by providing the gradient components of the objective function with respect to these factors. Regarding source/sink terms, the operator δ -V V in Eq. ( A6) is modified to include a perturbation c q δ -. As this perturbation is independent of that on the parameters δω , δ u , and δ D , the integration ( ) unmodified. After the adjoint state µ has been calculated by solving Eq. (A18), the gradient components of the objective function with respect to source/sink terms c q are calculated by employing:

0 T d dt δ µ Ω - Ω ∫ ∫ V V
0 0 c T T c c q c f q d dt q d dt q δ δ δ µ Ω Ω ∂ = Ω - Ω ∂ ∫ ∫ ∫ ∫ L (A24)
Regarding initial and boundary conditions, the operator δ -V V in (A6) is modified on its initial and boundary conditions in Eq. (A7). The latter become:

] ] ] ] 0 over at 0 over and 0, . . . . over and 0,

D D C c N N C c c t c c T c c q T δ δ δ δ δ δ δ Γ Γ  = Ω =  = ∂ Ω   -∇ - ∇ = ∂ Ω  D n D n (
Ω Ω Ω Ω Ω   ∂ ∂ Ω = - Ω + Ω   ∂ ∂   ∂ = - Ω + Ω ∂ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ (A26)
where the properties

Γ Γ ∂Ω ∂Ω Γ ∂Ω Γ Γ ∂ Ω ∂ Ω Γ ∂ Ω ∇ Γ - ∇ Γ - ∇ Γ = ∇ Γ + ∇ Γ + Γ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫ D n D n D n D n D n n (A27)
The contour integrals of Eq. (A27) take into account the following properties: Table 1. 17 reference studies dealing with coupled flow-transport inversions. Main characteristics of these studies.

  the local components of the fluid velocity, u is the velocity the longitudinal and transverse dispersities [L], respectively. Dirichlet boundary conditions apply to contours D H ∂ Ω for flow and D C ∂ Ω for transport as Neumann boundary conditions apply to contours N H ∂ Ω for flow and N C ∂ Ω for transport. Γ n represents the outer vector normal to the contours where Neumann boundary conditions are applied.

  that both the objective function and the constraints are gathered in a Lagrangian operator in the form:

(

  restriction of the parameter ( )

  differently and knowing that µ has been calculated beforehand in a single step very similar to that of the forward model (see below), the access to gradient components is very rapid, irrespective of the type and the number of sought parameters. This renders the adjoint state technique suited to highly parameterized problems, when methods based on model sensitivities (the calculation of each sensitivity being very similar to the forward problem) are plagued by computation costs.

  , who developed the continuous adjoint states for single-phase flow, two-phase flow, transport under various conditions, and specific applications targeting the identification of a single type of parameter, a specific form of the objective function, or steady-state problems.[START_REF] Delay | A comparison of discrete versus continuous adjoint states to invert groundwater flow in heterogeneous dual porosity systems[END_REF] proposed a unified development (and presentation) of the continuous adjoint state for flow in dual porosity systems, which became available for the identification of all types of parameters, initial and boundary conditions, and source/sink terms. The continuous adjoint state equations for steady-state single-phase flow in Eq. (1) are an adaptation from[START_REF] Delay | A comparison of discrete versus continuous adjoint states to invert groundwater flow in heterogeneous dual porosity systems[END_REF]. The interested reader is referred to that work for further details. Regarding the transient solute transport in Eq. (2), the unified development of the adjoint state equations do not yet appear in the literature and is reported in Appendix A.

  iterations handling both head and concentration data; this swap is repeated until convergence. In this work, we set Mh=20 and Mc=5. Alternated conditioning suggests that a highly flawed flow field could result in simulated concentrations far from data and eventual problems of convergence. Thus, iterations handling head data only should first rough-out the flow field, then refined by using both heads and concentrations. The number of iterations for each type of conditioning is problem dependent. This number was chosen by trials counting solutions that converge for 10 runs started at different locations in the parameter space, and also counting the mean number of iterations required per solution to converge. Notably, we did not delineate any strategy seeking hydraulic conductivity on the basis of concentration measurements only. In general, when concentrations are measured in open wells over a groundwater system, measures of heads in those wells are also available. It would be counter-productive to discard head data. On the other hand, there exist many actual systems where concentrations are never monitored. Another reason is that seeking hydraulic conductivity fields only via concentration measurements is usually not a convergence problem.Concentration distributions over time and space depend on the velocity field, a quantity defined as the product of conductivity and head gradient. Lack of information on heads can result in many hydraulic conductivity fields rendering the same model outputs. In other words, concentrations mainly inform on transit times between two locations along a flow line. This time is insensitive to various distributions of the velocities along the line, provided their harmonic means are similar. Therefore, with smoothly varying head gradients, many distribution of hydraulic conductivities along a flow line would result in the same transit times.That being said, heads along a flow line are also mainly controlled by the harmonic mean of conductivities. The point is that the flow equation is diffusive. Widespread head data inform on mean conductivities along all the segments joining the pairs of observations with a consequence of a better deciphering of hydraulic conductivity distributions.

3 2×10 m 3

 33 synthetic examples investigate a rectangular flow domain of 250 m in length and 150 m in width. The aquifer thickness is constant at 20 m. The domain is discretized by a mesh with 1,566 triangular elements of 8 m average length. Four test cases are designed with increasing complexity in the flow field and the solute transport conditions (Table 2). The boundary conditions for flow are a prescribed uniform head of 100 m at the west boundary and a constant over time uniformly distributed outflow rate of -.s -1 over the east boundary. The north and south boundaries are considered as impervious (see Fig 1). The hydraulic conductivity field follows a log-normal statistical distribution generated with an exponential isotropic covariance function with an effective correlation length of 60 m and a variance of ( )

  m 3 .s -1 . The mathematical model for flow is solved using mixed finite elements[START_REF] Younes | Mixed finite elements for solving 2-D diffusion-type equations[END_REF] with a code developed in the lab and checked for being rigorously massconservative at the scale of each cell of the mesh.It is worth noting that all the settings employed to generate the various flow fields are assumed to render a well-posed problem for the inversion of hydraulic conductivities. The existence of source/sink terms or Neumann boundary conditions in the flow equations avoids the evaluation of conductivity up to the multiplication by a constant. Nevertheless, with only 20 head values available for the identification of more than 1500 parameters (in essence, one per cell of the discretized domain, noting those 1500 values are not independent), good inversion results are not guaranteed. This motivates to complement the head data set with concentrations from transport scenarios. Incidentally, rendering the flow field tortuous by adding injection/pumping wells in the system (tests cases C2 and C3) is supposed to complicate the inversion using head data only. In these test cases, concentrations following the main flow paths are expected to increase the degree of improvement of the inverse solutions conditioned on heads only.Regarding the simulation of solute transport, a null dispersive flux is set at all boundaries except at the line x= 0 m, of 1 kg.m -3 is prescribed for test cases C1 and C2, and at the upstream (western) flow boundary of test case C3 where the solute concentration is set to zero. For test case C3, a pulse of solute is injected over a rectangular domain of delineation kg.m -3 .s -1 . The transport parameters are homogeneous and assumed to be known for test cases C1, C2, and C3. The effective porosity is equal to 0.05 and the longitudinal and

  computations show that the confidence intervals did not change significantly by running more simulations. With these 50 realizations, an average value of the decimal logarithm of the hydraulic conductivity

  K are located close to the eastern boundary of the domain assigned with a prescribed flux for flow. At these locations, fluxes remain almost similar irrespective of the local conductivity values. The reliability of the estimated parameters, measured by the local standard deviation of ( ) log K , is poor (high standard deviation) in this eastern area compared with the rest of the domain. It is also worth noting that after 120 days of transport, solute concentrations did not reach the eastern boundary and thus poorly inform the inverse problem on local hydraulic conductivity experienced during transport. This could be an additional explanation for the poor evaluation of conductivities close to the eastern boundary of the domain.

µ

  A15)The expression in Eq. (A15) has been integrated by parts twice to extract the perturbation c δ from the operator .∇ ∇. The simplification of contour integrals results from 0 appear in (A16) results in a null term over N C ∂ Ω (see boundary conditions (A7)).

0µ=

  over the contour D C∂ Ω , and the perturbation δω is zero over the whole contour ∂Ω . Gathering all terms in c δ over the expressions in Eqs. (A8) to (A17), then reintroducing them into the variation of the Lagrangian operator δ L in Eq. (A3) and cancelling out the whole, poses the continuous equations of the adjoint state, which render:

.

  A18) It is worth noting that Eq. (A18) has to be solved backward in time and with a reversed velocity field -u ; otherwise, the problem lacks physical meaning with a negative dispersion tensor in the advection-dispersion equation. The Neumann boundary condition of the forward problem in Eq. (A2) has been transformed into a Robin boundary condition (with a reversed velocity field -u ) along the contour N C∂ Ω . Usually in transport problems, a so-In that case, the Robin boundary condition for the adjoint state simplifies into of a "free" advective flux exiting the system with no modification beyond the boundary. In that case, the Robin boundary condition of the adjoint state with a reversed flow field simply consists in letting a prescribed value 0 µ = enter into the domain.

  the x direction over a subdomain i Ω of Ω and over the time step n t ∆ within the period ] ] 0,T . The gradient component of the objective function with respect to ,

  latter does not alter the development between Eqs. (A8) and (A17). Consequently, the continuous adjoint state equation (A18) and the gradient components of the objective function (precisely, the components of the operator δ L ) with respect to parameters ω , u, and D (Eqs. A20, A21, and A22) are

  the two successive integrations by parts of Eq. (A15) and the integration by parts of (A16). More precisely, only the contour integrals are modified; their sum from expressions in (A15) and (A16) renders:

  Figure and table captions

Fig. 1 .

 1 Fig. 1. Flow domains of the three synthetic test cases C1, C2, and C3. Dirichlet boundary conditions for flow at western (red bolt line, h = 100 m) and eastern (blue bolt line, h = 90 m) boundaries. North and south boundaries are of no-flow type. Green lines in test cases C1 and C2 delimit the boundary associated with continuous uniform injection of solute. The green rectangle in test case C3 is the location of a short-step injection of solute. Observations of hydraulic heads and concentrations are located at the square dots. Blue and red dots mark the locations of extraction and injection wells, respectively. Wells P7, P9, and P11 are locations where solute concentration breakthrough curves are compared between reference and inverse solutions.

Fig. 2 .

 2 Fig. 2. Reference log hydraulic conductivity (logK) and porosity (ω) fields.

Fig. 3 .Fig. 4 .

 34 Fig. 3.Reference and examples of simulated head fields (constant over time) and solute plumes after 120 days for flow and transport test case C1, 80 days for test case C2, and 50 days for test case C3. S1, S2, and S3 denote the strategies employed for inversion, and h or c refer to either head or concentration fields. Black lines through some fields mark the main streamlines of the flow field. The reference flow fields for test cases C2 and C3 are similar. The color scale ranks values for both head and concentration casted between 0 and 1 for minimal and maximal reference values, respectively Fig. 4. Solute concentration breakthrough curves at the three locations P11, P7, and P9 (ranked with increasing distance from solute source) for test case C1 and three inversion strategies S1, S2, and S3. Black lines denote the references; blue areas delineate the envelope

Fig. 5 .

 5 Fig. 5. Solute concentration breakthrough curves at the three locations P11, P7, and P9 (ranked with increasing distance from solute source) for test case C2 and three inversion strategies S1, S2, and S3. Black lines denote the references; blue areas delineate the envelope

Fig. 6 .

 6 Fig. 6. Solute concentration breakthrough curves at the three locations P11, P7, and P9 (ranked with increasing distance from solute source) for test case C3 and three inversion strategies S1, S2, and S3. Black lines denote the references; blue areas delineate the envelope

Fig. 7 .

 7 Fig. 7. Average log hydraulic conductivity

Fig. 8 .

 8 Fig. 8. Average log hydraulic conductivity

Fig. 9 .

 9 Fig. 9. Average log hydraulic conductivity

Fig. 10 .

 10 Fig. 10. Variogram of the reference log hydraulic conductivity (black lines) compared with the variogram envelopes

Fig. 11 .

 11 Fig. 11. Reference (left) and simulated (right) solute plumes after 280 days for flow and transport test case C4.

Fig. 12 .

 12 Fig. 12. Solute concentration breakthrough curves at the three locations P11, P7, and P9 (ranked with increasing distance from solute source) for test case C4. Black lines denote the references; blue areas delineate the envelope

Fig. 13 .

 13 Fig. 13. Average parameters -log hydraulic conductivity ( ) log K and porosity ω -(left), standard deviation of parameters (middle), and mean error in parameters compared to reference (right) for 50 inverse solutions of test case C4.

Fig. 14 .

 14 Fig. 14. Variograms of the reference (black lines) log hydraulic conductivity (left) and porosity (right) compared with the variogram envelopes

  Fig. 1.

  The detailed mathematical development of the continuous adjoint state for transport is reported in Appendix A. It is presented in a new, convenient manner (as done for flow in[START_REF] Delay | A comparison of discrete versus continuous adjoint states to invert groundwater flow in heterogeneous dual porosity systems[END_REF] that allows for the derivation of the gradient of the objective

function with respect to all types of factors, including transport parameters, initial and boundary conditions, and source/sink terms of the transport problem. The objective function and the minimization strategies are presented in Section 4. Finally, the numerical experiments are presented and discussed in Section 5; emphasis is put on the comparison between inverse solutions obtained from calculations conditioned by head data only and by the joint use of head and concentration data.

  . For its part, the independent implementation of the continuous adjoint state is able to work with any forward model. Both simply exchange information, such as the source/sink terms of the adjoint . For the part of the objective function handling measured and simulated heads, the gradient component with respect to K is given by (see[START_REF] Delay | A comparison of discrete versus continuous adjoint states to invert groundwater flow in heterogeneous dual porosity systems[END_REF]

	states, h f h ∂ ∂ and c f c ∂ ∂

, those inheriting from calculations performed by the forward model. Notably, Eqs. (12) and (13) show that initial and boundary conditions for the adjoint states are always null. Stated differently, we can say that adjoint states let local source/sink terms diffuse or transport over domains that are partly disconnected from the domain of the forward problem. One might, for example, invert part of the forward problem on a restriction of its domain of definition by simply calculating the continuous adjoint states on this restriction with appropriate types of boundary conditions. With steady-state flow and transient transport aimed at retrieving hydraulic conductivity fields with known porosity and dispersivity transport parameters as well as known initial and boundary conditions, the gradient components of the objective function are of two types: h dF dK linking heads with hydraulic conductivity, and c dF dK linking concentrations with hydraulic conductivity via the mean water velocity u and the dispersion tensor D

Table 2 .

 2 Main characteristics of synthetic tests cases C1 to C4 inverting coupled flow and transport to retrieve hydraulic conductivities (and porosity, only test C4).

Acknowledgments

The authors are indebted to the Qatar National Research Foundation (QNRF) for partly funding this study under the grant NPRP9-030-1-008 "Dual porosity modeling of Qatar karst aquifers".

This study is also a contribution to the GEOTREF project funded by the French government in the framework of the "Investissements d'Avenir" program and tutored by ADEME.

Appendix A

The development of the continuous adjoint state associated with the advectiondispersion equation for solute transport in porous media is presented below. The continuous expressions of the gradient components of the objective function are also provided.

The continuous operator that solves solute transport and that represents the constraint ( )

in the Lagrangian operator of an optimization problem is written as:

At this stage, perturbations δ p resulting in perturbations c δ to keep ( )

δ δω δ δ = p u D , noting that u is a parameter for the transport equation but also an output of a flow model. A later discussion will address perturbations regarding source/sink terms, initial conditions, and boundary conditions. A variation of the Lagrangian operator is employed to derive the equations of the adjoint state, which takes the form:

with f the continuous objective function, µ the continuous adjoint state, and the constraint ( )

. The latter is developed as:

In Eq. (A12), the contour integral over the Dirichlet-type boundary D C

∂ Ω has been cancelled out, considering that 0 c δ = along this type of boundary (see Eq. ( A7)).

( )

In Eqs. ( A13) and (A14), the contour integrals cancel out because the perturbations δ u and δω are null along the boundaries ( ω and u are not defined straight at the boundary of the modeled domain, they do not appear in the boundary conditions), and incidentally, the adjoint state µ , which is defined up to the addition of a constant, can be considered as null along a Dirichlet-type boundary.

Reference


Amelia Domínguez-Celorrio

Manuel Vilas-Varela

Diego Peña

Véronique Langlais

Symmetry driven formation of chiral boroxine-based organometallic oligomers on Ag(001)

Two dimensional molecular assemblies on surfaces have opened a way to design and control chirality, featuring promising electronic and chemical properties that depend on the local handedness of the layer. Yet, the mechanisms leading to spontaneous chiral resolution are not fully understood at every reaction stage. Here, starting from achiral 10-bromoanthracene-9boronic acid as molecular precursor, we demonstrate enantiomeric induction in the products during the stage of covalent bonding, stemming from the competing point symmetries of the Ag(001) supporting surface and the reaction products. Upon dehalogenation and dehydration of precursors, hexagonal boroxine rings linked by organometallic anthracene dimers form, which undergo a strong interaction with four-fold symmetric substrate. The pro-chiral structural character of the resulting oligomer and its impact over the spatial distribution of the electronic molecular states are revealed by high resolution scanning tunneling microscopy and spectroscopy.

INTRODUCTION

Chiral crystallization of a racemic or achiral set of precursors is a striking phenomenon.

Still, nature tends to depart spontaneously from random distributions of enantiomers. [START_REF] Avalos | Symmetry Breaking by Spontaneous Crystallization -Is It the Most Plausible Source of Terrestrial Handedness We Have Long Been Looking for? -A Reappraisal[END_REF][START_REF] Green | The Macromolecular Route to Chiral Amplification[END_REF] A necessary condition for chirality is the breaking of inversion symmetry, which can be achieved by carrying out the synthesis atop a surface. For instance, the surface of water droplets is thought to polymerize 2D crystallites of homochiral oligopeptides, a key ingredient of the prebiotic soup. [START_REF] Weissbuch | Chiral Amplification of Oligopeptides via Polymerization in Two-Dimensional Crystallites on Water[END_REF] The process of biomineralization, which often leads to chiral macrostructures, is explained in terms of asymmetric biomolecules interacting with the step edges of inorganic crystalline surfaces. [START_REF] Orme | Formation of Chiral Morphologies through Selective Binding of Amino Acids to Calcite Surface Steps[END_REF] Another advantage of using surfaces is the confinement of reactants to a 2D space, which lowers their degrees of freedom and allows minute energy differences between enantiomorphic geometries to induce spontaneous chiral resolution [START_REF] Karageorgaki | Chiral Surface from Achiral Ingredients: Modification of Cu(110) with Phthalic Acid[END_REF][START_REF] Forster | Tailoring Homochirality at Surfaces: Going Beyond Molecular Handedness[END_REF][START_REF] Mugarza | Orbital Specific Chirality and Homochiral Self-Assembly of Achiral Molecules Induced by Charge Transfer and Spontaneous Symmetry Breaking[END_REF] and chiral amplification. [START_REF] Fasel | Amplification of Chirality in Two-Dimensional Enantiomorphous Lattices[END_REF] However, the vast majority of experimental studies about homochiral segregation on surfaces at the single molecule level is restricted to molecular self-assemblies where the building units are not covalently linked to each other. [START_REF] Zaera | Chirality in Adsorption on Solid Surfaces[END_REF][START_REF] Raval | Molecular Assembly at Surfaces: Progress and Challenges[END_REF][11][12] At the level of chemical bonding, understanding the role of a surface in chiral resolution requires the use of a well-controlled 2D reaction environment.

For this purpose, the use of on-surface synthesis (OSS) 13 techniques is a very convenient strategy.

OSS is a powerful synthetic approach in vacuum where the supporting substrate catalyzes several steps and at the same time is used to steer the reaction towards low dimensional products confined to the surface plane. Molecular precursors with specific chemical functionality diffuse throughout the surface to form designated covalent molecular nanostructures. In this process, the constrains imposed by the molecule-surface and molecule-molecule interactions play a critical role. In the last decade OSS has shown its potential to produce organic macromolecules and graphene nanoribbons 14,[START_REF] Grill | Large Molecules on Surfaces: Deposition and Intramolecular STM Manipulation by Directional Forces[END_REF] with atomically precise control over both the internal structure 16,[START_REF] Ruffieux | On-Surface Synthesis of Graphene Nanoribbons with Zigzag Edge Topology[END_REF] and the edges. [18][19][START_REF] De Oteyza | Substrate-Independent Growth of Atomically Precise Chiral Graphene Nanoribbons[END_REF] Chiral precursors can naturally form planar prochiral (formally achiral molecules that become chiral by the desymmetrization introduced by the adsorption surface) molecular structures via OSS. [START_REF] De Oteyza | Substrate-Independent Growth of Atomically Precise Chiral Graphene Nanoribbons[END_REF][21][START_REF] Merino-Díez | Transferring Axial Molecular Chirality through a Sequence of On-Surface Reactions[END_REF] Another possibility is that a heterochiral adsorption footprint [START_REF] Forster | Tailoring Homochirality at Surfaces: Going Beyond Molecular Handedness[END_REF] breaks the inversion symmetry of achiral reactants, giving rise to purely electronic chirality in certain molecular orbitals [START_REF] Mugarza | Orbital Specific Chirality and Homochiral Self-Assembly of Achiral Molecules Induced by Charge Transfer and Spontaneous Symmetry Breaking[END_REF] or even to prochiral macromolecules via selective bond activation. [START_REF] Sánchez-Sánchez | Purely Armchair or Partially Chiral: Noncontact Atomic Force Microscopy Characterization of Dibromo-Bianthryl-Based Graphene Nanoribbons Grown on Cu(111)[END_REF] Here we report the synthesis of a structurally prochiral organometallic framework (MOF) starting from achiral precursors, where the prochirality is induced by the different point symmetry of the substrate lattice and the polymer itself. Individual 10-bromoanthracene-9boronic acid (BABA) molecules were sublimated onto a clean Ag(001) surface kept at room temperature. Simultaneous debromination and dehydration of the B(OH) 2 moieties [START_REF] Schlögl | Synthesis of Two-Dimensional Phenylene-Boroxine Networks through in Vacuo Condensation and on-Surface Radical Addition[END_REF][START_REF] Clair | Growth of Boronic Acid Based Two-Dimensional Covalent Networks on a Metal Surface under Ultrahigh Vacuum[END_REF] results in a polymeric network consisting in boroxine hexagons (B 3 O 3 ) interconnected by C-Ag-C organometallic (OM) bonds 26 (see Fig. 1). The chirality appears as a consequence of the competition between the threefold symmetry of the B 3 O 3 rings and the interaction of the Ag atoms with the four-fold symmetric substrate. This competition invariably leads to enantiomeric building units that undergo chiral recognition at the OM bonds. The proposed mechanism broadens the pool of chemical pathways to obtain chiral frameworks, and could be applied to the bottom-up synthesis of 2D chiral polymers with responsive behavior to chirality, such as aminoacid sensors, filters, or seeds to transfer the chirality into the third dimension for asymmetric catalysis.

EXPERIMENTAL METHODS

Experiments were performed with an ultrahigh vacuum low temperature scanning tunneling microscope, the Joule-Thomson scanning tunneling microscopy (JT-STM) system (SPECS), with a pressure better than 1 × 10 -10 mbar. The current (I t ) is measured in the tip, which is grounded, and we use sample bias (V s ) convention. The atomically clean Ag(001) surface was prepared by cycles of ionized argon sputtering at 1kV, 3 10 -6 mbar and subsequent annealing at 450˚C. The BABA molecules were sublimated at a temperature of 200 ºC from a homemade resistive evaporator onto Ag(001) surface held at room temperature, regulating the deposition flux to 0.1-0.4 ML/min by means of a quartz microbalance. After deposition, the samples were transferred to a SPECS JT-STM operating at 4.2K. Differential conductance dI/dV data were acquired using a lock-in amplifier at 913Hz and modulation voltage of 8.8 mV, in constant height mode for the maps and constant current mode for point spectra. All images were analyzed using WSxM. [START_REF] Horcas | A Software for Scanning Probe Microscopy and a Tool for Nanotechnology[END_REF] The X-ray photoelectron spectra (XPS) of the precursor powder were measured on a XPS Spectrometer Kratos AXIS with an monochromatic incident radiation from a Al Kα source (Supplemental Figure S1).

RESULTS AND DISCUSSION

BABA molecules were synthesized following ref. [START_REF] Filthaus | Supramolecular Structures and Spontaneous Resolution: The Case of Ortho-Substituted Phenylboronic Acids[END_REF] and then sublimated onto a Ag(001) surface held at room temperature. Fig. 1a shows a STM image of a representative sample with two distinguishable phases. The chemical bifunctionality of BABA precursor including a bromine and a boronic acid moiety per molecule, opens the door to two different chemical reactions: Ullmann coupling and dehydration of boronic acids. Correspondingly, we find on one hand organometallic dimers (OD phase) that self-assemble into crystalline islands, and on the other hand, nested organometallic structures of linked boroxines with different lengths (MOF phase). The presence of these different phases suggests that there are two coexisting reaction paths (see Fig. 1b). Path A yields the formation of the organo-silver dimer, which is a wellknown intermediate of the Ullmann coupling on Ag surfaces, [START_REF] Simonov | Synthesis of Armchair Graphene Nanoribbons from the 10,10′-Dibromo-9,9′-Bianthracene Molecules on Ag(111): The Role of Organometallic Intermediates[END_REF][START_REF] Eichhorn | On-Surface Ullmann Polymerization via Intermediate Organometallic Networks on Ag(111)[END_REF] and typically precedes the covalent C-C bonding after the reductive elimination step. In the OD phase the boronic acid is cleaved. Path B takes place when three BABA precursors or organo-silver dimers retaining the boronic acid group undergo a dehydration reaction, forming a 2D structure of boroxine rings 30,31 linked by organometallic dimers. Deposition onto samples held at low temperatures could shed further insight into the reaction order of dehydration and debromination. As illustrated by Fig. 1a, in the MOF phase we have found oligomers of up to 14 boroxine units linked by organometallic dimers. ensures that the condensation into boroxines observed on Ag(001) does not take place in the crucible or during sublimation, as occurs with smaller boronic acids. [START_REF] Schlögl | Synthesis of Two-Dimensional Phenylene-Boroxine Networks through in Vacuo Condensation and on-Surface Radical Addition[END_REF] In order to check whether the B-C bonds in the pristine BABA powder are broken while heating up the crucible, we performed X-ray photoemission spectroscopy. We find no changes in either the binding energy of the B1s edge or the relative boron/carbon ratio (see Supplementary Fig. S1) after repeated sublimations of the powder. Hence, we conclude that most of the BABA precursors reach the Ag(001) surface intact. Under these conditions, we could expect that all precursors dehydrate on surface to form a boroxine, but instead we observe anthryl units bonded to boroxines (38.2±1% as derived from statistical analysis of 6100 anthracenes) and no intact acid groups at all, which points towards a strong tendency to the C-B cleavage of the B(OH) 2 groups on coinage metal surfaces (see data on Ag(111) and Au(111) in Supplementary Note). This finding could be applied to implement Suzuki cross-coupling [START_REF] Miyaura | Stereoselective Synthesis of Arylated (E)-Alkenes by the Reaction of Alk-1-Enylboranes with Aryl Halides in the Presence of Palladium Catalyst[END_REF] in OSS using similar precursors to produce graphene nanoribbons. In our case, the smooth anthracene edges at the position previously occupied by boron atoms (Fig. 2 and Supplemental Fig. S3), suggest that the C bonds are passivated by hydrogen after cleavage (for example coming from the OH-groups). For BABA on Ag(001), no single anthracene with intact boronic acid was found on the surface, suggesting that the formation of MOF strongly depends on the availability of pristine precursors to produce the dehydration of boronic acids before bonding to a silver adatom or losing the acid. On the other hand, 94±0.1% of the boroxine terminal anthracenes are part of an organo-silver dimer, evidencing that the condensation of boronic acids into a boroxine does not affect the reactivity of brominated carbons.

The dehydration of aryl boronic acid derivatives with halogen groups has also been observed on Ag(111). [START_REF] Schlögl | Synthesis of Two-Dimensional Phenylene-Boroxine Networks through in Vacuo Condensation and on-Surface Radical Addition[END_REF]34 However, in the case of Ag(001) the organo-silver dimers obtained at room temperature do not transform into C-C bonds, as is characteristic of the Ullmann coupling.

Subsequent post-annealing does not help to initiate this second step because desorption takes place at 300 ºC, before any substantial change of the MOF is observed (see Supplemental Figure S4). The high stability of the C-Ag bonds in both phases suggests that the steric repulsion between the hydrogens in vicinal rings of our planarly adsorbed precursors impedes the C-C bonding to form bianthracene moieties. The MOF geometry is strongly influenced by the interaction of the substrate with the remaining Ag atoms between anthracenes. Fig. 2a shows a high-resolution image of a short MOF with three linked boroxines. The Ag atoms appear as round protrusions surrounded by two elongated anthryl units, which are arranged around the boroxine centers and exhibit bonding axes at approximately 120º. The seven Ag atoms contained in the structure do not have the same appearance. The two Ag atoms bonded to anthracenes along the [1-10] direction of the substrate lattice (Fig. 2b) are clearly brighter than the others. This is an electronic fingerprint for different adsorption sites. Because of the higher symmetry of the dimers along {110} directions, we assume that they correspond to the energetically most favorable adsorption geometry.

Considering that the sp 2 boron orbitals are oriented with an angle of 120º between each other, the bonding axes of the two other branches of the boroxine should exhibit a threefold symmetric arrangement. However, this is prohibited by the fourfold symmetry of the substrate they are strongly interacting with. The arrangement of the other two branches will have to minimize at the same time the adsorption energy of the remaining Ag adatoms, and the energy penalty due to the distortion of the C-B-O bonding angles away from 120º. The independent optimization of both terms is not allowed, forcing in this way a compromise for the adsorption site of the Ag atoms between the fourfold and the threefold symmetries. Thus, attending to the sign of the angle formed by the actual B-Ag segment and the axes at 120 º from the [110] direction (corresponding to the ideal hexagonal structure, white dotted lines in Fig. 2b), each building block composed by one B 3 O 3 ring and its three organo-silver dimers can be classified into a chiral rotational sense.

For example, the two buildings blocks found in Fig 2 exhibit R-handedness.

The minimal structural model that allows us to fit our high-resolution images consists in three lines pinned at the boron atoms of a perfect hexagonal boroxine. The Ag atoms are placed at the middle point of the lines. The only free parameters are defined in Figure 2a: (i) the relative angles between the dimer axis along {110} and the axes of the other two branches (, ); (ii) the B-Ag distance along the other two branches (d 2 and d 3 ). The length of organo-silver dimers is kept constant to a value of 0.76(9) nm as extracted from the analysis of individual dimers in the OD phase (Supplemental Figure S3b). The size of the boroxine hexagon is set to the value of one benzene ring. Finally, the inter-boron distance along {110} directions (2d 1 ) is also fixed in our model to 1.32(1) nm. This value was calculated as the measured inter-boroxine distance along {110} directions minus the length of the boroxine hexagon. Note that this model is equivalent to fixing the B-C bond length and C-B-O angle along {110} directions to their experimental values and letting free the other two. This model is fitted to STM constant-height current images of different sample areas where creep and thermal drift are relaxed below a negligible level of less than 1 pm during the scan of a MOF unit as the one shown in Fig. 2a.

In Fig. 2a we superimpose the fit of our structural model to one building block. The comparison clearly shows that ≠≠120º. Since high symmetry positions cannot be occupied by the Ag atoms in branches 2 and 3, several structural models could exist, each one corresponding to a metastable adsorption geometry (see for instance the Ag sites in the models sketched in Figure 2a). We have identified three of them (as well as their corresponding enantiomers, see an example at Supplementary Fig. S5) by analyzing a total of 86 boroxine building blocks in high resolution images of the same quality as Fig. 2a. The structural parameters for these models are given in Table 1. The building blocks are always prochiral, with deviations ranging 0.5(5) º<|-120|<8.5(5) º. In general, d 1 , d 2 and d 3 are different. Because our model assumption that the size of the organo-silver dimer yields a perfect fit of the images (see sketch in Figure 3), the different Ag-B distances imply that C-B bond length depends on the linking direction. As shown in Figure 3, MOFs self-assemble into homochiral domains of the order of 10 nm. The MOF phase island in this example is in contact with the OD phase, which is a general trend in our samples (see Fig. 1a). Here, we identify up to three levels of chiral recognition. First, we note that the C-Ag-C bonding direction in the branches with low adsorption symmetry is univocally determined by β and  (see Fig. 2a). Henceforth, these angles fix the growth directions in which a MOF can bind to another building block via OM bonds. The link only occurs when  or  in the reacting building blocks coincide, which permits a straight organometallic bond. This precludes reactions between building blocks of different chirality (see the large angle difference between enantiomers in Table 1) along directions away from {110}.

Trimer

This enantioselectivity is therefore imposed by geometric constrains. Second, in the case of C-Ag-C bonds along {110} directions, we also find asymmetric synthesis in all MOFs. In this direction, however, there is not any geometric burden to link building blocks with different chirality, because anthracenes and Ag atoms have a unique adsorption geometry that is totally independent from  and . This provides strong evidence that the enantioselectivity in the MOF growth is not only due to geometric constrains, but rather has an underlying electronic component, say, a formation energy that depends on the chiral character of the molecular orbitals and favors the OM bonding between boroxine building blocks with the same chiral character.

Although ab initio calculations would be necessary to understand the interaction mechanism, in our work this electronic effect is supported by the energy resolved spectroscopy discussed hereinafter. In addition, a similar effect has been previously reported for weak intermolecular electronic interactions driving homochiral self-assembly of Cu-phtalocyanines. 7 1) are overlaid on prochiral (R)organometallic polymer fragments. Blue lines represent dimers linking boroxines in directions different from {110}. In these cases, since 2d 2 does not match with full precision the particular inter-boroxine distance, the Ag-B distance adapts with some tensile strain of about 5 %.

Finally, we observe domains where all individual MOFs have the same handedness (Figure 3a).

This spontaneous segregation in homochiral domains is attributed to halogen-mediated selfassembly, previously reported for other molecules on noble metal substrates. [START_REF] Fasel | Amplification of Chirality in Two-Dimensional Enantiomorphous Lattices[END_REF]35 In our case, Br atoms are clearly imaged as round protrusions between MOFs, occupying equivalent positions relative to the building blocks in a regular pattern. This suggests a strong interaction between H terminated carbons and the free halogen atoms left over the surface after cleaving the Br-C bonds. This interaction, of electrostatic origin, may play a relevant role during the formation building blocks, favoring the availability of neighboring units with the same chirality. show some weaker intensity distributed across the OM dimers along the [110] axis. This delocalization is not observed across the OM dimers lying along low symmetry directions. On the other hand, the orbital at 2.45 V (Figure 4d) has a triangular shape with sides coincident with the anthracenes bound to a boroxine ring. The corners of these triangles are approximately located at the terminal hydrogens of the anthracenes, and exhibit different intensities, which could be attributed to some tilting with respect to the surface plane.

Additional orbitals pop up in the constant height maps at 1.5 eV (Figure 4e), which are not visible in point spectra due to the intrinsic Ag(001) contribution at this energy range (grey curve in Figure 4b). These orbitals can be described as three lobe structures around the boroxines, being reminiscent of the 2.45 eV molecular state. Two lobes lie just outside the phenyl rings, and another one on top of a B-C bond, all three belonging to the chiral branches of the building block. In this state, the previously observed chirality in the structural analysis becomes evident in the electronic structure too. The 1.5 V state overlaps in space with the boroxine state, which in turn has been found to connect boroxines along {110} directions. This indirect overlap of the chiral orbitals through the boroxine orbitals (which does not discard the existence of additional chiral selective chemical interactions which are not visualized in our spectroscopy data), provides a rationale for the enantiospecific links along these directions in the absence of geometric constrains. The 1.5 V map displays a very faint intensity over the Br atoms as well.

CONCLUSIONS

In conclusion, the competition between the internal symmetries of an on-surface synthesized macromolecule and the supporting substrate provides a mechanism to confer chirality on the otherwise achiral product. Chiral recognition during the MOF growth is ascribed to geometric constrains and to the sensitivity of orbital overlap to the chiral character of the basic building blocks. We also highlight the role of halogen bonds in the homochiral segregation of selfassembled MOF islands. This work shows a strategy to create chiral organometallic networks from achiral precursors directly on the surface. We also evidence how Ullmann coupling intermediates can adopt chiral configurations depending on their interaction strength with surface provided that there are different point symmetries.

Figure 1 .

 1 Figure 1. (a) Constant current STM image (V s = 0.5 V; I t = 0.05 nA) of a representative BABA

Figure 2 .

 2 Figure 2. (a) STM topography of atomically resolved Ag(001) surface (Vs = 20 mV; It = 50 nA)

  parameters have been extracted measuring on high resolution STM images. Angles are determined with an error of 0.5 º. Model 1 and 2 are sketched in Figure2a(red and wine respectively), model 3 is sketched in Figure3(pink).

Figure 3 .

 3 Figure 3. Constant height STM image of a self-assembled MOF island (V s = 5 mV). The

Figure 4 .

 4 Figure 4. (a) Constant height STM image of MOF with (R)-Model 1 and (R)-Model 2 as

Table 1 .

 1 Structural parameters of prochiral building blocks of MOF. All the models and the

		Handedness	β (º)	γ (º)	d 2 /d 1	d 3 /d 1
	Model 1	R	117	120.5 1.21(3)	1.22(5)
	Model 1	S	120.5	117	1.22(5)	1.21(3)
	Model 2	R	115	128.5 1.20(6)	1.05(1)
	Model 2	S	128.5	115	1.05(1)	1.20(6)
	Model 3	R	115	122.5 1.23(2)	0.98(1)
	Model 3	S	122.5	115	0.98(1)	1.23(2)

SUPPORTING INFORMATION DESCRIPTION

The following files are available free of charge as Supporting information.

-babaAg001_SI_v3.0 (PDF file): Peripherial data on Au(111) and Ag(111) and comments including XPS precursor powder characterization and image sequence illustrating the effect of annealing no Ag(001)

AUTHOR INFORMATION

The authors declare no competing financial interests.

ACKNOWLEDGMENT

Financial support was provided by the European Union through Interreg-POCTEFA (Grant TNSI/EFA194/16) and Horizon 2020 (Grant Agreement No. 863098 ("SPRING")) programs, including EFRD funds, as well as by the Spanish Plan Nacional de I+D+i (MAT2016-78293-C6 and PID2019-107338RB-C6). We also acknowledge funding from Gobierno de Aragón (grupo de referencia E13_20R "NANOMIDAS") and Xunta de Galicia (Centro singular de investigación de Galicia, accreditation 2016-2019, ED431G/09).