Corinna Köpke

James Irving
email: james.irving@unil.ch

Delphine Roubinet

Stochastic inversion for soil hydraulic parameters in the presence of model error: an example involving ground-penetrating radar monitoring of infiltration

Proxy forward solvers are commonly used in Bayesian solutions to inverse problems in hydrology and geophysics in order to make sampling of the posterior distribution, for example using Markov-chain-Monte-Carlo (MCMC) methods, computationally tractable. However, use of these solvers introduces model error into the problem, which can lead to strongly biased and overconfident parameter estimates if left uncorrected. Focusing on the specific example of estimating unsaturated hydraulic parameters in a layered soil from time-lapse ground-penetrating radar data acquired during a synthetic infiltration experiment, we show how principal component analysis, conducted on a suite of stochastic model-error realizations, can for some problems be used to build a sparse orthogonal basis for the model error arising from known forward solver approximations and/or simplifications. Projection of the residual onto this basis during MCMC permits identification and removal of the model error before calculation of the likelihood. Our results indicate that, when combined with an informal likelihood metric based on the expected behaviour of the 2 -norm of the residual, this methodology can yield posterior parameter estimates exhibiting a marked reduction in bias and overconfidence when compared to those obtained with no model-error correction, at reasonable computational cost. 1

Introduction

Stochastic parameter estimation and inversion have become increasingly popular in hydrology and geophysics over the past decade. In particular, it is now computationally feasible and common to solve many inverse problems in these domains in a Bayesian manner, whereby prior knowledge about the subsurface parameters of interest is combined with measured data to yield a posterior probability distribution. The latter is typically sampled using Markov-chain-Monte-Carlo (MCMC) methods [START_REF] Linde | On uncertainty quantification in hydrogeology and hydrogeophysics[END_REF]. Notable advantages of the Bayesian-MCMC approach are that (i) it is highly flexible and can incorporate any information that can be expressed as a probability density into the inverse problem; (ii) it provides a natural framework for data integration; and (iii) it has the potential to provide more accurate parameter uncertainty estimates than traditional methods based on linearization. This does, however, come at the cost of being highly computationally expensive. Indeed, many thousands if not millions of MCMC iterations, each requiring a numerical solution of the forward problem, are typically required to obtain a sufficient number of posterior samples for use in subsequent probabilistic forecasting and risk analysis (e.g., [START_REF] Ruggeri | Systematic evaluation of sequential geostatistical resampling within MCMC for posterior sampling of near-surface geophysical inverse problems[END_REF].

A critical component of framing an inverse problem in a Bayesian context is proper characterization of the expected statistical nature of the residual. This is, in order to formulate the likelihood, we must have detailed knowledge about the statistical distribution of the difference between the measured data and those calculated through the numerical solution of the forward problem on the "true" set of subsurface model parameters. In arguably most cases, the residual is attributed solely to data-measurement errors and described as multi-Gaussian, usually with independent and identically distributed elements (e.g., [START_REF] Bodin | Seismic tomography with the reversible jump algorithm[END_REF][START_REF] Gallagher | Markov chain Monte Carlo (MCMC) sampling methods to determine optimal models, model resolution and model choice for Earth Science problems[END_REF][START_REF] Irving | Stochastic inversion of tracer test and electrical geophysical data to estimate hydraulic conductivities[END_REF][START_REF] Linde | Distributed soil moisture from crosshole ground-penetrating radar travel times using stochastic inversion[END_REF][START_REF] Scholer | Examining the information content of time-lapse crosshole GPR data collected under different infiltration conditions to estimate unsaturated soil hydraulic properties[END_REF][START_REF] Vrugt | Treatment of input uncertainty in hydrologic modeling: Doing hydrology backward with Markov chain Monte Carlo simulation[END_REF]. This is despite the fact that, in order to improve the computational tractability of the Bayesian-MCMC approach, approximate versions of the forward solver, for example using coarsened discretizations and/or simplifications of the underlying physics, are typically employed (e.g., [START_REF] Christen | Markov chain Monte Carlo using an approximation[END_REF][START_REF] Cui | Bayesian calibration of a large-scale geothermal reservoir model by a new adaptive delayed acceptance Metropolis Hastings algorithm[END_REF][START_REF] Efendiev | Modified Markov Chain Monte Carlo method for dynamic data integration using streamline approach[END_REF][START_REF] Hinnell | Improved extraction of hydrologic information from geophysical data through coupled hydrogeophysical inversion[END_REF][START_REF] Ray | Bayesian calibration of the community land model using surrogates[END_REF][START_REF] Scholer | Examining the information content of time-lapse crosshole GPR data collected under different infiltration conditions to estimate unsaturated soil hydraulic properties[END_REF]. The use of such computationally efficient "proxy solvers" leads to model error, which if left uncorrected has the potential to overwhelm the effects of data measurement uncertainties and lead to strongly biased and overconfident posterior distributions [START_REF] Brynjarsdóttir | Learning about physical parameters: The importance of model discrepancy[END_REF].

In recent years, a number of studies in the hydrological and geophysical literature have attempted to address the issue of model error in Bayesian inversions, with the aim of making more effective use of proxy solvers when dealing with computationally expensive forward operators. In general, the approaches that have been presented can be divided into two categories.

In the first category, researchers have focused on the overall or "global" statistical characterization of the model error, with the goal of using this information to develop more appropriate parametric likelihood functions that better reflect the true nature of the residual. This has generally been accomplished through the analysis of stochastic model-error realizations, which are generated by running the full and approximate forward solvers for randomly drawn sets of model parameters. Typically, multi-Gaussian statistics for the model error are assumed, meaning that means and covariances estimated from the realizations can be incorporated into a Gaussian likelihood (e.g., [START_REF] Arridge | Approximation errors and model reduction with an application in optical diffusion tomography[END_REF][START_REF] Hansen | Accounting for imperfect forward modeling in geophysical inverse problems -exemplified for cross hole tomography[END_REF][START_REF] Kaipio | Statistical inverse problems: Discretization, model reduction and inverse crimes[END_REF][START_REF] Lehikoinen | Dynamic inversion for hydrological process monitoring with electrical resistance tomography under model uncertainties[END_REF][START_REF] Stephen | Scale and process dependent model errors in seismic history matching[END_REF], but other parametric likelihood functions have also been considered (e.g., Del [START_REF] Del Giudice | Improving uncertainty estimation in urban hydrological modeling by statistically describing bias[END_REF][START_REF] Schoups | A formal likelihood function for parameter and predictive inference of hydrologic models with correlated, heteroscedastic, and non-Gaussian errors[END_REF][START_REF] Smith | Development of a formal likelihood function for improved Bayesian inference of ephemeral catchments[END_REF][START_REF] Smith | Modeling residual hydrologic errors with Bayesian inference[END_REF]. Accounting for model error in this manner has been shown to lead to broadened posterior distributions and a reduction in parameter bias. One key issue, however, concerns the validity of the assumption that the errors can be adequately described by a given parametric distribution. In many inverse problems in geophysics and hydrology, for example, model errors will exhibit complex statistics and correlations that arise from the typically high dimension of the data and/or model-parameter spaces in these problems, combined with the non-linearity of the forward operators involved. Indeed, there has been much increased interest in "likelihood-free" inference methods such as approximate Bayesian computation (ABC) (e.g., [START_REF] Vrugt | Towards diagnostic model calibration and evaluation: Approximate Bayesian computation[END_REF] and generalized likelihood uncertainty estimation (GLUE) (e.g., [START_REF] Beven | The future of distributed models: Model calibration and uncertainty prediction[END_REF], to a large extent because of this issue.

In the second category of developed approaches for addressing model error, researchers have focused on building a parameter-dependent or "local" error model in order to describe the discrepancy between the full and approximate forward solvers. As with the approaches mentioned above, this is constructed based on computed realizations of the model error for different parameter sets. However, in this case the results are used to effectively correct the output of the approximate solver rather than to develop a more appropriate Bayesian likelihood function. Construction of the error model can be done in a number of different ways. This includes simple nearest-neighbour or linear interpolation between model-error realizations (e.g., [START_REF] Cui | Bayesian calibration of a large-scale geothermal reservoir model by a new adaptive delayed acceptance Metropolis Hastings algorithm[END_REF][START_REF] O'sullivan | Simulation error models for improved reservoir prediction[END_REF], representing the discrepancy as a Gaussian process conditioned to the points in the parameter space where the model error is known (e.g., [START_REF] Kennedy | Bayesian calibration of computer models[END_REF]Xu and Valocchi , 2015), or using statistical regression approaches (e.g., [START_REF] Doherty | Use of paired simple and complex models to reduce predictive bias and quantify uncertainty[END_REF][START_REF] Josset | Functional error modeling for uncertainty quantification in hydrogeology[END_REF]. In all of this work, the implicit assumption is that the full and approximate model-response surfaces are regular enough such that the model error for a set of parameter values where it is unknown can be effectively predicted through some kind of interpolation between the existing realizations.

While this may be the case for some inverse problems, difficulties can arise in the presence of strongly non-linear forward solvers and/or large numbers of model parameters. That is, it may not be possible to sufficiently sample the model-parameter space with model-error realizations such that interpolation between these realizations will provide a reliable model-error estimate at some new location.

Recently, [START_REF] Köpke | Accounting for model error in Bayesian solutions to hydrogeophysical inverse problems using a local basis approach[END_REF] presented a new approach to account for the model error arising from the use of proxy forward solvers in Bayesian-MCMC inversions, whereby information about the error is gathered during the inversion procedure through occasional runs of the approximate and full solvers together, the results of which are stored in a dictionary. In contrast to the existing methods mentioned above, the approach of [START_REF] Köpke | Accounting for model error in Bayesian solutions to hydrogeophysical inverse problems using a local basis approach[END_REF] focuses on the projection-based identification of the model-error component of the residual through the construction of a local, parameter-dependent, orthogonal model-error basis, rather than on attempting to fit the overall model-error statistics to a prescribed statistical distribution or develop an interpolation-based error model. The model error estimated by projecting onto the basis is then subtracted from the residual before computing the likelihood of a proposed set of model parameters in MCMC. Application of the approach of [START_REF] Köpke | Accounting for model error in Bayesian solutions to hydrogeophysical inverse problems using a local basis approach[END_REF] to a high-dimensional spatially distributed tomographic example was found to yield parame-ter estimates exhibiting a notable reduction in bias compared to those obtained when the model error was ignored. The presented method does, however, still require occasional runs of the full forward solver along the Markov chain as MCMC iterations progress, which can be computationally costly depending on the problem at hand.

In this paper, we build on the work of [START_REF] Köpke | Accounting for model error in Bayesian solutions to hydrogeophysical inverse problems using a local basis approach[END_REF] and show that, for some inverse problems, it may be possible to derive a suitable global basis for the model error over the entire parameter space through the application of principal component analysis to a large number of stochastic model-error realizations. These realizations can be conveniently computed in parallel prior to MCMC, and they must be properly organized in the data space before analysis to maximize similarity in their spatial characteristics. We begin in Section 2 with an overview of Bayesian-MCMC inversion, followed by a detailed description of our developed approach. This is followed in Section 3 with application to synthetic data corresponding to a vadosezone inverse problem that has been the subject of much investigation in previous work, which is the estimation of unsaturated soil hydraulic parameters from time-lapse zero-offset-profile (ZOP) ground-penetrating radar (GPR) data acquired during infiltration. We compare the results obtained with our methodology to those obtained when no model-error correction is applied, and conclude in Section 4 with an overall assessment of the method with regard to its advantages and limitations.

Methodological background

Bayesian-MCMC inversion

Consider to begin the forward problem linking a set of M subsurface model parameters of interest m true ∈ R M to a set of N measured or observed data d obs ∈ R N :

d obs = F (m true ) + e d , (1)
where forward operator F : R M → R N contains the physics and geometry of the measurements and e d is vector of data measurement errors. The goal of the corresponding inverse problem is to estimate m true given d obs , which requires knowledge of F and in most cases some prior information about the model parameters. From a probabilistic point of view, this can be formulated using Bayes' theorem, whereby an initial prior state of information for the model parameters ρ(m) is updated to a more refined posterior state of knowledge σ(m) based on the available data (e.g., [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF]. That is,

σ(m) = k L(m) ρ(m), (2)
where normalization constant k ensures that the posterior probability distribution integrates to unity, and likelihood L(m) expresses the conditional probability of model parameter set m

given the observed data d obs . Assuming that (i) the underlying physics are completely known and considered in the inverse problem; and (ii) the data measurement errors are independent and identically normally distributed having mean zero and standard deviation s d , L(m) takes on the simple multi-Gaussian form

L(m) = 1 (2πs 2 d ) N/2 exp - ||r(m)|| 2 2s 2 d , (3)
where || • || denotes the 2 -norm and r(m) is the residual or difference between the observed data and those predicted for some model parameter set m using F . The latter quantity is given by

r(m) = d pred -d obs = F (m) -[F (m true ) parameter-error component + e d ], (4)
where d pred denotes the predicted data. We see in equation (3) that L(m) will be maximized when the 2 -norm of r(m) is minimized, which corresponds to the case where m = m true and the parameter-error component as defined in equation ( 4) is zero. The spread of the likelihood distribution about the maximum value is controlled by the data measurement error standard deviation s d along with the number of data N , with larger errors and lesser amounts of data yielding broader likelihoods.

Equations ( 2) through (4) together provide a means of calculating the posterior probability of a set of model parameters m. This is commonly used within MCMC sampling algorithms to quantify posterior uncertainty, thereby solving the inverse problem, since it is not generally possible to perform the multi-dimensional integrations needed to obtain the statistical moments of σ(m). In this regard, a basic Metropolis-Hastings algorithm [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] that is guaranteed (after burn-in) to generate a Markov chain of samples {m 1 , ..., m k } from the Bayesian posterior distribution proceeds as follows:

1. Draw the first model in the Markov chain m 1 from the Bayesian prior distribution ρ(m).

Set i = 1.

2. Draw a perturbed model-parameter set m from the proposal distribution Q(m |m i ), whose width around m i is chosen so as to provide a balance between efficiently moving through the parameter space and generating proposals that have a reasonable probability of being accepted.

3. Calculate the probability of accepting m as the next model in the Markov chain using

P acc = min 1, σ(m ) Q(m i |m ) σ(m i ) Q(m |m i ) . (5)
4. Draw a random number x ∈ U (0, 1). If x ≤ P acc , then set m i+1 = m . Otherwise set

m i+1 = m i .
5. Set i = i + 1 and go to Step 2.

Note that, in the case where the proposal distribution is symmetric (i.e., Q(m

|m i ) = Q(m i |m )),
the above algorithm reduces to the original MCMC sampler of [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] where the acceptance probability is given by

P acc = min[1, σ(m )/σ(m i )]
. We consider the latter sampler for the example inversions presented in Section 3.

Accounting for model error

Likelihood equation ( 3) is perfectly theoretically valid for the case where the only contribution to the difference between the observed and predicted data, when considering the correct set of model parameters m true , is a set of Gaussian data-measurement errors having standard deviation s d . However, as mentioned previously, approximate forward solvers are typically used in hydrological and geophysical problems to improve the computational efficiency of the Bayesian-MCMC procedure, meaning that the residual more realistically takes the form

r(m) = F (m) -[F (m true ) + e d ] = F (m) -F (m) model-error component + F (m) -[F (m true ) parameter-error component + e d ], (6)
where F is the approximate forward operator. The presence of an additional model-error term in equation ( 6) as compared with equation (4), which is commonly of large magnitude, strongly correlated, and/or highly non-Gaussian [START_REF] Kaipio | Statistical inverse problems: Discretization, model reduction and inverse crimes[END_REF][START_REF] Schoups | A formal likelihood function for parameter and predictive inference of hydrologic models with correlated, heteroscedastic, and non-Gaussian errors[END_REF][START_REF] Smith | Development of a formal likelihood function for improved Bayesian inference of ephemeral catchments[END_REF], makes use of likelihood expression (3) inappropriate. In particular, it means that (i) the residual will not necessarily be minimized when m = m true , implying posterior parameter bias; and (ii) feasible model parameter sets may have an extremely low likelihood when considering realistic levels of data error. Although simple inflation of s d can be used to broaden the Gaussian likelihood and reduce the latter issue, it cannot address the former and be viewed as an effective solution for reliable posterior uncertainty quantification.

In order to address the model-error issue, we build on the work of [START_REF] Köpke | Accounting for model error in Bayesian solutions to hydrogeophysical inverse problems using a local basis approach[END_REF] in this paper and focus on learning about the nature of the model error through stochastic simulation such that it may be identified and removed from the residual during MCMC. The overall idea is that, for some problems, a representative set of stochastic model-error realizations, computed prior to MCMC for random model parameter sets using the full and approximate forward solvers, can be used to construct an orthonormal basis for the model error. Projection of the residual onto this basis in each MCMC iteration is used to isolate the model-error component, which is the subtracted from r(m) before calculating the likelihood. Note that, whereas [START_REF] Köpke | Accounting for model error in Bayesian solutions to hydrogeophysical inverse problems using a local basis approach[END_REF], used a dictionary-based K-nearest-neighbour (KNN) approach to construct a different local model-error basis for each proposed set of model parameters in MCMC, with runs of the full forward solver being required periodically along the entire Markov chain, we focus here on the development of a global basis (i.e., over the entire model parameter space) before posterior sampling begins. Although not appropriate for all problems, this methodology has the advantage that all expensive forward solver computations can be run in a simple parallel manner outside of the MCMC iterations. The corresponding set of modelerror realizations can also be directly reused in any subsequent inversions. Our approach proceeds as follows:

1. Generate k random sets of model parameters {m 1 , ..., m k } from the Bayesian prior distribution ρ(m).

2. Compute the corresponding set of stochastic model-error realizations {E 1 , ..., E k }, where

E i = F (m i ) -F (m i ).
3. If necessary, organize the information in each realization to improve coherency for subsequent analysis (see Section 3.3).

R(m ) = r(m ) -B B T r(m ) (7) 6. Use R(m ) to determine L(m ) within MCMC (see Section 2.3).
It is important to note that the success of the modified MCMC approach described above, in terms of providing refined and unbiased posterior parameter estimates using an approximate forward operator, hinges on our ability to effectively separate the model-error component of equation ( 6) from (i) data-measurement errors, and (ii) parameter-related errors. The implicit assumption in our work is that these two other sources of error lie orthogonal to the elements of B, such that projection of the residual onto the basis will preserve only the modelerror component. With regard to (i), we have found that this is a reasonable expectation as the limited number of model-error basis vectors, which tend to possess a high degree of spatial correlation, are generally not capable of representing random data-measurement errors through a linear combination [START_REF] Köpke | Accounting for model error in Bayesian solutions to hydrogeophysical inverse problems using a local basis approach[END_REF]. With respect to (ii), although there is no guarantee that the basis cannot represent at least part of the parameter-error term through a linear combination, our experience has been that the model and parameter-related errors typically possess significantly different statistical characteristics meaning that the latter tend to be quite effectively attenuated through projection onto B . If this is not the case and a particular incorrect model-parameter set tested within MCMC happens to yield a parametererror component that resembles what was observed in the model-error realizations, this error will be removed and the parameter set will have a reasonably high chance of being accepted [START_REF] Köpke | Accounting for model error in Bayesian solutions to hydrogeophysical inverse problems using a local basis approach[END_REF]. This latter point is discussed in further detail in Section 3.5.

Likelihood evaluation

Ideally, the remainder R(m) in equation ( 7) should represent the residual in ( 6) with the model-error component perfectly removed, meaning that it should be identical to equation ( 4)

and thus suitable for inclusion into Gaussian expression (3) to evaluate the likelihood. In reality, however, small but correlated and non-Gaussian errors in the approximation of the model-error component of the residual, related to our inability to perfectly separate model error from data measurement and parameter uncertainty using the sparse basis B, mean that R(m true ) will deviate somewhat from multi-Gaussian and use of equation ( 3) can be problematic. Indeed, the strong ranking of models provided by a Gaussian likelihood function is well understood to pose difficulties for Bayesian inference when the underlying statistical assumptions regarding the residual are violated, in the sense that sets of model parameters that are perfectly acceptable may be mapped to extremely low likelihoods (e.g., [START_REF] Beven | The future of distributed models: Model calibration and uncertainty prediction[END_REF][START_REF] O'sullivan | Simulation error models for improved reservoir prediction[END_REF]. To address this issue, we evaluate the likelihood in this work using a statistically informal but more practical metric based on the expected univariate distribution of the 2 -norm of the remainder in equation ( 7), as opposed to the expected multivariate distribution of the vector R(m). Specifically, assuming for lack of better information that the elements of R(m) are uncorrelated and normally distributed having mean zero and standard deviation s R , it can be shown that the 2 -norm ||R(m)|| will follow a scaled chi distribution [START_REF] Forbes | Statistical Distributions, Wiley Series in Probability and Statistics -Applied Probability and Statistics Section Series[END_REF], leading to the following equation:

L(m) = 2 1-N 2 Γ N 2 , 0 s -N R ||R(m)|| N -1 exp - ||R(m)|| 2 2s 2 R , (8)
where Γ(•, •) is the incomplete gamma function.

For the typical case where s R is unknown and can only be bounded between lower and upper values s R 1 and s R 2 , respectively, equation ( 8) can be integrated over s R yielding

L(m) ∝ Γ N -1 2 , ||R(m)|| 2 2s 2 R s R 2 s R 1 . (9)
Figure 1 shows L(m) calculated using equation ( 9 9) is much less sensitive to small changes in R(m) compared with equation ( 3), and represents a significantly more relaxed and inherently conservative constraint than that provided by a formal Gaussian likelihood function. Indeed, the use of such informal likelihood measures within stochastic inverse methods has gained widespread acceptance in hydrology and other domains in recent years (e.g., [START_REF] Beven | Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the GLUE methodology[END_REF][START_REF] Beven | GLUE: 20 years on[END_REF][START_REF] Blasone | Generalized likelihood uncertainty estimation (GLUE) using adaptive Markov Chain Monte Carlo sampling[END_REF][START_REF] Nott | Generalized likelihood uncertainty estimation (GLUE) and approximate Bayesian computation: What's the connection?[END_REF]Sadegh and Vrugt, 2013a,b;[START_REF] Wilkinson | Approximate Bayesian computation (ABC) gives exact results under[END_REF], as researchers have realized the shortcomings of placing too much importance on the detailed statistical properties of the residual for many real-world problems. Equation ( 9) can in fact be considered as a slight variation of the generalized likelihood uncertainty estimation (GLUE) approach, originally proposed by [START_REF] Beven | The future of distributed models: Model calibration and uncertainty prediction[END_REF], where the distinction between "behavioural" and "non-behavioural" models is quantified by using the 2 -norm of the data misfit. Models whose remainder norm falls within the bounds indicated in Figure 1 will have a high chance of being accepted in MCMC, whereas those falling significantly outside these bounds will tend to be rejected.

[Figure 1 about here.]

3 Example: GPR monitoring of infiltration

We now apply the model-error methodology presented in Section 2 to a synthetic example involving GPR monitoring of an infiltration experiment. Zero-offset-profile (ZOP) GPR data, acquired between two boreholes over the course of the experiment, provide estimates of horizontally averaged soil water content as a function of depth and time (e.g., [START_REF] Annan | GPR methods for hydrogeological studies[END_REF]).

Together with a numerical model for the infiltration process, the latter results are then used to estimate unsaturated soil hydraulic properties in a layered subsurface. This particular problem has been the focus of much previous research in the field of hydrogeophysics (e.g., [START_REF] Binley | Vadose zone flow model uncertainty as conditioned on geophysical data[END_REF][START_REF] Cassiani | Modeling unsaturated flow in a layered formation under quasi-steady state conditions using geophysical data constraints[END_REF][START_REF] Looms | Identifying unsaturated hydraulic parameters using an integrated data fusion approach on cross-borehole geophysical data[END_REF][START_REF] Rucker | Parameter estimation for soil hydraulic properties using zero-offset borehole radar[END_REF][START_REF] Rucker | Inverse upscaling of hydraulic parameters during constant flux infiltration using borehole radar[END_REF], and was most recently investigated within the context of Bayesian-MCMC inversion by [START_REF] Scholer | Estimating vadose zone hydraulic properties using ground penetrating radar: The impact of prior information[END_REF][START_REF] Scholer | Bayesian Markov-Chain-Monte-Carlo inversion of time-lapse crosshole GPR data to characterize the vadose zone at the Arrenaes site, Denmark[END_REF][START_REF] Scholer | Examining the information content of time-lapse crosshole GPR data collected under different infiltration conditions to estimate unsaturated soil hydraulic properties[END_REF]. Here, we consider the model errors arising from a simplifying assumption common to all past work, which is that water movement occurs in a purely vertical direction through the subsurface.

Governing equations and model simplifications

The general movement of water through unsaturated soils is described by Richards' equation [START_REF] Richards | Capillary conduction of liquids through porous mediums[END_REF], given by

∂θ(h) ∂t = ∇ • [K(h)∇h] + ∂K(h) ∂z , ( 10
)
where θ is the volumetric water content, K is the unsaturated hydraulic conductivity, h is pressure head, t is time, and z is elevation. The relationships θ(h) and K(h) for different soils are commonly described using the van Genuchten -Mualem (VGM) model [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF][START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soil[END_REF]. With this model, the soil water retention, expressed in terms of effective saturation S e , is given by

S e (h) = θ(h) -θ r θ s -θ r =      (1 + |αh| n ) -m , for h ≤ 0 1 , for h > 0 (11)
where θ r and θ s are the residual and saturated water contents, respectively, and α, m, and n are empirical shape factors with m = 1 -1/n. The unsaturated hydraulic conductivity is described by

K(h) = K s S e (h) 1/2 1 -(1 -S e (h) 1/m ) m 2 , ( 12
)
where K s is the hydraulic conductivity value at full saturation. A total of five parameters (K s , θ r , α, n, and θ s ) therefore describe a soil's hydraulic properties using the VGM model.

Equations ( 10) through ( 12) provide a link between a set of subsurface VGM parameters and the corresponding spatiotemporal distribution of water content in response to infiltration.

That is, knowing the distribution of soil VGM parameters along with the boundary and initial conditions of the infiltration experiment, we can calculate the evolution of water content in the subsurface. This forward link provides the basis for inverting for the soil hydraulic properties given a set of dynamic GPR-derived water-content measurements. However, in the context of stochastic inversion, repeated solution of a fully 3D unsaturated flow model based on (10) can be extremely computationally demanding. As a result, previous work in this domain has typically assumed that flow occurs only in the vertical direction (e.g., [START_REF] Binley | Vadose zone flow model uncertainty as conditioned on geophysical data[END_REF][START_REF] Cassiani | Modeling unsaturated flow in a layered formation under quasi-steady state conditions using geophysical data constraints[END_REF][START_REF] Looms | Identifying unsaturated hydraulic parameters using an integrated data fusion approach on cross-borehole geophysical data[END_REF][START_REF] Scholer | Bayesian Markov-Chain-Monte-Carlo inversion of time-lapse crosshole GPR data to characterize the vadose zone at the Arrenaes site, Denmark[END_REF], such that the following 1D version of Richards' equation can be utilized in the inversion procedure:

∂θ ∂t = ∂ ∂z K(h) ∂h ∂z + 1 . ( 13
)
The vertical flow assumption may hold in layered subsurface environments under natural loading conditions (e.g., [START_REF] Binley | Vadose zone flow model uncertainty as conditioned on geophysical data[END_REF], but it will be clearly violated during infiltration experiments where the area over which loading occurs is spatially restricted and loading rates are significantly higher. This will be particularly the case where there exist large contrasts in subsurface hydraulic properties (e.g., [START_REF] Looms | Identifying unsaturated hydraulic parameters using an integrated data fusion approach on cross-borehole geophysical data[END_REF][START_REF] Rucker | Inverse upscaling of hydraulic parameters during constant flux infiltration using borehole radar[END_REF]. As a result, the 1D flow assumption represents a significant source of model error. Although such errors and their potential for posterior parameter bias have been acknowledged in previous research efforts [START_REF] Scholer | Examining the information content of time-lapse crosshole GPR data collected under different infiltration conditions to estimate unsaturated soil hydraulic properties[END_REF], they have never before been examined and accounted for in the inversion procedure.

Infiltration experiment and data

Figure 2 shows the overall setup considered for our synthetic infiltration experiment. Infiltration at a rate of 2 cm/h is applied to a circular region on the Earth's surface having a diameter of 3 m. The infiltration is carried out for a period of 11.6 d, during which GPRderived estimates of horizontally averaged water content are considered to be available every 2.8 h. The water-content measurements are considered between boreholes 2-m apart and 8-m deep, with a depth sampling interval of 0.1 m. The subsurface consists of two layers whose VGM parameters are given in Table 1, which texturally describes a sandy soil underlain by a less permeable silt loam. The boundary between the layers is located at 3-m depth.

[Figure 2 about here.]

[Table 1 about here.]

To determine the spatiotemporal distribution of water content corresponding to the experimental setup described above, we used the code VS2D [START_REF] Lappala | Documentation of computer program VS2D to solve the equations of fluid flow in variably saturated porous media[END_REF] Figure 3 shows snapshots of the modeled subsurface water-content distribution over the course of the infiltration experiment for times t = 0, 1, 2, 5, 7, and 11 d. We see that, once the experiment begins, the infiltration front moves approximately vertically through the sandy soil layer, making its way to the boundary with the silt loam in just under 2 d. From this point onwards, although the front continues to move downwards, a strong lateral component to the flow is observed because of the lesser permeability of the lower layer. Indeed, as the water is not able to infiltrate as quickly into the silt loam, it begins to build up at the interface between the two soils and spread horizontally. Such behaviour cannot be captured using a 1D flow model based on equation ( 13), which is discussed in further detail below.

[Figure 3 about here.]

We next simulated the GPR-derived water-content measurements acquired during the infiltration experiment, which again represent the data to be inverted for the VGM parameters in each soil layer. To this end, every 2.8 h, the horizontal average of the water-content field between the boreholes was calculated from the VS2D results using a depth discretization interval of 0.1 m. This yielded 81 measurements in depth across 101 GPR acquisition times, to which zero-mean Gaussian random noise with a standard deviation of 0.01 (roughly 5%)

was added to simulate the effects of measurement error. It is important to note that, for the sake of simplicity in this example, we did not explicitly model the propagation of GPR energy between the transmitter and receiver antennas in the two boreholes based on the VS2D results, but rather assumed that the ZOP GPR experiment provided a measure of the horizontal average of soil water content as a function of depth. Although, in doing this, we admittedly neglect several aspects of the physics that would be encountered in a field setting such as critical refractions of GPR energy and frequency-dependent resolution limitations (e.g., [START_REF] Rossi | A stochastic analysis of cross-hole groundpenetrating radar zero-offset profiles for subsurface characterization[END_REF][START_REF] Rucker | Parameter estimation for soil hydraulic properties using zero-offset borehole radar[END_REF]), these aspects were not considered essential for this numerical study into the effects of model error arising from the 1D flow assumption.

Figure 4a shows the simulated GPR-derived water-content data, organized into a matrix with depth on the vertical axis and measurement time on the horizontal axis. To gain insight into the importance of model error for this example, Figure 4b shows the corresponding water content calculated as a function of depth and time assuming purely vertical flow, such that the 1D Richards' equation ( 13) could be applied. The results obtained using the 3D and 1D models for the same set of VGM parameters and boundary conditions are clearly and significantly different, most notably with respect to: (i) the speed at which the infiltration front travels through the lower layer, which is greater using the 1D model; and (ii) the evolution of water content in the upper layer after the infiltration front reaches the soil interface, in that the upper layer is seen to "fill up" to full saturation in the 1D case rather than pool and spread laterally. Figure 4c shows the difference between and 1D and 3D simulation results, equal to the sum of the model error and Gaussian measurement uncertainties. Here we see that there are parts of the data space where the magnitude of the error is almost 50%, and that the model error exhibits a high degree of correlation. All of this means that using a 1D flow model to stochastically invert the data in Figure 4a, without accounting for model error, will result in a strong bias in the estimated VGM parameters and unreliable posterior statistics (see Section 3.4). Finally, Figure 4d shows the error image from Figure 4c with the results reorganized such that they are plotted relative to the arrival time of the infiltration front as a function of depth observed in the data (Figure 4a). The importance of this data arrangement step is explained in the following section.

[Figure 4 about here.]

Model-error realizations and analysis

The first step in our approach to dealing with a known source model error in this paper involves generation of a set of stochastic model-error realizations corresponding to parameter sets randomly drawn from the Bayesian prior distribution. Again, this is done so that we can learn about the overall characteristics of the model error, with the goal of using this information to identify the model-error component of the residual during MCMC. Table 2 shows the lower and upper bounds of the uniform prior distributions that were assumed for the different VGM parameters in our synthetic study. Note that these distributions are rather broad and encompass a wide range of soil types (e.g., [START_REF] Carsel | Developing joint probability distributions of soil water retention characteristics[END_REF], and that the same priors were assumed for each soil layer. In this way, relatively little information about the hydraulic properties is provided to the inversion procedure and we rely strongly upon the data to resolve them. Also note that the prior bounds for K s are specified in terms of its logarithm, which is consistent with previous work and reflects the wide range of natural variability of this parameter (e.g., [START_REF] Scholer | Bayesian Markov-Chain-Monte-Carlo inversion of time-lapse crosshole GPR data to characterize the vadose zone at the Arrenaes site, Denmark[END_REF][START_REF] Scholer | Examining the information content of time-lapse crosshole GPR data collected under different infiltration conditions to estimate unsaturated soil hydraulic properties[END_REF].

[Table 2 about here.]

Each model-error realization was generated by: (i) drawing a random set of VGM parameters for each soil layer from the prior distributions in Table 2; (ii) computing the corresponding GPR-derived water-content data as a function of depth and time based on the general 3D

Richards' equation ( 10); (iii) computing the GPR-derived water-content data under the assumption of purely vertical flow using the 1D Richards' equation ( 13); and (iv) calculating the difference between the 3D and 1D simulation results. It is important to reiterate that this part of our model-error approach is easily parallelized in the sense that different model-error realizations can be computed on different processors of a cluster, thereby greatly reducing the time needed to run the relatively large number of expensive 3D unsaturated flow simulations required. In this regard, runs of the 3D solver for our example took approximately 100 s on a standard desktop computer, whereas runs of the 1D solver were over 60 times faster at 1.5 s.

Figure 5 shows an example of 18 model-error realizations, each of which has been plotted relative to the arrival time of the infiltration front observed in the 3D simulation results, as was done for Figure 4d. This latter step, whereby the realizations are effectively "aligned" on the curve representing the 3D infiltration-front arrival in depth, is important for this problem because, without it, the realizations would be highly dissimilar in the data space and not amenable to any kind of global analysis. In contrast, after alignment, the model-error realizations are seen to take on a similar form which is described by: (i) a triangular region below 3-m depth that results from the difference in the speed of propagation of the infiltration front in the lower layer between the 3D and 1D simulations; and (ii) another triangular region above 3-m depth that results when the upper layer "fills up" in the 1D simulation for cases where the lower layer is less permeable. Although the widths and amplitudes of these triangular regions are significantly different across the various realizations in Figure 5, the images, due to their similarity in form, are generally well suited to PCA analysis with the aim of generating a compact orthonormal basis for the model error. At the same time, however, it is important to note that the strong variations between the realizations in Figure 5 in terms of width and amplitude mean that the model error is not well described using a simple parametric distribution, and thus not amenable to the global statistical approaches for model error mentioned earlier. Indeed, detailed analysis of the 6500 model-error realizations indicates that the model-error values are highly non-Gaussian-distributed with complex correlation patterns in the data space.

[Figure 5 about here.]

To construct the model-error basis, a total of 6500 realizations were analyzed using PCA, the results of which showed that only the first 50 principal components (out of 6561) were necessary to capture 98% of the variance of the input. Note that the number of principal components required to capture this percent of variance tends to increase with the number of model-error realizations considered, as smaller sets of realizations will generally exhibit a lesser range of variability that can be represented by a smaller basis (Figure 6). Our choice of 6500 realizations represents a point after which this trend stabilizes and the addition of further realizations does not require more principal components to capture 98% of the variance.

Figure 7 shows the first 15, and last 3, vectors in the orthonormal model-error basis, ordered with respect to their decreasing contribution to the total variance and plotted as images in the data space. As expected, we see a gradual increase in the spatial frequency content of each vector as its index increases, with the first few vectors tending to capture the overall large-scale trends seen in the realizations in Figure 5 and the higher-order basis elements being necessary to resolve the finer details. Again, under the assumptions of orthogonality stated in Section 2.2, projection of the residual onto this basis during MCMC should adequately identify the model-error component, which can then be removed prior to computing the likelihood.

[Figure 6 about here.]

[Figure 7 about here.]

Stochastic inversion results

We now present the results of three different Bayesian-MCMC inversions to estimate the "true" VGM parameters in Table 1 from the GPR-derived water-content data in Figure 4a, all of which are based on use of a simplified 1D flow model. We begin by presenting the results of inverting using a "standard" Gaussian likelihood given by equation ( 3), where no correction for model error is considered and the standard deviation of the data errors s d is artificially inflated in order to compensate for the additional error source. This is followed by inverting using the informal 2 -norm-based likelihood measure given by equation ( 9), again with no correction for model error, such that the results obtained using this measure and using the Gaussian likelihood can be directly compared. Finally, we show the posterior results obtained for the case where the 2 -norm-based likelihood is combined with the correction for model error described in Section 2.2. For each inversion, a uniform MCMC proposal distribution Q(m |m i ), centered on the current state of the Markov chain and whose width was chosen to provide a model acceptance rate of approximately 30% [START_REF] Gelman | Efficient Metropolis jumping rules[END_REF], was employed.

A total of 800,000 MCMC iterations were run in each case, from which the first 10,000 samples were discarded as burn-in. These latter values were deemed appropriate based on visual inspection of each model parameter, along with its mean and variance, as a function of iteration (e.g., [START_REF] Hassan | Using Markov Chain Monte Carlo to quantify parameter uncertainty and its effect on predictions of a groundwater flow model[END_REF].

Gaussian likelihood, no model-error correction

Figure 8 shows the marginal posterior histograms obtained for the VGM parameters in each soil layer for the case where the data in Figure 4a were inverted using a standard Gaussian likelihood function. The error standard deviation in equation ( 3) in this case was arbitrarily set to s d = 0.2, which is 20 times the level of the random noise added to the data, in order to compensate for the model-error contribution to the residual and counteract the strong ranking of models provided by a Gaussian likelihood when the true residual statistics do not agree precisely with those that are assumed. Without such error inflation, the use of equation ( 3) would result in a highly peaked posterior distribution that could only be sampled with an extremely narrow proposal distribution and unreasonably large number of MCMC iterations. Indeed, [START_REF] Brynjarsdóttir | Learning about physical parameters: The importance of model discrepancy[END_REF] point out that, when model errors are present and not accounted for in Bayesian inference, the posterior tends to become narrowly focused around the wrong set of model parameters, and this only gets worse as more data are considered. Error inflation permits, at the very least, for the biased parameter set(s) to be identified at the expense of the posterior parameter uncertainties being arbitrary.

[Figure 8 about here.]

We see in Figure 8 that, because model error is present but has not been accounted for in the inversion procedure, the posterior VGM-parameter histograms are consistently focused on the wrong values. That is, there exists a set of incorrect parameter values whose predicted data, obtained using a 1D flow model, are a better match to the observed data than the true parameters in Table 1. The most significant bias in parameters occurs for the saturated hydraulic conductivity in both layers and the saturated water content of the upper layer, where the true values are seen to fall outside of the limits of the posterior distributions despite that fact that the error inflation imposed in this example is significant. As infiltration occurs significantly more rapidly in a 1D simulation than in 3D for the same set of model parameters (Figure 4), an inversion based on the 1D flow model will tend to select lower values for K s in both layers in order to best match the observed, 3D-generated data. Further, the 1D flow model predicts a greater accumulation of water at the interface between the two layers, which can be reduced by selecting lower values for θ s in the upper layer.

To gain insight into how such model-error-related biases translate into quantities relevant to flow and transport, Figure 9 shows the water retention and unsaturated hydraulic conductivity functions for the two soil layers corresponding to (i) the posterior VGM-parameter sets (color); (ii) the prior parameter ranges (grayscale); and (iii) the true parameter set in Table 1 (blue curve). Here we observe that the true curves often fall either at the limits of the posterior ranges or outside of them, meaning that the posterior parameter sets do not well reflect the soil hydraulic behaviour. Clearly, the model errors arising from the 1D vertical flow assumption cannot be neglected if we wish to have reliable predictions of flow and transport through this system.

[Figure 9 about here.]

3.4.2 L2-norm likelihood, no model-error correction

Figure 10 shows the marginal posterior histograms obtained for the case where the data in Figure 4a were inverted using the informal L 2 -norm-based likelihood measure developed in Section 2.3. Again, the advantage of using this measure is that the likelihood is determined based on the expected behaviour of a summary measure of the residual (i.e., its L 2 -norm), rather than on the residual vector itself, thereby avoiding an overly strong preference for model parameter sets whose corresponding residual statistics fit exactly the assumed Gaussian model.

As in the Gaussian likelihood case, no attempt was made to remove the effects of model error in this inversion. To account for the increased residual energy due to model error and allow for effective MCMC sampling, a modest amount of error inflation was made by setting the residual standard deviation in equation ( 9) to lie between s R 1 = 0.01 and s R 2 = 0.04.

[Figure 10 about here.]

We observe in Figure 10 that, as was the case with the standard Gaussian likelihood, a strong bias exists in the posterior results because of the model error coming from the 1D flow assumption. Indeed, the marginal posterior histograms look similar to those in Figure 8, with the true parameter values for K s and θ s often falling far outside of the limits of the posterior distributions. In terms of the water retention and unsaturated hydraulic conductivity functions, Figure 11 shows results that are almost identical to those in Figure 9. Note, however, that because of the use of the informal likelihood measure, the results presented here were obtained with significantly less error inflation than in the Gaussian likelihood case. That is, in using the L 2 -norm-based likelihood, we greatly increase the probability of acceptance of model parameter sets whose residual norm fits our expectations, but whose residual vector may deviate slightly from Gaussian.

[Figure 11 about here.]

L2-norm likelihood, correction for model error

Finally, Figure 12 shows the marginal posterior histograms obtained for the case where the water-content data in Figure 4a were inverted using our informal L 2 -norm-based likelihood measure combined with the proposed correction for model error described in Section 2.2.

Only a small amount of error inflation was done in this inversion by setting s R 1 = 0.01 and s R 2 = 0.015 in order to account for the fact that, even with the correction, it is unlikely that the model-error component of the residual will be perfectly removed. As a result, the energy in the remainder should be slightly larger than the level of noise added to the data.

[Figure 12 about here.]

We see in Figure 12 that, as a result identifying and subtracting the model-error component of the residual before evaluation of the likelihood in MCMC, the posterior VGM-parameter histograms are no longer biased, with the true parameter values falling in most cases near the middle of the posterior ranges. With regard to the corresponding hydraulic behaviour, we observe in Figure 13 that the true water retention and unsaturated hydraulic conductivity functions now lie well within the extent of the posterior curves. It is important to point out that, despite the fact that a minimal amount of error inflation was done for this inversion compared to the Gaussian-and informal-likelihood inversions, the posterior distributions are broader, most notably for parameters K s and θ s . This results from the fact that (i) a significant amount of residual energy is removed with our PCA-based correction before calculating the likelihood; and (ii) there exist some incorrect model-parameter sets whose corresponding parameter-error component of the residual will resemble (and will thus be identified as) model error, leading to the parameter sets being accepted in the MCMC inversion procedure. This latter important point is discussed in further detail in the following section.

[Figure 13 about here.]

Discussion

It is clear from the previous results that, in the context of the considered example problem, our proposed correction for model error offers an effective means of overcoming the posterior parameter bias related to use of a simplified forward model, thereby providing more accurate and useful uncertainty estimates. We now attempt to gain further insight into the reason why, with this correction, particular sets of incorrect model parameters may be accepted in the MCMC procedure, which contributes to the broadening of the obtained posterior distributions. Figure 14 presents the results of an analysis of three different parameter sets, the first row corresponding to the true subsurface VGM parameters (Table 1) and the last two rows corresponding to random "test" sets of VGM parameters drawn from the prior distribution (Table 3). In the columns of the figure we show (i) the predicted GPR-derived water-content data assuming 1D vertical flow; (ii) the residual obtained by subtracting the "observed" data in Figure 4a and expressing the results relative to the arrival time of the infiltration front;

(iii) the projection of this residual onto the model-error basis, which represents our estimate of the model-error component of the residual; and (iv) the corresponding remainder, obtained by subtracting the projection from the residual.

[Figure 14 about here.]

[Table 3 about here.]

We see in Figure 14 that, when the true set of VGM parameters is considered and thus when the parameter-error component of the residual is zero (see equation ( 6 The fact that certain sets of incorrect model parameters, whose residuals under the 1D

flow assumption appear similar to the stochastic model-error realizations, are given a high likelihood in our modified inversion procedure may be initially disconcerting. However, it must be emphasized that, in any situation where parameter-related errors cannot be distinguished from model errors, the corresponding model-parameter set cannot be rejected as a possibility.

Indeed, in this regard, our proposed algorithm should be viewed as a conservative stochastic inversion approach in the presence of model error; if a residual appears to resemble model error based on the generated model-error realizations, then the corresponding parameter set should not be excluded from the Bayesian posterior distributions. The strong advantage of our approach compared to not accounting for model error is that bias is strongly reduced and the true parameter set becomes well represented by these distributions.

Conclusions

Building on the recent work of [START_REF] Köpke | Accounting for model error in Bayesian solutions to hydrogeophysical inverse problems using a local basis approach[END_REF], we have presented in this paper a methodology for accounting for model errors in Bayesian-MCMC inversions that is geared towards the common case where such errors arise from the use of an computationally efficient simplified forward model in place of a more accurate but computationally burdensome numerical solution. Our approach is based on the analysis of a suite stochastic model-error realizations, created before the MCMC iterations by running the simplified and full forward solvers together for randomly drawn model-parameter sets from the prior distribution, which leads to the development of an orthonormal basis for the model error. Under the assumption that the model errors for the considered problem can be well described by this basis and that the model-error component of the residual lies orthogonal to the parameter-error and datameasurement-error components, projection of the residual onto the basis identifies the model error, which is then subtracted from the residual before evaluating the likelihood.

We saw through the considered example problem that application of our model-error cor-rection, combined with an informal likelihood measure based on the expected behaviour of the L 2 -norm, leads to a strong reduction in bias and notably better characterization of posterior uncertainties. This comes at the cost of needing to perform a number of full forward model runs (in our case a few thousand) to generate the model-error basis prior to MCMC.

Note again, however, that these full numerical simulations can be conducted in parallel. Our approach represents a remarkable computational savings when compared to MCMC based entirely on the full forward solver, in which hundreds of thousands of expensive model evaluations, conducted in series, would be necessary. For the specific example presented in this paper involving 800,000 Metropolis iterations, use of the fully 3D Richards' equation solver within MCMC would require over 900 days on a standard desktop computer. In contrast, running our algorithm based on the 1D forward solver with model-error correction took less that two weeks.

A critical assumption in our proposed approach is that of orthogonality between the model-and parameter-error components of the residual. As much as our experience until now suggests that this will be approximately true in many cases, hence explaining the success of our method, it cannot be proven and we observed that some incorrect model-parameter sets may produce a residual that looks like model error. In the latter cases, the model-and parameter-error components of the residual cannot be distinguished by projecting onto the basis and the parameter sets will stand a good chance of being accepted. In our view, this is not a concern as it simply means that the posterior parameter distributions will be broadened to include such parameters; i.e., our approach will conservatively include the parameters as possibilities. However, if this behaviour is undesirable, a two-stage MCMC algorithm could be proposed in which our approach would be used in a first accept/reject phase to effectively filter out unreasonable parameter sets from being tested with the full numerical solution, albeit at greatly increased computational cost.

It must be emphasized that the approach described herein is only intended for known sources of model error, for which random realizations of the error can be generated and used to help identify the model-error component of the residual. Although this will often be the major source of bias for inverse problems in geophysics and hydrology, there are situations where even our best forward solution will not provide a good enough description of the physical process involved. In the case of such unknown or unspecified model errors, our methodology can still be expected to effectively deal with the model errors for which it was intended, thereby providing more reliable posterior uncertainty estimates. Another related issue is the fact that, for many problems, particularly those of high dimension with spatially distributed parameters, the nature of the model errors may change significantly over the model parameter space and it may not be possible to effectively describe them using a single global basis. In this case, the work of [START_REF] Köpke | Accounting for model error in Bayesian solutions to hydrogeophysical inverse problems using a local basis approach[END_REF] shows that a KNN dictionary-based approach to modelerror identification, whereby the basis is constructed locally at each MCMC iteration, can be a highly effective means of obtaining reliable posterior parameter distributions when using an approximate forward solver. It is also likely that the computational efficiency of the approach 2) and the posterior distribution obtained using an inflated Gaussian likelihood function for the residuals with no model-error correction (color; Figure 8). The blue lines represent the curves corresponding to the true parameter set in Table 1. The prior and posterior results are expressed in terms of curve densities. 2) and the posterior distribution obtained using an L 2 -norm-based likelihood measure for the residuals with no model-error correction (color; Figure 10). The blue lines represent the curves corresponding to the true parameter set in Table 1. The prior and posterior results are expressed in terms of curve densities. Figure 14: For the true set of VGM parameters (Table 1) and two sets of incorrect "test" model parameters (Table 3)

  4. Perform principal component analysis (PCA) on the model-error realizations {E 1 , ..., E k } in order to obtain a sparse orthonormal basis B = [b 1 , ..., b b ] for the model error. The number of basis vectors b should be chosen to be the minimum required to capture a high percentage of the variance of the realizations, typically around 98%. In this way, the basis will be able to capture the model-error behaviour, but will have minimal ability to represent contributions to the residual that do not resemble model error such as data measurement uncertainties. 5. For each set of model parameters m tested within MCMC, calculate the best approximation of the residual r(m ) = F (m ) -d obs using the model-error basis, obtained in a least-squares sense using B B T r(m ), and remove this result from the residual. This yields the remainder

  to solve the general 3DRichards' equation (10) under the assumption of rotational symmetry about the vertical axis, meaning that the model domain was parameterized in terms of radius (r) and depth (z), with r = 0 corresponding to the center of the infiltration region. A specified-flux boundary condition was imposed at the Earth's surface (z = 0 m) with no-flow conditions assumed outside of the infiltration region (r > 1.5 m). No-flow conditions were also assumed along the outside of the model domain, the latter of which was set at r = 4 m. At the bottom of the domain (z = 10 m), a fixed-pressure-head value of h = -0.5 m was specified in order to simulate the presence of the water table at 10.5-m depth. The initial distribution of soil water content prior to running the infiltration experiment was obtained using a 1D steady-state infiltration code based on the work of[START_REF] Rockhold | An analytical solution technique for one-dimensional, steady vertical water flow in layered soils[END_REF] assuming a constant infiltration rate of 0.036 cm/h.

  )), projection of the residual onto the model-error basis correctly identifies the model-error component, which after subtraction leaves a low-amplitude remainder that is mostly comprised of Gaussian datameasurement uncertainties. As the L 2 -norm of the remainder determines the likelihood, the true parameter set stands a high chance of being accepted in MCMC. For the first (incorrect) set of test model parameters, we observe that the corresponding residual, which now is comprised of non-zero model-and parameter-error components, closely resembles the model-error realizations presented in Figure5. Ideally, projection of this residual onto the PCA-derived basis would isolate only the model-error component. However, in this case the entire residual is identified as model error, which again results in a low-amplitude remainder and a correspondingly high probability of acceptance. In other words, when the sum of the model-and parameter-error components of the residual tends to look similar to the stochastic model-error realizations, both of these components will be subtracted in our correction procedure, leading to a high likelihood of an incorrect parameter set. Finally, for the last set of test model parameters, we observe the intended functioning of the algorithm; the projection of the resid-ual onto the model-error basis correctly identifies the model-error component, but leaves the parameter-error component which then forms part of the remainder. The high amplitudes observed in the remainder yield a low probability of the parameter set being accepted.

ofFigure 1 :

 1 Figure 1: Normalized 2 -norm likelihood given by equation (9) as a function of ||R(m)|| for N = 1000 with s R 1 = 0.1 and s R 2 = 0.2.

Figure 2 :

 2 Figure 2: Setup for the synthetic infiltration experiment considered in this study. The white dots represent transmitter and receiver antenna positions for the ZOP GPR measurements.

Figure 3 :

 3 Figure 3: Spatial distribution of water content in the subsurface at various times throughout the infiltration experiment. The GPR boreholes are shown for reference.

Figure 4 :

 4 Figure 4: (a) Average soil water content between the boreholes as a function of depth and measurement time, computed using a 3D infiltration model with the addition of Gaussian measurement noise, representing the synthetic data to be inverted for the VGM parameters in Table 1. The arrival time of the infiltration front as a function of depth is indicated with a black dashed line. (b) Corresponding water-content distribution obtained assuming purely vertical (1D) flow. (c) Difference (b)-(a), which is equal to the sum of the model error and measurement uncertainties. (d) Error image from (c) expressed relative to the arrival time of the infiltration front in (a).

Figure 5 :

 5 Figure 5: Example stochastic realizations of the model error corresponding to random sets of VGM parameters drawn from the prior distributions in Table2. For greater coherency between the images, each has been expressed relative to the arrival time of the infiltration front in depth as observed in the 3D flow simulation. A total of 6500 realizations were generated to construct the model-error basis.

Figure 8 :

 8 Figure 8: Marginal posterior histograms for the VGM parameters in each soil layer, obtained through MCMC sampling using an inflated Gaussian likelihood function for the residuals with no model-error correction. The red dots indicate the true parameter values. The limits of the horizontal axis on each plot represent the prior uniform parameter bounds.

Figure 9 :

 9 Figure 9: Water retention (left) and unsaturated hydraulic conductivity (right) functions for each soil layer corresponding to the prior distribution (gray; Table2) and the posterior distribution obtained using an inflated Gaussian likelihood function for the residuals with no model-error correction (color; Figure8). The blue lines represent the curves corresponding to the true parameter set in Table1. The prior and posterior results are expressed in terms of curve densities.

Figure 10 :

 10 Figure 10: Marginal posterior histograms for the VGM parameters in each soil layer, obtained through MCMC sampling using an L 2 -norm-based likelihood measure for the residuals with no model-error correction (see text for details). The red dots indicate the true parameter values. The limits of the horizontal axis on each plot represent the prior uniform parameter bounds.

Figure 11 :

 11 Figure 11: Water retention (left) and unsaturated hydraulic conductivity (right) functions for each soil layer corresponding to the prior distribution (gray; Table2) and the posterior distribution obtained using an L 2 -norm-based likelihood measure for the residuals with no model-error correction (color; Figure10). The blue lines represent the curves corresponding to the true parameter set in Table1. The prior and posterior results are expressed in terms of curve densities.

Figure 12 :

 12 Figure 12: Marginal posterior histograms for the VGM parameters in each soil layer, obtained through MCMC sampling using an L 2 -norm-based likelihood measure for the residuals after correcting for model error (see text for details). The red dots indicate the true parameter values. The limits of the horizontal axis on each plot represent the prior uniform parameter bounds.

  Figure14: For the true set of VGM parameters (Table1) and two sets of incorrect "test" model parameters (Table3): (a) Predicted water content assuming 1D vertical flow; (b) Residual obtained by subtracting the synthetic data in Figure 4a from the results in (a) and expressing relative to the arrival time of the infiltration front observed in the data; (c) Projection of the residual in (b) onto the model-error basis; (d) Corresponding remainder (b)-(c).

Table 1 :

 1 log 10 (K s [m/s]) θ r [-] α [1/m] n [-] θ s [-] VGM parameters considered for the 2-layer synthetic example. 48 log 10 (K s [m/s]) θ r [-] α [1/m] n [-] θ s [-]

	Layer 1	-4.074	0.036 12.552 2.830 0.405
	Layer 2	-5.826	0.077	3.165	1.819 0.462

Table 2 :

 2 Lower and upper bounds of the uniform Bayesian prior distributions assumed for the VGM parameters in each layer.

Acknowledgments

This work was supported by a grant to J. Irving from the Swiss National Science Foundation, grant number 200021-140864. The authors also wish to thank Richard Healey from the U.S. Geological Survey for support with the VS2D software during the initial stages of this research project.

the assumption of model error, Statistical applications in genetics and molecular biology, 12 (2), 129-141. Xu, T., and A. J. Valocchi (2015), A Bayesian approach to improved calibration and prediction of groundwater models with structural error, Water Resources Research, 51 (11), 9290-9311. 2) and the posterior distribution obtained using an L 2 -norm-based likelihood measure for the residuals after correcting for model error (color; Figure 12). The blue lines represent the curves corresponding to the true parameter set in Table 1. The prior and posterior results are expressed in terms of curve densities.

Number of model-error realizations


Lucie Cez

Régis Vallet

Lionel Darras

The Hydraulic System of Larsa (South Iraq)

Since the reopening of the mission in 2019, the water supply of Larsa has been a central focus of field investigations on this urban site of South Mesopotamia (Iraq). The first map of its hydraulic network can be established thanks to the cross-analysis of photogrammetric and geomagnetic images, to which is added the data of the geo-archaeological study and archaeological operations. In its final state, the intra-mural network had more than a hundred canals including several large ones such as the Grand canal, more than 22 m wide, basins of storage, the largest reserve in the north-west had a 26,000 m 2 , irrigation and drainage canals, but also a port area of 31,645 m 2 to the south. This approach made it possible to establish the connection with the canals works of the surrounding plain structuring the peri-urban network and to apprehend the major developments during the 2 nd millennium, between the

Introduction

Although it was crucial to the development of Larsa (Tell es Senkereh) (Fig. 1), the 200 ha city dominating the Mesopotamia at the beginning of the 2 nd millennium BC, its water supply has never been specifically studied. The existence of large canals crossing the Uruk plain, some of them located in the vicinity of the site and oriented towards it, is known thanks to the archaeological and geomorphological surveys as well as the historical sources [START_REF] Adams Mc | The Uruk Countryside: The Natural Setting of Urban Societies[END_REF][START_REF] Mc | Heartlands of Cities Surveys of Ancient Settlement and Land Use on the Central Floodplain of the Euphrates[END_REF][START_REF] Geyer | Nouvelle contribution à l'étude géomorphologique de la région de Larsa-'Oueili (Iraq)[END_REF][START_REF] Pournelle | Marshland of Cities: Deltaic Landscapes and the Evolution of Early Mesopotamian Civilization[END_REF][START_REF] Steinkeller | New Light on the Hydrology and Topography of Southern Babylonia in the Third Millenium[END_REF]. Until the resumption of fieldwork in 2019 [START_REF] Vallet | Larsa-'Uwaili Annual Report[END_REF], Vallet et al. 2020[START_REF] Vallet | Larsa-'Uwaili Annual Report 2021-2022[END_REF][START_REF] Vallet | Larsa and Tell el-'Uwaili (Iraq), Preliminary Results[END_REF], the available data was limited to large traces visible on the series of aerial photographs from G. Gerster (1973) to the north of the site. The imprint of large structures will give rise to divergent interpretations: as large roadways or large canals according to the authors [START_REF] Huot | La structure urbaine de Larsa, une approche provisoire[END_REF]Margueron 2013: fig. 111), none however relayed by archaeological data.

The study of the hydraulic system of Larsa consisted first of all in reconstructing the organization of the hydraulic works inside the city, that is the network of canals allowing its supply and including possible basins allowing the water management. At a larger scale, it is also a question of establishing the connexion between internal hydraulic network and large canals crossing the lowland of the Euphrates River where Larsa was located. Finally, the objective pursued was to understand the general functioning of the system (water source, flow direction), specifically due to its insertion within the urban fabric, as well as to establish the history of its evolution over time.

A cross-analysis for identifying the canals network

To meet these objectives, an interdisciplinary and multi-scalar research was implemented in 2019. This is based on the contributions of geophysics, photogrammetry and geo-archaeological field study and archaeological operations [START_REF] Vallet | Larsa-'Uwaili Annual Report[END_REF]Vallet, Cez and Darras 2022: 34-36).

In order to highlight the layout of ancient canal works and other hydraulic structures, the carto-photointerpretation work relied on the ortho-photos by drone [START_REF] Suire | Note on Photogrammetric Surveying and Mapping of Larsa[END_REF], satellite image and ancient aerial photographs of the site, crossed with the geophysical image from the geomagnetic survey. 70 hectares of the site and its immediate surroundings were explored over the period 2019-2023 with the geomagnetic method, which provides an image of the built and dug structures of the basement, up to 2 m deep [START_REF] Aspinall | Magnetometry for Archaeologists[END_REF]. All of this high-resolution and georeferenced documentation has been integrated into a Geographic Information System (GIS) built using QGis software, in which the map of the hydraulic system has been drawn.

On the field, the geo-archaeological survey enabled to identify many canals, detectable by their typical morphology of large depression delimited by marked levees, and whose upper sandy filling is punctuated by scattered xerophyte shrubby vegetation (Fig. 2). Following their relatively rectilinear route on the ground over a certain distance makes it possible to distinguish constructed canals from wadi beds, especially if branches are still visible. The significant silting up of the site, mainly in the form of barchan-type dunes and sand toppings of varying thickness, however, makes annual photogrammetric coverage essential. With the dunes displacement, structures can be thus detected in previously silted sectors.

The map of the Larsa hydraulic system

This methodological approach allowed the identification of more than a hundred of canal layouts of variable sizes and eighteen storage basins on which the hydraulic system of Larsa, in its final state, was based (Fig. 3). Such a very complex network extended over approximately one 1,200 ha including the intra-mural network and the peri-urban one, no less dense. In total, 124 canals were numbered in the first map produced in 2022 (Vallet, Cez and Darras 2022: 35) but this is bound to evolve because some sectors are still unexplored in geomagnetic prospecting. Remember that less than half the size of the Larsa site, of 200 ha of extension, has been covered to date by geophysics survey. The main limit to the detection of fossil urban and hydraulic structures, however, remains significant silting up of the eastern part of the city, which have provided much less data, and that of the area south of the site where coalescent dune fields cover the plain.

The north-east sector

Originating from the northeast, the canal labelled C1 in the general map of the hydraulic system was the most important structure for water supply of the city. Its layout has been identified over a distance of 1.2 km thanks to the drone coverage but previous archaeological and palaeohydrological surveys highlighted from the Tigris near Zabalam, 62 km from Larsa, to Bad Tibira, 17 km away, and beyond.

But its final section, the connection to Larsa, was missing (Adams and Nissen 1972: 36-39;Pournelle 2003: fig. 76). This canal was established as the main supply of the city during the reign of Siniddinam (1849-1843 B.C.) (Sigrist 1990: 24;Frayne 1990: 158-160;[START_REF] Steinkeller | New Light on the Hydrology and Topography of Southern Babylonia in the Third Millenium[END_REF]. This latter work probably replaced an older canal of lesser importance, part of which was dated to the Ur III period (2112-2004 B.C.) (Adams and Nissen 1972: 39). Its entry point into the city was located about 100 m north of the Parrot Gate (East Gate B1), where a massive gate allowed passage through the rampart and here called watergate B65 (Fig. 3). The C1 canal work, 30 to 20 m wide in the plain, was reduced to 5 m wide at the crossing point in the city as it is shown by the geomagnetic image (Fig. 4).

Reducing the size made it possible to increase the flow speed and the water height. The powerful earth dykes framing the canal all along its course on the plain (D1-D2) were reinforced before it entered the city and associated with a lateral embankment to the south, itself topped by a baked brick wall. Such a device was intended to support the pressure exerted by the water in the canal work due to its curved layout. The Watergate B65 itself was fitted out for the water circulation, notably the semi-circular shape of the two internal faces of the baked-brick piers, 43 m in length for a maximum of 20 m large.

The C5 canal was designed as a relief canal for part of the flow of the C1 canal before its entry into the city, while the valve C6 had to evacuate the overflow of water towards the plain, the latter also having a defensive role since it allowed the surroundings of the eastern gate of the city to be flooded4 .

Once through the Watergate B65, water was transported to the center of the city by a succession of large canals forming a loop (C1-C4) (Fig. 3). The canal C2 formed a large curve to the north over 1.2 km, then continued south with the canal C3 until the ziggurat in the city centre and, finally, turns to the southeast by the narrower canal C4, to leave the city. The 'Grand Canal' (C3) is the main artery of the internal system, 22 to 25 m wide, and was flanked on the west by a large landing stage and a high dyke, revealed by geomagnetic prospecting (Fig. 5). The geophysical image clearly shows the piers of two bridges at its northern and southern ends, extended by causeways of which some portions of their rectilinear lines have been found (see roads 51, R2, R3 and R10 in Fig. 3). This shows that the hydraulic works were perfectly integrated into the urban structure of Larsa. The ongoing excavation of the east pier of the south bridge B69 has already revealed 8 courses of baked bricks for a depth of nearly 1 m (Fig. 6a) (Zingarello, Bachelot and Hsony 2022: fig. 2). The upper's courses of the pillar yielded an inscribed brick of king Sin-iddinam (Charpin and Jacquet 2022: fig. 6) (Fig. 6b).

The north-west sector

The canal C2 supplied the northwestern part of the city by the C7 canal of diversion reinforced with two side dykes due to the large quantities of water redistributed here (Fig. 5). The water of C7 was also diverted towards the west by other canals of variable sizes. The other work linking the Grand Canal and the western sector is the C9 canal, which runs along the Nur-Adad Palace and joins the P1 basin at its northwest extremity (Fig. 3).

The distinctive feature of the western zone from the point of view of its hydraulic works is that it houses a set of water storage basins. The P1 storage system was composed of several and connected basins covering a 6,640 m 2 area (Fig. 3). Its modular plan probably results from successive additions over time. The most important storage work of the city was the huge P2 reserve of 26,000 m 2 , whose clear signature on ortho-photos by drone (Fig. 7) is partly explained on the field by the density of carbonate concretions and mollusc. The location of such a structure takes into account the flat topography of this western sector of the city. These basins were indirectly fed by the Grand Canal via the C7 secondary canal for the basin P1 and via C9-C15 canals for the basin P2.

The geomagnetic prospecting also revealed an irrigated plot on the edge of buildings B48 and B49 and connected to this P2 reserve by canals C20 and C22 (Darras 2020: fig. 5; Vallet, Cez and Darras: fig.

14) (Fig. 3). A network of irrigation canals intersecting at right angles fed this geometric plot, while a drainage canal C23 ran from its southwest corner in order to evacuate waters outside of the city. Many incomplete sections of canals of this type around the P2 reserve suggest that this must have been an area of the city dedicated to irrigated agriculture5 .

The south sector

Spatial organization and functioning of the south hydraulic network is not yet fully understood. High silting up and numerous gullies of wadi blurred the ancient hydraulic works traces or has widened and scoured their original profile. The canal C26 is the major work in the southern zone, starting from the central basin P5 and flowing southeast to the Watergate B17 designed according to the same scheme as the Watergate B65 to the northeast, but with more modest dimensions of 37 m in length for 10 m large (Fig. 3). The central basin P5 probably collected rainwater, surrounding flows and the waters drainage of the E-babbar, but it should be also connected with the city network from the northeast. The other important work was the C27 canal, with a still incomplete layout to this day, but which has branches of secondary canals.

At the periphery of the southern area, a vast port zone of 30,000 m 2 was identified on either sides of the rampart (Fig. 8). The external basin followed the enclosure contour, protected by a 17 m wide dyke. Access to the internal smaller basin was by the west through the rampart, by two passages of 4 m wide and 25 m apart. Access to the port from the south was by a canal ten meters wide. The limits of the inner basin are less clear but had probably two docks, a main and a secondary. Associated with this port basin in the city, the function of the large rectangular building B57, of more of 4,000 m 2 , known for a long time and that yielded several inscribed bricks of king Sin-iddinam, finds a new interpretation (Huot, Rougeulle and Suire 1989: 50;Vallet and Suire 2020: 58). It was probably devoted to storage in the view of its partitioned internal organization, and more specifically a warehouse of the port.

The peri-urban hydraulic network

Around the city, the density of the peri-urban network is also impressive, especially since the barchan fields prevents the visibility of a large number of ancient canals in the south of the plain. Everywhere else in the close vicinity of the city, the deflation exerted on the former topography by the prevailing winds reveals their fossil layouts. Except the C1 canal to the northeast which continues into the city, the large peri-urban canals bring water to key points, where it is collected and redistributed. The many observable overlaps between the canal works show that not all of them were operational at the same time. An analysis of overlaps and doubling of canals is therefore to be done to periodize and characterize the different states of the network over time (Fig. 9). Four phases could be identified, corroborated by historical data, enhancing successively the different parts of the city during the 2 nd millennium B.C. (Fig. 3).

The oldest system (phase 1) established west of the city is attributable to king Gungunum (1932-1906 B.C.), sovereign which initially erected the city wall enclosure in the 21 st year of its reign [START_REF] Sigrist | Larsa Year Names[END_REF]: 9-10) and who carried out hydraulic works, which is reflected in the large number of 'year names' celebrating the digging of canals (Charpin 2002: tabl. 1). Two large canals organized the network, each 15 m wide, C37 and C88, fed by the Iturungal and the Euphrates (Fig. 3). The water was thus brought to the closest to the city at the Watergate C19. Inside the city, the canals C19 and C16 should correspond to this most ancient system. The extension of this first network to a larger territory (phase 2) was carried out under the reigns of Abi-sare (1905Abi-sare ( -1895 B.C.) B.C.), son and successor of Gungunum, then Sumu-El (1894-1866 B.C.). Historical sources attribute the strengthening of the city wall and the digging of the 'hiritum' water ditch to the king Abi-sare (Frayne 1990: 149;Sigrist 1990: 12-13). According to our hydraulic system reconstruction, the hydraulic works of this phase were the large basin P9 redistributing water to the west by the canal C39 and to the north by the canal C95 (Fig. 3). Several mapped hydraulic works could correspond to the 'hiritum, among them the internal canal C19 or the peri-urban ones C34, C36 or C446 .

Further developments of the hydraulic network (phase 3) took place during the reign of the king Nur-Adad (1865-1850 B.C.), especially in the northwest and southwest parts, which have been split. The large canal C38 was built about 200 m west and parallel to the former one C88, fed by the Iturungal or further by Euphrates River (Fig. 3). The same process can be observed with the canal C50 coming from the northwest and which splits the previous C49. Such developments made it possible to increase the water supply for the city in the political context of rivalry between Isin and Larsa cities and which had led to a break in the water supply of the latter, under the reign of Nur-Adad (Charpin 2002: 556-557;Charpin 2004: 102). The most recent system (phase 4) corresponds to the large hydraulic works program of the king Siniddinam (1849-1843 B.C.) previously described (infra). It specifically develops the eastern part of the network, including the C1-C4 loop, which involves a new water source from the Tigris River (Fig. 3). Such a reorganization of the hydraulic system is understandable in light of the difficulties of its predecessor Nur-Adad to insure the fresh water to its city from the Euphrates (Charpin and Jacquet 2022: 148). The canal C1 coming straight from Bad Tibira, about 17 km northeast of the site, is the backbone of the whole system: not only does the C1 canal supply the city's internal network but numerous canals divert its water to the north before entering Larsa. This northern network developed from several branches derived from C1 canal (C94, C95, C92) and integrated old sections of palaeochannels, more or less regularized (C90, C91 and a section west-east of the C92).

The entire system became considerably denser over time required large-scale water management. This was the objective of the huge reservoirs located in the plain, about 1.2 km north of the city wall. They was located in the heart of a vast agricultural area where other smaller reserves of water, canals and fossil plots have been identified during the last field campaigns of 2022. These structures had a geometric and modular plan resulting from successive enlargements of its outer limits to meet growing water needs.

Conclusion

Our knowledge of the general urban organization of Larsa and of its hydraulic network has made great strides since the resumption of fieldworks in 2019. The first map of its water supply system presented to the 13 th annual conference of ICAANE was produced with support of photogrammetry, geophysics and geoarchaeology, an effective methodological approach enabling to identify fossil traces of hydraulic structures and their functional interpretation. The coexistence and high density of hydraulic works for water conveyance (canals of large templates and secondary), storage (basins of variable sizes, large reserve), irrigation and drainage (canals and canalizations) and protection (dykes) gives the

Fig. 1 .

 1 Fig. 1. The Tell of Larsa with contour lines (1 m) in the lower plain of the Euphrates River (southern Iraq).© Mission of Larsa-Uwaili, 2022.

Fig. 2 .

 2 Fig. 2. The 'Grand Canal' C3 looking north as a depression filled with sand and shrubs. © L. Cez, November 2021.

Fig. 3 .

 3 Fig. 3. General plan of the urban structure of Larsa and its hydraulic system. CAD R. Vallet, L. Cez, A. Reiff (IDAAR). © Mission of Larsa-Uwaili, 2022.

Fig. 4 .

 4 Fig. 4. Geomagnetic image showing the canal C1 before entering the city, with lateral embankments and the Watergate B65. © L. Darras and Mission of Larsa-Uwaili, 2022.

Fig. 5 .

 5 Fig. 5. Geomagnetic image showing the hydraulic network of the north of Larsa. White box highlights the magnetic anomalies of the baked bricks piers of the bridge B25 (north). © L. Darras and Mission of Larsa-Uwaili, 2022.

Fig. 6 .

 6 Fig. 6. The east pillar of the south bridge B69 of the Grand canal (a). Inscribed brick of king Sin-iddiman (b). © Mission of Larsa-Uwaili, 2022.

Fig. 8 .

 8 Fig. 8. The port area of Larsa in the south of the city. Drone image by J. Lisein. © Mission Larsa-Uwaili, 2019.

Fig. 9 .

 9 Fig. 9. Superposition of three generation of canals of the peri-urban network west of Larsa.Drone image by J. Lisein. © Mission Larsa-Uwaili, 2019.

Pantheon-Sorbonne Paris 1 University, UMR 7041 ArScAn, Nanterre, France. Corresponding author : luciecez@gmail.com

CNRS, UMR 8068 TEMPS, Nanterre, France, regis.vallet@cnrs.fr

CNRS, UMR 5133 Archéorient, Lyon, France, lionel.darras@mom.fr

Fan-shaped flood traces originating at the outlet of the valve and the plugging of the Parrot Gate entrance with raw brick as well as abundant projectiles found on-site support this hypothesis(Parrot 1933: 177 ; Huot, Rougeulle and Suire 1989: 40). The technical characteristics of this hydraulic system and in particular of the north-eastern part of the city articulated around the Watergate B65 will be detailed in a special publication to come.

The nature of artefacts uncovered by the pedestrian survey[START_REF] Giraud | Rapport final d'activités scientifiques 2019-2022[END_REF], with sickle blades and lithic material, is consistent with an agricultural function dominant in the sector.

It is not established that this hydraulic structure is necessarily external to the rampart like a moat (seeCharpin and Jacquet 2022: 147-148).