\makeatletter
\@ifundefined{HCode}
{\documentclass[screen, CRGEOS, Unicode, Thematic]{cedram}
\usepackage[square,comma,authoryear]{natbib}
\renewcommand\bibsection{\section*{\refname}}
\newenvironment{noXML}{}{}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{}%\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.7\linewidth}{#1}}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\ndash{\text{--}}
\def\og{\guillemotleft}
\def\fg{\guillemotright}
\def\0{\phantom{0}}
\def\xoverrightarrow#1{\overrightarrow{#1}}
\makeatletter
\g@addto@macro{\UrlBreaks}{\UrlOrds}
\gappto{\UrlBreaks}{\UrlOrds}
\RequirePackage{etoolbox}
\usepackage{upgreek}
\def\jobid{crgeos20230417}
%\graphicspath{{/tmp/\jobid_figs/web/}}
\graphicspath{{./figures/}}
\newcounter{runlevel}
\usepackage{multirow}
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\let\MakeYrStrItalic\relax
\def\refinput#1{}
\let\ubreak\break
\csdef{Seqnsplit}{\\}
\def\botline{\\\hline}
\def\cbotline{\\\hline}
\def\back#1{}
\def\rmmu{\upmu}
\DOI{10.5802/crgeos.247}
\datereceived{2023-04-26}
\daterevised{2023-10-18}
\dateaccepted{2023-11-24}
\ItHasTeXPublished
}
{
\PassOptionsToPackage{authoryear}{natbib}
\documentclass[crgeos]{article}
\usepackage[T1]{fontenc}
\def\CDRdoi{10.5802/crgeos.247}
\def\xmorerows#1#2{\morerows{#1}{#2}}
\def\cbotline{\\\botline}
\let\ubreak\relax
\makeatletter
\def\CDRsupplementaryTwotypes#1#2{}
\def\xoverrightarrow#1{\vec{#1}}
}
\makeatother

\usepackage{upgreek}

\dateposted{2023-12-18}
\begin{document}

\begin{noXML}

%\makeatletter
%\def\TITREspecial{\relax}
%\def\cdr@specialtitle@english{Geo-hydrological Data \& Models}
%%\def\cdr@specialtitle@french{G\'eo-hydrologie, donn\'ees \& mod\`eles}
%\makeatother

\title{Cokriging of transmissivity from head measurements revisited in
the case of two quasi steady state flows}

\alttitle{Cokrigeage de la transmissivit\'{e} \`{a} partir de mesures
de la charge, revisit\'{e} dans le cas de deux \'{e}coulements en
r\'{e}gime quasi permanent}

\author{\firstname{Chantal} \lastname{de Fouquet}\CDRorcid{0000-0002-0008-4801}\IsCorresp}
\address{Mines Paris, PSL, Centre de g\'{e}osciences, 77300
Fontainebleau, France}
\email[C. de Fouquet]{chantal.de\_fouquet@minesparis.psl.eu}

\author{\firstname{J\'{e}r\'{e}my} \lastname{Nos}\CDRorcid{0000-0002-9315-7554}}
\address{Mines Paris, PSL, France}
\email[J. Nos]{jeremy.nos@mines-paris.org}

\keywords{\kwd{Inverse problem}\kwd{Flow equation}\kwd{Steady
state}\kwd{Cokriging}\kwd{Geostatistics}}

\altkeywords{\kwd{Probl\`{e}me inverse}\kwd{\'{E}quation
d'\'{e}coulement}\kwd{R\'{e}gime
permanent}\kwd{Cokrigeage}\kwd{G\'{e}ostatistique}}

\thanks{The presented study was funded by Schlumberger Water Service,
France.}

\begin{abstract}
The theoretical relationships between simple and cross-covariances of
the transmissivity and the hydraulic head, deduced from the
linearisation of the flow equation are used for cokriging
transmissivity field with data from one then two quasi steady-state
flows. First a non-exhaustive historical review of geostatistical
modelling based on flow equations in a simplified case is given. Then
original results show that the estimation of the logtransmissivity by
cokriging with head data from two flows is improved compared with the
case of head data from a single flow. The estimation in the presence of
a transmissivity barrier is thoroughly investigated. Practical
conclusions are drawn and perspectives based on the recent numerical
variogram are proposed.
\end{abstract}

\begin{altabstract}
Les relations th\'{e}oriques entre les covariances simples et
crois\'{e}es de la transmissivit\'{e} et de la charge hydraulique,
d\'{e}duites de la lin\'{e}arisation de l'\'{e}quation
d'\'{e}coulement, sont utilis\'{e}es pour le cokrigeage du champ de
transmissivit\'{e} avec des donn\'{e}es provenant d'un puis de deux
\'{e}coulements en r\'{e}gime quasi permanent. Une revue historique
non-exhaustive de la mod\'{e}lisation g\'{e}ostatistique bas\'{e}e sur
les \'{e}quations d'\'{e}coulement dans un cas simplifi\'{e} est
d'abord donn\'{e}e. Ensuite, des r\'{e}sultats originaux montrent que
l'estimation de la logtransmissivit\'{e} par cokrigeage avec des
donn\'{e}es de charge provenant de deux \'{e}coulements est
am\'{e}lior\'{e}e par rapport au cas des donn\'{e}es de charge
provenant d'un seul \'{e}coulement. L'estimation en pr\'{e}sence d'une
barri\`{e}re de transmissivit\'{e} est \'{e}tudi\'{e}e de mani\`{e}re
approfondie. Des conclusions pratiques sont tir\'{e}es et des
perspectives bas\'{e}es sur le r\'{e}cent variogramme num\'{e}rique
sont propos\'{e}es.
\end{altabstract}

\maketitle
\vspace*{.25pc}

\twocolumngrid

\end{noXML}

\section{Introduction} \label{sec1}

In \textit{four decades on inverse problem in hydrogeology}, \citet[see
also \citealp{deMarsilyetal1999}]{deMarsilyetal2000} observe that ``in
hydrogeology, there are a great number of competing methods for solving
the inverse problem''. \citet{Kitanidis1995} notes: ``although the
consensus is that the linear method is appropriate only when the
variance of log conductivity is less than~1, in my experience the
method may work reliably for larger variances provided that the
variation happens to be gradual or in the direction perpendicular to
the streamlines''. However, referring to the comparison of methods for
solving the inverse problem \mbox{carried} out by several teams on four
synthetic data sets,{\break} \citet{Zimmermanetal1998} point out ``the
evident superiority of the nonlinear methods as compared with linear
ones''. In a more recent review of the inverse problem,
\citet{Illman2014} presents the evolution of approaches for mapping
heterogeneity. 

In this paper, we come back to the early cokriging of the
transmissivity from transmissivity and head measurements, based on the
linearisation of the flow equation linking head and transmissivity
fluctuations. First a non-exhaustive historical review of
geostatistical modelling based on the linearisation of flow equations
is given. The relationships between simple and cross-covariances of
transmissivity and head fluctuations for two quasi-steady state flows
established by \citet{ZhangNeuman1996} in another context are applied
to the inverse problem.{\break} On synthetic cases original results
show the interest of crossing two quasi-steady-state flows to improve
the estimation of transmissivity. The restrictive hypotheses used for
the linearisation of the flow equation are then discussed, as well as
the practical application of cokriging of log-transmissivity. 

The present paper does not give an extensive review of ``geostatistical
inversion'', even for methods based on the linearisation of the flow
equation. For a substantiated analysis of the difficulties posed by the
spatial variability or the heterogeneity of subsurface reservoir
properties for flow and transport modelling, the reader is invited to
refer to \citet{deMarsilyetal2005}. Neuman's feedback
\citeyearpar{Neuman2020} also states the importance of spatial
variability of parameters. A recent non-exhaustive review of methods
for constraining geostatistical estimation by flow physics can be found
among others in \citet{deFouquetetal2023}.

\section{Steady state linear flow} \label{sec2}
\subsection{Linearisation of the flow equations} \label{sec2.1}

For steady state flow without recharge, the diffusivity equations in
two dimensions is written
{\begin{equation}\label{eq1}
\text{div}(T \cdot \text{Grad}\, H)=0
\end{equation}}\noindent
where $T$ represents the transmissivity and $H$ the hydraulic head.

\citet{DelhommedeMarsily2005} observe that Equation~(\ref{eq1}) can be
written
{\begin{equation}\label{eq2}
\Delta H = -\text{Grad}(\ln T)\cdot \text{Grad}\, H
\end{equation}}\noindent
making ``it clear that the gradient of $\ln T$ was intrinsically
involved in the flow equation, regardless of whether $T$ is lognormally
distributed or not''. 

Let now $x=(x_{1}, x_{2})$ denote a point in the plane. In the
simplified case with a macroscopic head oriented in the
$\xoverrightarrow{Ox_{1}}$ direction with constant modulus $J$, and
constant scalar transmissivity $T(x)=T_{0}$, the hydraulic head $h$ is
linear:
{\begin{equation*}
h(x)=h_{0} - Jx_{1}
\end{equation*}}\noindent
where $h_{0}$ is an arbitrary constant (since the hydraulic head is a
potential).

For a transmissivity field with spatial variability, for example when
the transmissivity is lognormal $T(x)=T_{0}\mathrm{e}^{\Theta(x)}$, let
us denote $\Phi$ the associated head perturbation:
$H(x)=h(x)+\Phi(x)$.

When only the first order perturbations are retained,
Equation~(\ref{eq1}) becomes:
{\begin{equation}\label{eq3}
\Delta\Phi = J\partial\Theta/\partial x_1.
\end{equation}}\noindent
The hypothesis that the transmissivity is lognormal is here not
necessary, Equation~(\ref{eq3}) being based on the first order
expansion of the transmissivity in \mbox{Equation}~(\ref{eq1}). 

The initial models of \citet{Mizelletal1982} and
\citet{KitanidisVomvoris1983} linking the covariances of the
perturbations $\Theta$ and $\Phi$ were based on stationarity
assumptions. \citet{Dagan1985} calculated a second-order approximation
of the covariance of the hydraulic head, observing that ``the
first-order approximation is very robust and even for a
log-conductivity variance equal to unity, the second-order correction
of the head variances is smaller than 10\% of the first-order
approximation''.

\citet{Dong1990} wrote the solution of Equation~(\ref{eq3}) as follows,
since transmissivity cannot generally be assumed to be
differentiable:
{\begin{equation}\label{eq4}
J\Theta = \Delta Z\quad\text{and}\quad \Phi = \partial Z/\partial x_1.
\end{equation}}\noindent
This equation shows that head and transmissivity do not have the same
degree of stationarity. Indeed, \citet{Matheron1973} showed that, if
$\Theta$ is a stationary random function of order two, the solution $Z$
of the \mbox{Poisson} equation
{\[
J \Theta = \Delta Z
\]}\noindent
is an Intrinsic Random Function of order 1 (IRF-1), twice
differentiable. More generally, if $\Theta$ is an IRF-$k$, $Z$ is an
IRF-($k+2$).

Matheron also showed that if $Z$ is an IRF-1, the solution $\Phi$ 
of the partial differential equation
{\[
\Phi = \frac{\partial Z}{\partial x_{1}}
\]}\noindent
is an anisotropic intrinsic random function (IRF-0). According to these
results, if the transmissivity is stationary of order two (i.e. mean
and spatial covariance are stationary), the head perturbation $\Phi$ is
an anisotropic IRF-0, and the head is an IRF-0 with a linear drift
(following the first coordinate in this \mbox{particular} case).

The generalised simple and cross covariances between transmissivity and
head perturbations which are solution of Equation~(\ref{eq4})
verify the following \mbox{relationships}, with
$h=(h_{1}, h_{2})$:{\vspace*{-.5pc}}
{\begin{equation}\label{eq5}
J^{2}C_{\Theta}=\Delta^{2}C_{Z};\quad C_{\Phi}=-\frac{\partial^{2}}
{\partial {h}_{1}^{2}} C_{Z};\quad JC_{\Theta\Phi}=\Delta
\frac{\partial}{\partial h_{1}}C_{Z}.
\end{equation}}\noindent
\par\vspace*{.5pc}\noindent
The simple (generalised) covariances are even. The cross-covariance is
odd, and thus equal to zero for $h=0$. The anisotropic variogram of the
head perturbation is generally unbounded.

\citet{Dong1989, Dong1990} and \citet{Roth1995} give the explicit
expressions of the generalised covariance of head perturbations as well
as the cross-covariance between transmissivity and head perturbations,
for usual various models of the variogram of transmissivity. These
expressions allow an estimate of the head consistent with
Equation~(\ref{eq4}). Above all they allow estimating transmissivity by
co-kriging from head and transmissivity measurements, providing a
linearised approximation of the solution of the inverse \mbox{problem}.

In this model, boundary conditions (fixed head, zero flow) are rejected
to infinity and occur via the macroscopic gradient of the head. In the
cokriging these conditions are introduced as data only at the locations
where they are known, and not necessarily all around the modelled
field.

The literature presents several applications of this cokriging for
solving the inverse problem. \citet{RubinDagan1987} consider that the
spatial simple and cross covariances of transmissivity and head depend
on a parameter vector $\theta$ which cannot be determined with
certainty from the data. They develop the analytical expressions for
the conditional covariances when $\theta$ is determined by a maximum
likelihood \mbox{procedure}.

\citet{AhmeddeMarsily1993} perform the cokriging with an approximated
but consistent model (ensuring the positivity of variances) of the
simple and cross-covariances of head and transmissivity, calibrated to
the sample simple and cross-covariances. On a synthetic case
approaching a real one, they show that cokriging improves
transmissivity estimation compared with kriging. Taking into account
the local head gradient rather than the regional one improves the
estimation of the transmissivity.

\citet{Hernandezetal2003} confirm that geostatistical inversion by
co-kriging hydraulic head and log-conductivity improves predictions of
head and flux compared with conditioning on conductivity or head data
alone.

Moreover the expressions~(\ref{eq5}) of simple and cross covariances
between transmissivity and head perturbations can be used for the
direct geostatistical simulation of transmissivity, head and flow when
\mbox{linearisation} assumptions are valid \citep{deFouquet2000,
Mariottietal2000}.

\subsection{Case of two flows} \label{sec2.2}

\citet{ZhangNeuman1996} observed that ``even locally uniform gradients
tend to show seasonal fluctuations in magnitude and direction''. An
example is given in \citet{Verleyetal2003} for the \textit{nappe de la
Beauce}, locally.

Zhang and Neuman examine the effects, on the longitudinal and
transverse dispersion respectively, of these temporal fluctuations in
the magnitude and direction of the mean velocity by developing the auto
and cross-covariances of velocity, head and log conductivity. The
previous first-order approach is generalised considering quasi
steady-state flow at two times. These expressions of auto and cross
covariances of head and log conductivity can also be used for
geostatistical inversion by ``crossing'' two flows, in order to improve
the estimation of transmissivity. 

Let us revisit their developments. Two macroscopically linear flows,
indexed by 1 and 2, are first considered in the same transmissivity
field, following two orthogonal directions $Ox_{1}$ and $Ox_{2}$. To
simplify the presentation, the module $J$ of the macroscopic gradient
is assumed to be the same in the two cases. However this hypothesis can
be relaxed. Based on the relationships~(\ref{eq4}), the solutions
$\Phi_{i}$ are written:
{\begin{equation}\label{eq6}
J\Theta = \Delta Z;\quad \Phi_{i}=\frac{\partial Z}{\partial x_{i}}.
\end{equation}}\noindent
The simple and cross covariances between transmissivity and head
perturbations are deduced from~(\ref{eq6}) 
{\begin{equation}\label{eq7}
\begin{array}{c}
\displaystyle J^{2}C_{\Theta}=\Delta^{2}C_{Z};\quad C_{\phi_{i}}=
-\frac{\partial^2} {\partial h_{i}^{2}}C_{Z};\\
\displaystyle JC_{\Theta \phi_{i}}=\Delta \frac{\partial}{\partial
h_{i}}C_{Z}\quad \text{and}\quad C_{\Phi_{1}\Phi_{2}}=
-\frac{\partial^{2}}{\partial h_{1}\partial h_{2}}C_{Z}.
\end{array}
\end{equation}}\noindent

For a flow with the same macroscopic head gradient $J$ and any angle
$\alpha$ with the first coordinate axe, the head perturbation is
deduced from the linearity of Equation~(\ref{eq6}), putting
{\begin{equation}\label{eq8}
\Phi_{\alpha}=\phi_{1}\cos \alpha +\phi_{2}\sin\alpha.
\end{equation}}\noindent
The simple and cross covariances of transmissivity and head
perturbations are derived from \mbox{Equations}~(\ref{eq7}).

However it is simpler to choose a base such that
$\xoverrightarrow{Ox_{1}} = \xoverrightarrow{Ox_{\alpha}}$. Denoting
$\beta$ the angle between the two macroscopic flows, the expressions of
the simple and cross covariances are simplified:
{\begin{equation}\label{eq9}
\begin{array}{c}
C(\Phi_{1},\Phi_{\beta}) = \cos\beta C(\Phi_{1},\Phi_{1})+\sin\beta
C(\Phi_{1},\Phi_{2});\vspace*{2pt}\\
C(\Theta,\Phi_{\beta}) = \cos\beta C(\Theta, \Phi_{1})+\sin\beta C
(\Theta,\Phi_{2})\vspace*{2pt}\\ 
\begin{array}{r@{~}c@{~}l}
\text{and}\quad C(\Phi_{\beta},\Phi_{\beta}) &=& \cos^{2}\beta
C(\Phi_{1},\Phi_{1})\vspace*{2pt}\\
&& +\,2\cos\beta\sin\beta C(\Phi_{1},\Phi_{2})\vspace*{2pt}\\
&& +\,\sin^{2} \beta C(\Phi_{2},\Phi_{2}).
\end{array}
\end{array}
\end{equation}}\noindent
Because of the linearisation, when the head data points are identical
for the two flows (assumed not to be parallel) the cokriging results
(transmissivity or head) do not depend on the angle $\beta$ between the
two flows.

Always because of the linearity of the equations, the head of a third
flow is expressed as a linear combination of the first two. The
co-kriging system with head data from three flows is therefore
singular, unless measurement errors are introduced in the variogram of
the head data, typically \textit{via} nugget effects.

\section{Example: comparison of cokriging the{\hfill\break}
transmissivity from one or two flows} \label{sec3}
\subsection{Construction of synthetic cases} \label{sec3.1}

Cokriging tests are classically performed as follows:
(i)~geostatistical simulation of a transmissivity field; (ii)~flow
calculation on this field (using a flow simulator, here Visual
Modflow), the fixed macroscopic gradient is imposed via conditions on
the boundaries of a larger field; (iii)~data extraction far enough from
these boundaries; (iv)~kriging and cokriging with data from one or two
flows.

Several areas are specified (Figure~\ref{fig1}):
\begin{itemize}
\item the transmissivity is simulated on a ``simulated area'' of 1600~m
${\times}$ 1600~m;
\item the ``study area'' is 500~m ${\times}$ 500~m from which the data
are taken.
\end{itemize}

The first zone is large enough to fix the boundaries conditions of the
various flows. These boundaries conditions define an ``intermediate
area'' of 1100~m ${\times}$ 1100~m, in which the heads are calculated.
However, in order for the head data to be far enough from the
boundaries, the ``study area'' is smaller than the ``intermediate
area'' (Figure~\ref{fig1}).

\begin{figure}[t!]
\includegraphics{fig01}
\caption{\label{fig1}Definition of the three areas: simulated area for
transmissivity simulation, intermediate area for setting boundary
conditions, and study area. The arrow indicates the main flow direction
for the represented boundary conditions.}
\end{figure}

The log-transmissivity is simulated with a spherical variogram with a
range equal to 70~m and a sill of~0.5. Three flows with a macroscopic
gradient with unit module are calculated
\begin{itemize}
\item North to South (NS)
\item Northwest to Southeast (NW--SE)
\item West to East (W--E)
\end{itemize}

Initial tests have shown that isolated head data reduce only slightly
the estimation variance of transmissivity, while in the vicinity of two
close head data, this estimation variance decreases significantly. This
is consistent with the known observation that ``head variations due to
heterogeneity are small, whereas those of velocities and travel time
are large'' \citep{MatherondeMarsily1980}.

The co-kriging by two flows is therefore examined with a significant
number of close head data (20), cf.~Figure~\ref{fig2}.

\begin{figure}[t!]
\vspace*{4pt}
\includegraphics{fig02}
\vspace*{4pt}
\caption{\label{fig2}Location of the transmissivity ($+$) and head
($\bullet$) data.}
\vspace*{4pt}
\end{figure}

The different estimations of transmissivity are compared according to
the following four criteria, where $N$ denotes the number of meshes
\begin{itemize}
\item the mean absolute error MAE:
{$$
({1}/{N})\sum_{i=1}^{N}|(\log (T_{i}))^{*}-\log (T_{i})|
$$}\noindent
\item the mean absolute standardised error MASE:
{$$
({1}/{N}) \sum_{i=1}^{N}({|(\log (T_{i}))^{*}-\log
(T_{i})|}/{\sigma_i^K})
$$}\noindent
\item the root-mean-square error RMSE:
{$$
\sqrt{({1}/{N})\sum_{i=1}^{N} ((\log (T_{i}))^{*}-\log (T_{i}))^{2}}
$$}\noindent 
\item and the root-mean-square standardised error
RMSSE:
{$$
\sqrt{({1}/{N})\sum_{i=1}^{N}({((\log (T_{i}))^{*}-\log
(T_{i}))^{2}}/{(\sigma_{i}^{K})^{2}})}
$$}\noindent
\end{itemize}

\noindent
where $(\log(T))^{*}$ denotes the estimate of $\log(T)$ and $\sigma^K$
the standard deviation of the estimation (kriging or cokriging) error.
The RMSSE criterion based on the standardised error allow verifying
the expected improvement of accuracy by cokriging.

\subsection{Comparison between one and two flows} \label{sec3.2}

The values in Tables~\ref{tab1} and~\ref{tab2} are given on the decimal
logarithm.

\begin{table*}[t!] %tab1
\tabcolsep=8pt
\caption{\label{tab1}Estimation of log-transmissivity on the 500~m
$\times$ 500~m study area (2601 meshes)\vspace*{4pt}}
\begin{tabular}{cccccc}
\thead
&& MAE & MASE & RMSE & RMSSE\\
\endthead
Kriging && 0.238 & 0.760 & 0.306 & 0.971{\vspace*{6pt}}\\
\morerows{2}{Cokriging 1 flow} & N--S & 0.223 & 0.760 & 0.288 & 0.970\\
& NW--SE & 0.220 & 0.754 & 0.286 & 0.967\\
& W--E & 0.223 & 0.766 & 0.290 & 0.983{\vspace*{6pt}}\\
\morerows{2}{Cokriging 2 flows} & N--S; NW--SE & 0.218 & 0.805 & 0.285 & 1.034\\
& N--S; W--E & 0.214 & 0.790 & 0.281 & 1.015\\
& NW--SE; W--E & 0.218 & 0.804 & 0.284 & 1.028
\botline
\end{tabular}
\tabnote{Values are given on the decimal logarithm.\vspace*{6pt}}
\end{table*}

\begin{table*}[t!] %tab2
\tabcolsep=8pt
\caption{\label{tab2}Estimation of the log-transmissivity in the
heterogeneous case, (a)~in the local area (121 meshes), (b)~on the
extended area (961 meshes) and (c)~in the study area\vspace*{3pt}}
\begin{tabular}{cccccc}
\thead
\multicolumn{6}{l}{(a)~Local area 100~m $\times$ 100~m}\\
&& MAE & MASE & RMSE & RMSSE\\
\endthead
Kriging && 0.916 & 2.667 & 1.077 & 3.138{\vspace*{6pt}}\\
\morerows{2}{Cokriging 1 flow} & N--S & 0.703 & 2.434 & 0.870 & 3.013\\
& NW--SE & 0.612 & 2.184 & 0.758 & 2.716\\
& W--E & 0.365 & 1.281 & 0.450 & 1.576{\vspace*{6pt}}\\
\morerows{2}{Cokriging 2 flows} & N--S; NW--SE & 0.347 & 1.486 & 0.416 & 1.781\\
& N--S; W--E & 0.305 & 1.296 & 0.374 & 1.579\\
& NW--SE; W--E & 0.326 & 1.393 & 0.406 & 1.747\cbotline
\multicolumn{6}{l}{(b)~Extended area 300~m $\times$ 300~m}\cbotline
Kriging && 0.300 & 0.905 & 0.462 & 1.372{\vspace*{6pt}}\\
\morerows{2}{Cokriging 1 flow} & N--S & 0.317 & 1.079 & 0.436 & 1.494\\
& NW--SE & 0.289 & 0.986 & 0.393 & 1.358\\
& W--E & 0.260 & 0.879 & 0.336 & 1.134{\vspace*{6pt}}\\
\morerows{2}{Cokriging 2 flows} & N--S; NW--SE & 0.319 & 1.211 & 0.398 & 1.513\\
& N--S; W--E & 0.316 & 1.193 & 0.408 & 1.532\\
& NW--SE; W--E & 0.247 & 0.938 & 0.321 & 1.225\cbotline
\multicolumn{6}{l}{(c)~Study area 500~m $\times$ 500~m}\cbotline
Kriging && 0.274 & 0.809 & 0.379 & 1.114{\vspace*{6pt}}\\
\morerows{2}{Cokriging 1 flow} & N--S & 0.301 & 0.964 & 0.390 & 1.264\\
& NW--SE & 0.270 & 0.871 & 0.351 & 1.149\\
& W--E & 0.247 & 0.794 & 0.314 & 1.013{\vspace*{6pt}}\\
\morerows{2}{Cokriging 2 flows} & N--S; NW--SE & 0.317 & 1.114 & 0.400 & 1.417\\
& N--S; W--E & 0.301 & 1.059 & 0.382 & 1.352\\
& NW--SE; W--E & 0.239 & 0.840 & 0.305 & 1.078
\botline
\end{tabular}
\tabnote{Values correspond to decimal logarithm.\vspace*{3pt}}
\end{table*}

The mean deviation MAE and the RMSE show that cokriging improves the
estimation accuracy of transmissivity compared with kriging, and even
more when crossing two flows. However, the standardised deviations are
rather stable between kriging and cokriging with one flow, and slightly
increase for two-flows cokriging. The RMSSE are close to unity, which
shows that the uncertainty is correctly quantified, and somewhat
overestimated by kriging or cokriging with one flow.

\subsection{Detection of a transmissivity barrier} \label{sec3.3}

The second test field presents a low transmissivity barrier of
$10^{-6}~\text{m}^{2}{\cdot}\text{s}^{-1}$ compared to a previous
average value of $10^{-4}~\text{m}^{2}{\cdot}\text{s}^{-1}$. The
variogram model of transmissivity remains unchanged, because no
transmissivity data is located in the barrier.

The test is performed with only four transmissivity data. The quality
criteria are first calculated on a small area of 100~m $\times$ 100~m,
centred on the anomaly (Figure~\ref{fig3}), in order to verify if
cokriging helps detect the heterogeneity.

\begin{figure}[t!]
\vspace*{4pt}
\includegraphics{fig03}
\vspace*{4pt}
\caption{\label{fig3}Location of the heterogeneity, represented by an
ellipse. Central square: local area; dashed square: extended area.}
\vspace*{4pt}
\end{figure}

The MAE and the RMSE (Table~\ref{tab2}a) show that the N--S or \mbox{NW--SE} flows, which
are yet little impacted by the transmissivity barrier, improve the
estimation. Logically the W--E flow, which crosses the barrier, informs
the best on the transmissivity field: the mean \mbox{deviation} MAE and
the RMSE decrease by almost 60\%. Note that, with a mean deviation of
0.347, the combination of the N--S and NW--SE flows gives an estimation
accuracy comparable to that of the W--E flow, with a mean deviation of
0.365. The lowest errors (MAE, RMSE) are obtained by cokriging using
the \mbox{N--S} and W--E flows, which both cross the\break barrier.

The standardised deviations are large: the method allows detecting the
heterogeneity, but with only four transmissivity data located outside
it, the estimation accuracy remains poor. However, the increase in the
variance of the transmissivity field was not taken into account in the
modelling.

These results are consistent with Kitanidis remark
\citeyearpar{Kitanidis1995} cited in the introduction according to
which ``the method may work reliably [${\ldots}$] provided that the
variation happens [${\ldots}$] in the direction perpendicular to the
streamlines''.

The results are less satisfactory on a larger scale. In the extended
area of 300~m $\times$ 300~m, bounded by dashes lines on the figure,
the head data improve or on the contrary deteriorate the MAE, depending
on the configurations, but in all cases the RMSE decreases compared
with kriging. The estimations are poorly improved by the head data of
the N--S flow, which is visible on the estimation maps
(Figure~\ref{fig4}): the low transmissivity barrier is detected, but
the surrounding area shows large transmissivities that do not exist.
Once again, the head data from the W--E flow improve greatly the
precision, and to a lesser extent, it is the same for the NW--SE flow.
The MAE and RMSE criteria are minimal when crossing these two flows
(Table~\ref{tab2}b).
In all cases, the RMSSE remains greater than one, for the same reason
as on the local area.\looseness=-1

\begin{figure*}[t!]
\vspace*{4pt}
\includegraphics{fig04}
\vspace*{4pt}
\caption{\label{fig4}Decimal logarithm of Transmissivity field.
(a)~Estimation from N--S flow; (b)~estimation from NW--SE flow;
(c)~estimation from the two N--S/NW--SE flows; (d)~``real'' field.}
\vspace*{6pt}
\end{figure*}

The results remain similar on the study zone (Table~\ref{tab2}c): with
head data from the N--S flow, the MAE and the RMSE increase compared
with kriging and they decrease sharply with head data from the W--E
flow only. Once again, MAE and RMSE are minimal by crossing W--E and
NW--SE flows. For all configurations, the RMSSE is lower than for the
homologous case for the extended area. 

However on these three heterogeneous areas, the available data do not
allow estimating the transmissivity accurately.

\subsection{Interpreting the results} \label{sec3.4}

The improvement of the precision of the transmissivity with head data
from two flows can be explained by returning to the stream lines. It is
known \citep{EmsellemdeMarsily1971, deMarsilyetal2000} that the
solution of the inverse problem for the steady-state flow equation,
including a recharge or pumping, is written
{\[
\ln \frac{T(s)}{T_{0}} = -\int_{h_{0}}^{h(s)}\frac{\Delta h}{|\nabla
h|^{2}}\,\mathrm{d}h,
\]}\noindent
where the integration is made along a stream-tube, $T_{0}$ is the
transmissivity and $h_{0}$ the hydraulic head at point $s_{0}$ of the
tube, and $T(s)$ the transmissivity at point $s$ where the head $h(s)$
is measured.

The transmissivity at any point of this stream-tube can thus be deduced
from the transmissivity at $s_{0}$. Formally, knowing the
transmissivity along a line intersecting all the stream-tubes allows
deducing the transmissivity everywhere in the field. In the case of two
flows with nonparallel stream lines, knowing the transmissivity in a
single point allows deducing the transmissivity everywhere in the
field; see for example \citet{Castelier1995}.

In practice, the head is not known exhaustively but estimated, and the
same for the stream-tubes. The estimate of transmissivity is therefore
uncertain.

\subsection{Practical conclusions from these synthetic{\hfill\break}
cases} \label{sec3.5}

This synthetic example leads to practical conclusions~on estimating
transmissivity field:
\begin{enumerate}[(i)]
\item sufficiently close head measurements allow estimating the head
gradient, which is more useful for the accuracy than isolated head
values;
\item ``crossing'' two steady-state flows improves the accuracy of
transmissivity estimation;
\item this flow crossing is useful for detecting heterogeneities of a
transmissivity field. However, transmissivity data have to be numerous
enough for precisely delineating the heterogeneities; 
\item since the accuracy of the estimation depends on the main flow
direction, the same applies to the ``optimal'' implantation of the
data. 
\end{enumerate}

\section{Conclusive remarks} \label{sec4}
\subsection{Towards more complex flow conditions} \label{sec4.1}

In many cases, the hypothesis of a macroscopically linear
two-dimensional flow is unsuitable. 

\citet{Rothetal1998} developed a flexible and general approach in order
to avoid excessive simplifications. The principle is to come back to
conventional hydrogeological modelling, using a flow simulator. A set
of geostatistical (non-conditional) simulated transmissivity fields is
used as input of a flow simulator. The set of associated output results
(head, concentration) is used to calculate the explicitly nonstationary
covariance $C(x,x')$ or variogram $\gamma(x,x')$ of the variable of
interest. These nonstationary ``numerical'' simple and cross
covariances or variograms of transmissivity and head are used to
estimate by cokriging the transmissivity from transmissivity and head
data. However an anamorphosis or at least a logarithmic transformation
allows avoiding negative transmissivity estimates.
\citet{Pannecoucke2020} generalised this ``numerical variogram'' to
transient flow, for mapping radionuclide activity.

To account for the non-linearity between concentration and hydraulic
conductivity field, \citet{SchwedeCirpka2010} show that a Monte-Carlo
approach is preferable for conditioning hydraulic conductivity fields
by steady-state concentration measurements.\looseness=1

\subsection{Back transformation of estimated{\hfill\break}
log-transmissivity} \label{sec4.2}

According to \citet{deMarsilyetal2005} ``One also needs to apply
Kriging with some rigor on the log-transformed transmissivity values,
in order to estimate geometric mean values and not arithmetic means.
Back-transforming the kriged $\ln T$ into $T$ values must also be done
correctly, i.e. simply as $T=\exp (\ln T)$ without any additional
factor using the variance of the estimation error to supposedly correct
an assumed estimation of the median rather than of the mean, as is
sometimes erroneously done''.

Another question is that of the support of the estimate: is kriging
``punctual'', that is to say that the transmissivity has to be
estimated on the same support as the transmissivity data, or is a
larger support chosen, in which case an equivalent transmissivity is
sought? However, except in the very special case of a stratified medium
with a flow parallel to the stratification, the permeability is not
additive. The equivalent block-permeability lies within the \mbox{harmonic}
and the arithmetic mean of sub-block permeabilities
\citep{Matheron1993, DelhommedeMarsily2005}. In the special case of
lognormal transmissivity with isotropic variogram, and with some
additional hypotheses, \citet{Matheron1967} has shown that the
equivalent transmissivity is the geometric mean. Coming back to the
estimated transmissivity simply consists thus to exponentiate the
estimated log-transmissivity.

In practice, the regularisation similarly has to be carried out on the
log-transmissivities, and the result is simply exponentied. The bounds
of the uncertainty interval can be calculated in the same way from the
bounds of this interval on logarithmic estimates. 

In a more thorough approach, tensor transmissivity fields should be
estimated.

\subsection{Biased flow predictions computed on{\hfill\break} estimated
transmissivity fields} \label{sec4.3}

However it is known that the (co)kriging ``smoothes'' map of the
estimated values, whose spatial variability therefore differs from that
of the real values. The non-linear calculations performed on an
estimated transmissivity field (by (co)kriging or any other estimator,
linear or not) is biased. A simulation based approach is then required.

The cokriging of log-transmissivity by log-transmissivity and head
measurements however allows the direct conditioning of
log-transmissivity simulations to these data. When the validity of the
linearisation of the flow equation is not verified, especially when the
head gradient cannot be assumed to be constant at the macroscopic
scale, the simple and crossed ``numerical variogram'' between
log-transmissivity and head (perturbations) provides a consistent
bivariate model, which is compatible with the boundary conditions where
these are known 
\citep{SchwedeCirpka2010, Pannecouckeetal2020}. However, the results
shown here in the simplified case of the linearized flow equation
remain useful in the general case: ``crossing'' two quasi steady-state
flows (two seasonal situations for example) should improve
characterising transmissivity heterogeneities.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds
from any organization that could benefit from this article, and have
declared no affiliations other than their research organizations.

\section*{Acknowledgments}

The authors warmly thank Jean-Pierre Delhomme for initiating and
contributing to the study of flow crossing. They also thank Marco De
Lucia for his initial contribution to review the theoretical part of
this study and the reviewers for their comments that helped to better
highlight the results.

\back{}

\bibliographystyle{crgeos}
\bibliography{crgeos20230417}
\refinput{crgeos20230417-reference.tex}

\end{document}
