GSTools v1.3: a toolbox for geostatistical modelling in Python

Geostatistics as a subfield of statistics accounts for the spatial correlations encountered in many applications of, for example, earth sciences. Valuable information can be extracted from these correlations, also helping to address the often encountered burden of data scarcity. Despite the value of additional data, the use of geostatistics still falls short of its potential. This problem is often connected to the lack of user-friendly software hampering the use and application of geostatistics. We therefore present GSTools, a Python-based software suite for solving a wide range of geostatistical problems. We chose Python due to its unique balance between usability, flexibility, and efficiency and due to its adoption in the scientific community. GSTools provides methods for generating random fields; it can perform kriging, variogram estimation and much more. We demonstrate its abilities by virtue of a series of example applications detailing their use.

1 Introduction

Geostatistics emerged as a distinct branch of statistics in the early 1950s through the pioneering work by Krige (1951). Krige's goal of estimating the abundance of mineral resources led him to develop some of the first methods, but it was the French mathematician Georges Matheron who developed the mathematical foundations (Matheron, 1962). Today, geostatistics is applied in fields like geology (Hohn, 1999), hydrogeology (Kitanidis, 2008), hydrology or soil sciences (Goovaerts, 1999), meteorology (Cecinati et al., 2017), ecology (Rossi et al., 1992; Sales et al., 2007), oceanography (Monestiez et al., 2004), and epidemiology (Schüler et al., 2021), and a large number of textbooks make the theory available to practitioners (Pyrcz and Deutsch, 2014; Rubin, 2003; Diggle and Ribeiro, 2007; Kitanidis, 2008; Banerjee et al., 2014).

Yet, the rate of adoption of geostatistics by practitioners has been slow and uneven (Zhang and Zhang, 2004; Rajaram, 2016). One reason is the perceived lack of ready-made geostatistical software (Zhang and Zhang, 2004; Neuman, 2004; Winter, 2004; Rajaram, 2016; Cirpka and Valocchi, 2016; Fiori et al., 2016). Although a decent number of geostatistical software solutions are available (Bellin and Rubin, 1996; Deutsch and Journel, 1997; Brouste et al., 2008; Rubin et al., 2010; Pebesma, 2004; Savoy et al., 2017; Heße et al., 2014; Vrugt, 2016), user-friendliness and licensing can hamper their adoption as pointed out by Rubin et al. (2018).

Addressing these challenges, we present GSTools – an extensive Python package for geostatistical analysis (Müller and Schüler, 2021). To the best of our knowledge, no open-source Python package that provides such a comprehensive collection of random field generation, forward modelling, kriging and data analysis currently exists.

We believe that the choice of Python has the potential to address several of the challenges for geostatistical applications. First, a script language like Python allows for striking a balance between ease of use (as provided by graphical user interfaces) and flexibility (as provided by command-line-based tools). The presence of a graphical user interface (GUI) is sometimes seen as an indication of usability (Remy, 2005; Rubin et al., 2018). However, a GUI does not necessarily make software more user-friendly, almost always limits flexibility, increases the programming effort, and makes reproducibility of results and workflows harder or even unfeasible (Queiroz et al., 2017). Also note that in recent years powerful tools emerged in data science to represent data along with code snippets, documentation and results in interactive computational notebooks provided by Jupyter (Perkel, 2018; Beg et al., 2021). Second, Python is known as a glue language, being able to combine independent software solutions to achieve complex workflows. This is particularly important since geostatistics often relies on ready-made solvers for data generation or partial differential equation (PDE)-based model solvers. Third, Python is a simple yet powerful language with an increasing user base and community support for scientific computing and data analysis. It thus has a wide appeal and excellent prospects for the foreseeable future. This guarantees that engineers and scientists with only a moderate background in computer science are able to apply the toolbox and to make the necessary application-specific adjustments. Finally, the licensing should be as permissive as possible to guarantee adoption and even further development by interested users.

We introduce GSTools and present its main features with a general overview of its functionality and abilities in Sect. 2. We focus on the covariance model, field generation kriging and variogram estimation. In Sect. 3, we discuss the wider context of GSTools, namely a number of Python packages connected with GSTools which can be used to seamlessly model geostatistical workflows. Section 4 presents a number of workflows to showcase the abilities of GSTools and demonstrate its usage. We close out with a short summary of the main advantages of GSTools and concluding remarks.

2 GSTools features

2.1 Covariance models and variography

The powerful CovModel class represents covariance and variogram models. Methods provided by this class are the basis for most of the functionality of GSTools, such as variography, spatial random field generation and kriging.

2.1.1 Covariance models

GSTools implements a CovModel class to define covariance models of weakly stationary (spatial) processes. Here, weak stationarity means that the associated semi-variogram is bounded, since we assume a constant mean and a finite variance. To approximate unbounded variograms such as the power-law model (Webster and Oliver, 2007), we provide a set of truncated power-law models following Di Federico and Neuman (1997). The internal representation of a (semi-)variogram γ is given by

$$\begin{array}{}\text{(1)}& \mathit{\gamma }\left(r\right)={\mathit{\sigma }}^{\mathrm{2}}\cdot \left(\mathrm{1}-\mathrm{cor}\left(s\cdot {\displaystyle \frac{r}{\mathrm{\ell }}}\right)\right)+n,\end{array}$$

where r is the (isotropic) lag distance, ℓ is the (main) correlation length, s is a rescaling factor to adjust model representation (default is 1), σ² is the variance or partial sill, n is the nugget or sub-scale variance, and cor(h) is the model-defining, normalized correlation function depending on the non-dimensional distance $h=s\cdot \frac{r}{\mathrm{\ell }}$.

The associated covariance and correlation functions are given by the following:

$$\begin{array}{}\text{(2)}& {\displaystyle }C\left(r\right)& {\displaystyle }={\mathit{\sigma }}^{\mathrm{2}}\cdot \mathrm{cor}\left(s\cdot {\displaystyle \frac{r}{\mathrm{\ell }}}\right),\text{(3)}& {\displaystyle }\mathit{\rho }\left(r\right)& {\displaystyle }=\mathrm{cor}\left(s\cdot {\displaystyle \frac{r}{\mathrm{\ell }}}\right).\end{array}$$

Note that covariance and correlation are neglecting the nugget effect at the origin. Thus, the variance is interpreted as the variation above the nugget, which is sometimes referred to as the partial sill of the semi-variogram or the correlated variability (Rubin, 2003). Consequently, the sill or limit of the semi-variogram is calculated as the sum of variance and nugget.

Figure 1 Initialization of an exponential covariance model given by $\mathrm{cor}\left(h\right)=\exp (-h)$ (Rubin, 2003). Note that the rescaling factor is 1 by default. The right panel shows the plot of the variogram, covariance and correlation functions of the model, which can be created with convenience plotting methods.

The (semi-)variogram, covariance and correlation functions of a model are accessible through model.variogram, model.covariance and model.correlation, respectively. Every covariance model is defined by at least six parameters: dimension dim, variance var, main length scale len_scale, rescale factor rescale, anisotropy ratios anis and rotation angles angles, with the last two being dimension dependent. Figure 1 shows an example code for instantiating an exponential model and the resulting model functions exemplifying the parameters. Table 1 provides an overview of the predefined models in GSTools.

Webster and Oliver (2007)Webster and Oliver (2007)Wackernagel (2003)Rasmussen and Williams (2005)Rasmussen and Williams (2005)Chilès and Delfiner (2012)Webster and Oliver (2007)Webster and Oliver (2007)Webster and Oliver (2007)Matérn (1960)Matérn (1960)Chilès and Delfiner (2012)Wendland (1995)Di Federico and Neuman (1997)Di Federico and Neuman (1997)Müller et al. (2021a)

Table 1 Predefined covariance models in GSTools.

Formulas including the subscript (h<1) are piecewise-defined functions being constantly zero for h≥1. Here, h is the non-dimensional distance, d is the dimension, Γ(x) is the gamma function, K_ν(x) is the modified Bessel function of the second kind, J_ν(x) is the Bessel function of the first kind, ${}_{\mathrm{2}}{F}_{\mathrm{1}}(a,b,c,x)$ is the ordinary hyper-geometric function and E_ν(x) is the exponential integral function (Abramowitz et al., 1972). All other variables are shape parameters of the respective models.

In addition to the predefined covariance models, users can specify their own model functions by providing a normalized correlation function. Figure 2 shows a reimplementation of the exponential model in only three lines of code.

Figure 2 Initialization of a user-defined exponential covariance model. The only thing that needs to be defined is the normalized correlation function cor.

Figure 3 Spatial covariance structure of an anisotropic exponential model in 3D plotted with the built-in interactive routines of GSTools. The example shows an eighth turn on the xy plane with anisotropy factors $(\mathrm{1}/\mathrm{2},\mathrm{1}/\mathrm{4})$. Rotation angles are given in radians.

The dimension-dependent spectrum of an isotropic covariance model can be called with model.spectrum. It is directly calculated from the covariance function by the following:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}S\left(\mathit{k}\right)& ={\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\right)}^d\cdot \underset{{\mathbb{R}}^d}{\int }C\left(\right|\mathit{r}\left|\right)\cdot e^{i\cdot \mathit{k}\cdot \mathit{r}}\mathrm{d}\mathit{r}\\ & ={\displaystyle \frac{\left|\mathit{k}\right|}{\left(\mathrm{2}\mathit{\pi }\right|\mathit{k}|{)}^{\frac{d}{\mathrm{2}}}}}\cdot {\mathcal{H}}_{\frac{d}{\mathrm{2}}-\mathrm{1}}\left\{r^{\frac{d}{\mathrm{2}}-\mathrm{1}}C\left(\cdot \right)\right\}\left(\right|\mathit{k}\left|\right).\end{array}\end{array}$$

Here, $r=\left|\mathit{r}\right|$ and $k=\left|\mathit{k}\right|$ are the norms of the corresponding vectors, and ℋ is the Hankel transform, which provides a mathematically self-contained and numerically robust formulation of the radially symmetric Fourier transformation. GSTools makes use of an implementation of ℋ provided by the Python package hankel (Murray and Poulin, 2019; Ogata, 2005). For models with a known analytical solution, GSTools uses them for improved computations.

A prerequisite for kriging or random field generation is that the applied covariance function is positive (semi-)definite. That can be checked through the spectral density which is derived by the following (note that S only depends on the norm of k):

$$\begin{array}{}\text{(5)}& {\displaystyle }E\left(k\right)={\displaystyle \frac{S\left(k\right)}{{\mathit{\sigma }}^{\mathrm{2}}}}=k(\mathrm{2}\mathit{\pi }k{)}^{-\frac{d}{\mathrm{2}}}\cdot {\mathcal{H}}_{\frac{d}{\mathrm{2}}-\mathrm{1}}\left\{r^{\frac{d}{\mathrm{2}}-\mathrm{1}}\mathit{\rho }\left(\cdot \right)\right\}(k).\end{array}$$

From Bochner's theorem (Rudin, 1990), it follows that the spectral density is a probability density function if and only if the underlying covariance function is positive (semi-)definite, which all predefined models in GSTools satisfy. As a consequence, the error variance during kriging is always non-negative.

2.1.2 Anisotropy and rotation

Variograms are typically defined based on the lag distance r, resulting in an isotropic model. However, many natural processes involve anisotropy with varying correlation ranges in different (orthogonal) directions. An example is hydraulic conductivity, where anisotropy typically arises from the geologic stratification. The implementation of anisotropy in GSTools is based on the non-dimensional distance (Rubin, 2003):

$$\begin{array}{}\text{(6)}& {\displaystyle }h=\sqrt{\sum _{i=\mathrm{1}}^d{\left({\displaystyle \frac{r_i}{{\mathrm{\ell }}_i}}\right)}^{\mathrm{2}}}={\displaystyle \frac{s}{\mathrm{\ell }}}\sqrt{\sum _{i=\mathrm{1}}^d{\left({\displaystyle \frac{r_i}{e_i}}\right)}^{\mathrm{2}}}=s\cdot {\displaystyle \frac{\stackrel{\mathrm{̃}}{r}}{\mathrm{\ell }}},\end{array}$$

where $\mathrm{\ell }=s\cdot {\mathrm{\ell }}_{\mathrm{1}}$ is the main length scale incorporating the rescale factor s, $e_i=\frac{{\mathrm{\ell }}_i}{{\mathrm{\ell }}_{\mathrm{1}}}$ represents the anisotropy ratios and $\mathit{r}=(r_{\mathrm{1}},r_{\mathrm{2}},\mathrm{\dots })$ represents the distances along the main axis of correlation resulting in the isotropic distance $\stackrel{\mathrm{̃}}{r}$. Consequently, GSTools uses a main length scale, a set of anisotropy ratios and a set of rotation angles to fully describe an anisotropic model.

In practice, the main directions of correlation do not necessarily follow the principal axis. The CovModel accounts for rotation through rotation angles, where their number m depends on the dimension d: $m=\frac{d\cdot (d-\mathrm{1})}{\mathrm{2}}$. In two dimensions, rotation is fully described by a single angle for rotation in the xy plane, and in three dimensions three angles are applied to the xy plane, xz plane and yz plane, respectively. In three dimensions these are often referred to as Tait–Bryan angles yaw, pitch and roll (Goldstein, 1980); see Fig. 3 for an example.

Figure 4 Initialization of a Yadrenko covariance model. We use the Earth's radius as the rescaling factor to have a meaningful length scale. The routine vario_yadrenko still depends on the central angle given in radians.

Figure 5 Estimating an empirical variogram of synthetic unstructured data and fitting an exponential model. The number of bins was calculated to be 21 with a maximum bin distance of ca. 45.

One unique feature of GSTools is the support of arbitrary dimensionality in all provided routines. For rotation in higher dimensions, we apply the following scheme: the first angles coincide with those of the next lower dimension, and the added d−1 angles describe rotations in the planes of the added dimension (in 3D: xz and yz). Thus, there are 6 rotation angles in 4D, 10 in 5D, etc. Rotation in higher dimensions is only relevant for spatio-temporal modelling with three spatial dimensions and application to other fields of research with high-dimension data. The scheme was chosen for metric spatio-temporal models to account for spatial anisotropy in a similar way as a simple spatial model.

Rotation in the x_ix_j plane is described by the matrix $\mathbf{G}(\mathit{\alpha },[i,j\left]\right)\in {\mathbb{R}}^{d\times d}$.

$$\begin{array}{}\text{(7)}& {\displaystyle }\mathbf{G}(\mathit{\alpha },[i,j]{)}_{kl}=\left\{\begin{array}{ll}\cos \mathit{\theta }& k=l=i,j\\ \sin \mathit{\theta }& k=i,l=j\\ -\sin \mathit{\theta }& k=j,l=i\\ \mathrm{1}& k=l\ne i,j\\ \mathrm{0}& \text{else}\end{array}\right.\end{array}$$

The order of rotating planes is determined by the described scheme, i.e. $I_{\mathrm{1}}=[\mathrm{1},\mathrm{2}]$ (xy plane), $I_{\mathrm{2}}=[\mathrm{1},\mathrm{3}]$, (xz plane) $I_{\mathrm{3}}=[\mathrm{2},\mathrm{3}]$ (yz plane), etc. These values define a rotation matrix Rot to transform principal axes to the directions of correlation and the back-rotation matrix $\mathbf{DeRot}={\mathbf{Rot}}^{-\mathrm{1}}$ for the inverse:

$$\begin{array}{ll}{\displaystyle }\mathbf{Rot}& {\displaystyle }=\stackrel{\mathit{\curvearrowleft }}{\prod _{i=\mathrm{1}}^m}\mathbf{\text{G}}\left(\right(-\mathrm{1}{)}^{i-\mathrm{1}}{\mathit{\alpha }}_i,I_i)\\ \text{(8)}& {\displaystyle }& {\displaystyle }=\mathbf{\text{G}}\left(\right(-\mathrm{1}{)}^{m-\mathrm{1}}{\mathit{\alpha }}_m,I_m)\cdot \mathrm{\dots }\cdot \mathbf{\text{G}}({\mathit{\alpha }}_{\mathrm{1}},I_{\mathrm{1}}).\end{array}$$

The alternating signs of the rotation angles $(-\mathrm{1}{)}^{i-\mathrm{1}}{\mathit{\alpha }}_i$ were chosen to match Tait–Bryan angles in 3D.

Figure 6 Estimation of directional variograms for given main axes: the code snippet shows the setup for estimating and fitting the variogram to an anisotropic field. The graphs show the main axes of the rotated model and the fitting results. Plotting commands have been omitted.​​​​​​​

Figure 7 Estimating an empirical variogram (bottom left) of synthetic unstructured data (top left) after Box–Cox normalization of skewed input values. Panels on the right show the histogram of the data values before (top) and after (bottom) the normalization. For demonstration purposes, a Matern model was fitted to the empirical variogram. Plotting commands have been omitted.

For applying or removing anisotropy, we define the isotropify matrix $\mathbf{Iso}=\mathrm{diag}(e_{\mathrm{1}}^{-\mathrm{1}},e_{\mathrm{2}}^{-\mathrm{1}},\mathrm{\dots })$ and anisotropify matrix $\mathbf{AnIso}={\mathbf{Iso}}^{-\mathrm{1}}$. Combining these two types of matrices allows us to isometrize (i.e. isotropify and derotate) and anisometrize (i.e. rotate and anisotropify) spatial points via the following:

$$\begin{array}{}\text{(9)}& {\displaystyle }\mathbf{Isom}& {\displaystyle }=\mathbf{Iso}\cdot \mathbf{DeRot},\text{(10)}& {\displaystyle }\mathbf{AnIsom}& {\displaystyle }=\mathbf{Rot}\cdot \mathbf{AnIso}.\end{array}$$

GSTools provides the routine CovModel.isometrize to convert spatial positions to their derotated and isotropic counterparts as required by Eq. (6) and the routine CovModel.anisometrize to invert this:

$$\begin{array}{}\text{(11)}& {\displaystyle }{\mathit{x}}_{\mathrm{isom}}& {\displaystyle }=\mathbf{Isom}\cdot \mathit{x},\text{(12)}& {\displaystyle }{\mathit{x}}_{\mathrm{anisom}}& {\displaystyle }=\mathbf{AnIsom}\cdot \mathit{x}.\end{array}$$

Figure 8 Comparison of parametric normalizers in GSTools.

2.1.3 Geographical coordinates

Earth's surface is a non-Euclidean manifold, and all large-scale, geographically referenced data will necessarily reflect that. We deal with the non-Euclidean nature of this kind of data by assuming the Earth to be a perfect sphere and then using the fact that the distance between two points $p_{\mathrm{1}}=({\mathit{\varphi }}_{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{1}})$ and $p_{\mathrm{2}}=({\mathit{\varphi }}_{\mathrm{2}},{\mathit{\lambda }}_{\mathrm{2}})$ is given by their latitude (ϕ) and longitude (λ) and can be described by a central angle calculated from the great circle distance:

$$\begin{array}{}\text{(13)}& \mathit{\zeta }(p_{\mathrm{1}},p_{\mathrm{2}})=\arccos \left(\sin {\mathit{\varphi }}_{\mathrm{1}}\sin {\mathit{\varphi }}_{\mathrm{2}}+\cos {\mathit{\varphi }}_{\mathrm{1}}\cos {\mathit{\varphi }}_{\mathrm{2}}\cos \left(\mathrm{\Delta }\mathit{\lambda }\right)\right)\end{array}$$

A huge family of valid models on the sphere can be derived from 3D models by inserting the chordal distance which results in the associated Yadrenko covariance model C_Y (Lantuéjoul et al., 2019):

$$\begin{array}{}\text{(14)}& {\displaystyle }{C}_{\mathrm{Y}}\left(\mathit{\zeta }\right)=C_{\mathrm{3}\mathrm{D}}\left(\mathrm{2}\cdot \sin \left({\displaystyle \frac{\mathit{\zeta }}{\mathrm{2}}}\right)\right).\end{array}$$

The underlying manifold introduces new restrictions for covariance models to be positive definite. The manifold structure of the sphere only allows for isotropic models. For small-scale applications, it is valid to assume anisotropy. An appropriate adaption is the use of a 2D projection like Gauss–Krüger coordinates. We provide Yadrenko models as a unified representation for non-Euclidian coordinates since they facilitate all presented models to be used with geographical coordinates as demonstrated in Fig. 4.

2.1.4 Empirical variogram, data preparation and model fitting

The empirical variogram is an important tool for analysing spatially correlated data. It is estimated with the subpackage gstools.variogram, which provides two estimators for the empirical variogram: Matheron and Cressie (Webster and Oliver, 2007). The default Matheron's estimator for a variogram γ of a spatial random field U is given by

$$\begin{array}{}\text{(15)}& \mathit{\gamma }\left(r\right)=\frac{\mathrm{1}}{\mathrm{2}}\cdot \left|M\right(r){|}^{-\mathrm{1}}\sum _{M\left(r\right)}{\left(U\left({\mathit{x}}_i\right)-U\left({\mathit{x}}_j\right)\right)}^{\mathrm{2}},\end{array}$$

where M(r) is the set of all pairwise spatial random field points, separated by the distance r and a certain tolerance ε>0.

Cressie's estimator, which is more robust to outliers, is given by

$$\begin{array}{}\text{(16)}& \mathit{\gamma }\left(r\right)={\displaystyle \frac{\frac{\mathrm{1}}{\mathrm{2}}\cdot {\left(\left|M\right(r){|}^{-\mathrm{1}}{\sum }_{M\left(r\right)}\sqrt{\left|U\right({\mathit{x}}_i)-U({\mathit{x}}_j\left)\right|}\right)}^{\mathrm{4}}}{\mathrm{0.457}+\mathrm{0.494}/\left|M\right(r\left)\right|+\mathrm{0.045}/\left|M\right(r){|}^{\mathrm{2}}}}.\end{array}$$

Both estimators require predefined bins M(r) to group the pairwise point distances of the given field. GSTools provides a standard binning routine, where the maximal bin width is set to one-third of the diameter of the containing box of the field, the number of bins is determined by Sturges' rule (Sturges, 1926) and all bins have equal width. Figure 5 provides an example of the variogram estimation of an unstructured spatial random field with automatic binning.

GSTools accounts for anisotropy by providing estimating routines for directional variograms along a specified direction with a certain angle tolerance and bandwidth. When providing orthogonal axes, it is possible to fit a theoretical model and its anisotropy ratios as shown in Fig. 6. Determining the main rotation axes from given data, however, is up to the user and beyond the scope of the presented GSTools version.

Figure 9 Comparison of simple, ordinary and universal kriging. All three routines have a similar setup, where simple kriging needs an estimated mean and where universal kriging needs additional drift functions. Plotting commands have been omitted.

Figure 10 A simple setup for linear regression kriging. Although the interpolation coincides with a piecewise linear function, we gain information about the error variances between the conditioning points as shown in the right plot.

Field data often do not follow a normal distribution, which is a crucial assumption for variogram estimation. For example, transmissivity is usually assumed to be log-normally distributed (Dagan, 1989), while rainfall data are normalized by applying the Box–Cox transformation (Cecinati et al., 2017). GSTools provides a set of normalizers based on power transforms, which can be fitted to a given data set using a maximum likelihood approach (Eliason, 1993): LogNormal, BoxCox (Box and Cox, 1964), YeoJohnson (Yeo and Johnson, 2000), Modulus (John and Draper, 1980) and Manly (Manly, 1976). An example application is shown in Fig. 7, and a comparison of all provided normalizers can be seen in Fig. 8.

GSTools also provides routines to detrend data. For example, temperature could decrease with elevation or conductivity could decrease with depth. Another application is analysing spatial correlation of residuals after application of a regression model to the data. All routines dealing with data have the keywords trend, normalizer and mean, where the last keyword describes the mean of the normalized data.

2.2 Kriging, random fields and conditioned random fields

2.2.1 Kriging

The subpackage gstools.krige provides routines for Gaussian process regression, also known as kriging, which is a method of data interpolation based on predefined covariance models (Wackernagel, 2003). Kriging aims to derive the value of a field z at some point x₀, when there are fixed conditioning values z(x₁)…z(x_n) at given points x₁…x_n. The resulting value z₀ at x₀ is calculated as a weighted mean $z_{\mathrm{0}}={\sum }_{i=\mathrm{1}}^n{w}_i\cdot z_i$, where the weights, $\mathit{w}=(w_{\mathrm{1}},\mathrm{\dots },w_n)$, are determined by the specific kriging routine.

We provide multiple kriging routines derived from the Krige class. (i) Simple: the data are interpolated with a given mean value for the kriging field. (ii) Ordinary: the mean of the resulting field is unknown and estimated alongside the interpolation (unbiasedness). (iii) Universal: in addition to ordinary kriging, one can provide drift functions $f_{\mathrm{1}},\mathrm{\dots },f_{\mathrm{k}}$. (iv) ExtDrift: like universal kriging but the drift is provided by an external source.

The advantage of using the general Krige class is the combination of all described features, such as for instance using universal kriging with a functional drift together with additional external drifts. A typical scenario is a temperature interpolation with an assumed north–south drift (functional drift) and a linear correlation to altitude (external drift).

Figure 11 Generation of a structured random field following a Gaussian variogram.

Figure 12 An example for an ensemble of 1D random fields conditioned to five measurements (dots). Plotting commands have been omitted.

Since all variogram models in GSTools assume weak stationarity, the kriging system is always built on the associated covariance function:

$$\begin{array}{}\text{(17)}& \textcolor{red}{}\end{array}$$

with $\mathbf{C}=\left(C\right({\mathit{x}}_i,{\mathit{x}}_j){)}_{ij=\mathrm{1}\mathrm{\dots }n}$ being the covariance matrix, depending on the conditioning points and the given model. $C_{\mathrm{0}}=\left(C\right({\mathit{x}}_i,{\mathit{x}}_{\mathrm{0}}){)}_{i=\mathrm{1}\mathrm{\dots }n}^{\mathrm{T}}$ is the covariance vector for the target point x₀. $\mathbf{F}=\left(f_j\right({\mathit{x}}_i){)}_{i=\mathrm{1}\mathrm{\dots }n,j=\mathrm{1}\mathrm{\dots }k}$ is a submatrix containing the functional drift values at the conditioning points and $F_{\mathrm{0}}=\left(f_i\right({\mathit{x}}_{\mathrm{0}}){)}_{i=\mathrm{1}\mathrm{\dots }k}^{\mathrm{T}}$ at the target point, where k is the number of functional drifts. $\mathbf{E}=(e_{ij}{)}_{i=\mathrm{1}\mathrm{\dots }n,j=\mathrm{1}\mathrm{\dots }l}$ is a submatrix containing the external drift values at the conditioning points and $E_{\mathrm{0}}=(e_{i\mathrm{0}}{)}_{i=\mathrm{1}\mathrm{\dots }l}^{\mathrm{T}}$​​​​​​​ at the target point, where l is the number of external drifts. The parameters μ, $\mathit{\varphi }=({\mathit{\varphi }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\varphi }}_k{)}^{\mathrm{T}}$ and $\mathit{\psi }=({\mathit{\psi }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\psi }}_l{)}^{\mathrm{T}}$ are Lagrange multipliers for the unbiased condition, the functional drifts and the external drifts, respectively. The vector 1 and its Lagrange multiplier μ are given in brackets since their appearance depends on whether the system should be unbiased or not (ordinary vs. simple kriging). Note that the number of functional drifts k and external drifts l can be zero, depending on whether they are given or not.

GSTools also provides the possibility to incorporate measurement errors variances ${\mathit{\sigma }}_i^{\mathrm{2}}$ for each conditioning point by adjusting the covariance matrix (Wackernagel, 2003):

$$\begin{array}{ll}{\displaystyle }\stackrel{\mathrm{̃}}{\mathbf{C}}& {\displaystyle }=\mathbf{C}+\mathrm{diag}({\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}},\mathrm{\dots },{\mathit{\sigma }}_n^{\mathrm{2}})\\ \text{(18)}& {\displaystyle }& {\displaystyle }=\left[\begin{array}{ccc}C({\mathit{x}}_{\mathrm{1}},{\mathit{x}}_{\mathrm{1}})+{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}& \mathrm{\dots }& C({\mathit{x}}_{\mathrm{1}},{\mathit{x}}_n)\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ C({\mathit{x}}_n,{\mathit{x}}_{\mathrm{1}})& \mathrm{\dots }& C({\mathit{x}}_n,{\mathit{x}}_n)+{\mathit{\sigma }}_n^{\mathrm{2}}\end{array}\right].\end{array}$$

By default, the measurement error variances, ${\mathit{\sigma }}_i^{\mathrm{2}}$, are set to the model nugget. In order to get numerically stable results, we solve the kriging system with the pseudo-inverse matrix, which has the advantage that redundant data or multiple measurements at the same location are averaged out in the resulting field (Mohammadi et al., 2017).

One last feature is the capability of kriging the mean (Wackernagel, 2003), which allows for deriving the mean value estimated during ordinary kriging or estimating the mean drift determined from given functional and/or external drift terms as shown in Fig. 9. A minimal example for regression kriging is shown in Fig. 10.

Figure 13 Generation of a structured incompressible random vector field with exponential variogram.

Figure 14 Example of a log-transformed binary field with the low values being connected by applying the zinnharvey, binary and lognormal transformations successively.

2.2.2 Random fields

A core element of GSTools is the spatial random field generator class SRF. A covariance model (Sect. 2.1) is needed to instantiate a spatial random field. We provide two ways for field generation: structured or unstructured. In both cases, the positions at which the field will be evaluated are given by a pos argument. In the structured case, pos contains one tuple per dimension, each defining the subdivision of the corresponding axis resulting in a rectilinear grid. For unstructured grids, the pos tuple contains the x, y and z coordinates of every evaluation point. SRF allows for controlling the underlying pseudo-random number generation by a seed to reproduce field generation. A code example is given in Fig. 11. Field generation is performed through the randomization method (Kraichnan, 1970; Heße et al., 2014), which utilizes the spectral density (Eq. 5) of the variogram model to approximate a Wiener process in Fourier space by

$$\begin{array}{}\text{(19)}& U\left(\mathit{x}\right)=\sqrt{{\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{2}}}{N}}}\cdot \sum _{i=\mathrm{1}}^N\left(Z_{\mathrm{1},i}\cdot \cos \left({\mathit{k}}_i\cdot \mathit{x}\right)+Z_{\mathrm{2},i}\cdot \sin \left({\mathit{k}}_i\cdot \mathit{x}\right)\right),\end{array}$$

where N is the number of Fourier modes of the approximation. The random variables $Z_{\mathrm{1},i},Z_{\mathrm{2},i}\sim \mathcal{N}(\mathrm{0},\mathrm{1})$ are mutually independent and are drawn from a standard normal distribution. The k_i values are mutually independent random samples, which are drawn from the spectral density with the aid of emcee, a Python package for Markov chain Monte Carlo sampling (Foreman-Mackey et al., 2013).

Figure 15 A workflow to generate a spatio-temporal random field with one spatial dimension.

The randomization method is implemented in the RandMeth class and used by default. The RandMeth routines create isotropic random fields. Thus, the corresponding covariance is radially symmetric, and the spectral density can be calculated by the Hankel transformation. Anisotropy is realized by rescaling and rotating the input points. The workflow allows users to generate a random field only from a given correlation, covariance or variogram function.

A key advantage of the randomization method implementation is the possibility to extend a generated SRF seamlessly, which not only preserves its statistical properties but also preserves the actual realization of the SRF. Potential applications are the following. (i) Particle simulations, where random incompressible velocity fields can be generated exactly at the location of the individual particles (see workflow in Sect. 2.3.1). It avoids interpolation errors, arising from grid-based velocity fields. (ii) If concentration plumes are simulated on a large domain, the SRF can be calculated on demand only for the time dependent spatial extent of the plume. And (iii) for high-performance computing applications, the field generation can be directly coupled to the domain decomposition, and each task only generates the SRF for its part of the domain. There are two main classes of alternative methods to the randomization method (Heße et al., 2014). By decomposing the covariance function, small spatial random fields can be computed very fast. But the computational cost quickly becomes unfeasible as the field grows in size. A second and quite popular class is the sequential Gaussian method, which can also create conditioned spatial random fields. However, numerical problems can arise if the underlying correlation function is smooth at the origin, and the computational costs also become unfeasible for highly resolved random fields (Emery, 2004).

Just like the kriging routines, the spatial random field generator allows for incorporating predefined trend, normalizer and mean for a greater variety of distributions. A special SRF class feature is the capability to perform variance upscaling to respect generation of random fields on mesh cells with a certain volume. We hereby use the upscaling method coarse graining (Attinger, 2003) to rescale the variance in Eq. (19) at each target point based on a given filter volume size λ:

$$\begin{array}{}\text{(20)}& {\mathit{\sigma }}^{\mathrm{2}}\left(\mathit{\lambda }\right)={\mathit{\sigma }}^{\mathrm{2}}\cdot {\left({\displaystyle \frac{{\mathrm{\ell }}^{\mathrm{2}}}{{\mathrm{\ell }}^{\mathrm{2}}+{\left(\frac{\mathit{\lambda }}{\mathrm{2}}\right)}^{\mathrm{2}}}}\right)}^{d/\mathrm{2}},\end{array}$$

where ℓ is the correlation length, and $\mathit{\lambda }=\sqrt[d]{V}$ is the filter size derived from the cell volume V depending on the field dimension, assuming the cell element is a hypercube. This approach was derived from the groundwater flow equation, assuming a Gaussian covariance model and should therefore be used with caution in differing scenarios. An example is provided in the workflow in Sect. 4.3.

2.2.3 Conditioned random fields

When point measurements of a target variable are available, they need to be considered when generating random fields. GSTools provides a class CondSRF combining kriging and random field generation, where we first derive the kriged field and the error variance and then generate a random field with zero mean where the variance in Eq. (19) is replaced with the estimated error variance. This procedure is advantageous to classical sequential Gaussian simulation (Webster and Oliver, 2007), as (i) we make use of the randomization method to generate a single random field, and (ii) we only need to solve the kriging problem once and not sequentially.

Figure 12 shows an example of an ensemble of conditioned random fields in one dimension. Where measurements of the target variables are available, all realizations satisfy them. However, random fields behave as unconditional fields (i.e. of an ensemble with identical parameters, like mean, variance and correlation length) where no point measurements are available (x>6). Characteristics, such as the ensemble variance, significantly change given the distribution of measurements and conditioning. The ensemble mean and the kriging field coincide, proving that the kriging field is the best linear unbiased predictor for the given data.

2.3 Additional features

2.3.1 Incompressible random vector field generation

Kraichnan (1970) was the first to suggest a randomization method for studying the diffusion of single particles in a random incompressible velocity field. He came up with a randomization method which includes a projector, ensuring the incompressibility of the vector field.

When $\stackrel{\mathrm{‾}}{U}$ is the mean velocity (oriented towards the first basis vector e₁), we generate random fluctuations with a given covariance model around $\stackrel{\mathrm{‾}}{U}$. And at the same time, we ensure that the velocity field remains incompressible; that is, $\mathrm{\nabla }\cdot \mathbf{U}=\mathrm{0}$ by using the randomization method (Eq. 19) and adding a projector p(k_i) to every mode being summed:

$$\begin{array}{}\text{(21)}& {\displaystyle }\begin{array}{rl}\mathbf{U}\left(\mathit{x}\right)& =\stackrel{\mathrm{‾}}{U}{\mathit{e}}_{\mathrm{1}}-\sqrt{{\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{2}}}{N}}}\sum _{i=\mathrm{1}}^{N}p\left({\mathit{k}}_i\right)\cdot \left(Z_{\mathrm{1},i}\cos \left({\mathit{k}}_i\cdot \mathit{x}\right)\right.\\ & \left.+Z_{\mathrm{2},i}\sin \left({\mathit{k}}_i\cdot \mathit{x}\right)\right),\end{array}\text{(22)}& {\displaystyle }p\left({\mathit{k}}_i\right)={\mathit{e}}_{\mathrm{1}}-{\displaystyle \frac{k_{i\mathrm{1}}}{|{\mathit{k}}_i{|}^{\mathrm{2}}}}\cdot {\mathit{k}}_i.\end{array}$$

Calculating $\mathrm{\nabla }\cdot \mathbf{U}=\mathrm{0}$ verifies that the resulting field is indeed incompressible. An example is shown in Fig. 13. Things like boundary conditions cannot be modelled with this method, but it can be used, for example, in transport simulations of the saturated subsurface (Schüler et al., 2016) or for studying turbulent open water (Kraichnan, 1970).

2.3.2 Field transformations

GSTools generates Gaussian random fields while real data often do not follow a Gaussian distribution. This is typically addressed through data transformation. GSTools provides a number of appropriate transformations beyond power transformations provided by the normalizer submodule (Sect. 2.1.4): (i) binary, (ii) discrete, (iii) boxcox (Box and Cox, 1964), (iv) zinnharvey (Zinn and Harvey, 2003), (v) normal_force_moments, (vi) normal_to_lognormal, (vii) normal_to_uniform, (viii) normal_to_arcsin and (ix) normal_to_uquad.

Figure 16 Generating spatial random fields on finite element method (FEM) meshes: either on cell centroids (middle) or mesh points (right). Plotting commands have been omitted.

Transformations can be combined sequentially to create more complex scenarios as in Fig. 14. Note that, in contrast to normalizers, transformations cannot be fitted to given data, which leaves the choice of the best transformation to the user.

2.3.3 Spatio-temporal modelling

Spatio-temporal modelling provides insights into time-dependent stochastic processes like rainfall, air temperature or crop yield, which is of high practical relevance. GSTools provides the metric spatio-temporal model (Cressie and Wikle, 2011) for all covariance models by enhancing the spatial with a time dimension, resulting in the spatio-temporal dimension d_st:

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{C}_m(\mathit{r},\mathrm{\Delta }t)& =C\left(\sqrt{\sum _{i=\mathrm{1}}^d{\left({\displaystyle \frac{r_i}{e_i}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{\Delta }t}{\mathit{\kappa }}}\right)}^{\mathrm{2}}}\right)\\ & =C\left(\sqrt{{\stackrel{\mathrm{̃}}{r}}^{\mathrm{2}}+\mathrm{\Delta }{\stackrel{\mathrm{̃}}{t}}^{\mathrm{2}}}\right),\end{array}\end{array}$$

where $\stackrel{\mathrm{̃}}{r}$ is the isotropified spatial distance as defined in Eq. (6), Δt is the temporal distance and $\mathrm{\Delta }\stackrel{\mathrm{̃}}{t}$ the isotropified temporal distance. The parameter κ accounts for a spatio-temporal anisotropy ratio and is the last entry of the anis array appended to the spatial anisotropy ratios. The implementation in GSTools enables the direct incorporation of spatial anisotropy and rotation in a spatio-temporal model. It further supports the use of arbitrary spatial dimensions in spatio-temporal models. Figure 15 shows the generation of a spatio-temporal random field with one spatial dimension.

2.3.4 Working on meshes

For improved handling of spatial random fields as input for PDE solvers like the finite element method (FEM), GSTools provides an interface for a number of mesh standards, such as meshio (Schlömer et al., 2021), PyVista (Sullivan and Kaszynski, 2019) and ogs5py (Müller et al., 2020). When using meshio or PyVista, the generated fields will be stored immediately in the mesh container. There are two options to generate a field on a given mesh: either on the points (points="points") or on the cell centroids (points="centroids"), which is important depending on the specification of the variable in the numerical scheme. Figure 16 provides an example.

Figure 17 Estimating the temperature variogram with geographic coordinates using the spherical Yadrenko model. Estimated length scale is ca. 0.9 (radians) and sill is around 13.

3 GSTools within the ecosystem of the GeoStat Framework

GSTools is part of a larger suite of Python packages, collectively hosted on GitHub under https://github.com/GeoStat-Framework (last access: 31 March 2022​​​​​​​​​​​​​​). The other packages in the GeoStat Framework complement some of the abilities of GSTools and form a comprehensive framework for geostatistical applications. We introduce some packages and demonstrate how they interact with, enhance and leverage the abilities of GSTools.

3.1  ogs5py

ogs5py (Müller et al., 2020) provides a Python API for the FEM-based OpenGeoSys 5 (Kolditz et al., 2012) scientific software suite for hydrogeological processes like groundwater flow and transport modelling where data scarcity is a typical shortcoming. Examples are point measurements of hydraulic head from observation wells and breakthrough curves from tracer experiments. Inferring hydraulic conductivity from these data requires a modelling framework with integrated stochastic data generation. The combination of GSTools with ogs5py enables a user to integrate the geostatistical modelling of an aquifer with hydrogeological simulations. Such an example for pumping test simulations is provided in Sect. 4.3.

3.2 welltestpy and AnaFlow

A pumping test is a cost-effective subsurface observation method typically used by hydrogeologists for aquifer characterization. The package welltestpy (Müller et al., 2021b) provides tools to handle, process, plot and analyse data from pumping test campaigns. It assists practitioners in identifying hydrogeological parameters by fitting measured drawdowns to some conceptual flow model. The package contains a number of examples that illustrate these abilities.

AnaFlow (Müller et al., 2021a) provides a wide range of analytical expressions for pumping tests under various conditions. Classical examples are Thiem's and Theis' solution assuming homogeneous aquifer properties. In addition, AnaFlow provides extended versions of both solutions, which account for aquifer heterogeneity and allow for estimating higher-order geostatistical parameters like variance and correlation length (Zech et al., 2012, 2016).

Figure 18 Setting up universal kriging with a drift function.

Figure 19 Workflow for regression kriging with a linear regression model.

3.3  PyKrige

GSTools provides an interface to the stand-alone package PyKrige (Murphy et al., 2021) for more specialized kriging applications. After 10 years of independent development, PyKrige has recently been migrated to the GeoStat Framework, and its functionality is currently integrated with the other packages. So far, the covariance model can be exchanged between the packages. In the future, PyKrige will become the kriging toolbox for the GeoStat Framework, providing advanced methods.

Figure 20 Plot of temperature measurements (a), universal kriging interpolation (b) and regression kriging results (c).

Figure 21 Scatterplot of latitude–temperature values (grey dots), the linear regression result (dashed orange line), universal kriging mean drift (dashed blue line) and the cross-sections of the respective kriging interpolation (solid lines).

3.4 Development, documentation and installation

GSTools is compatible with Python versions ≥3.6​​​​​​​, although previous releases support older versions of Python. Performance-critical parts, like variogram estimation (Sect. 2.1.4), kriging summation (Sect. 2.2.1) and the summation of the randomization method (Sect. 2.2.2), are implemented in Cython (Behnel et al., 2011). GSTools mainly depends on the SciPy ecosystem with its mandatory dependencies numpy (Harris et al., 2020) and scipy (Virtanen et al., 2020). The source code is maintained under a GitHub organization for optimizing team efforts. Users have the opportunity to communicate with developers by asking questions in a discussions forum, raising issues or improving code by making pull requests. All packages come with detailed documentation via https://readthedocs.org (last access: 31 March 2022​​​​​​​), which contains a range of tutorials explaining the features and a full API documentation created by Sphinx. Continuous integration is established through GitHub actions where Python wheels are pre-built for the most common operating systems (Windows, Linux, MacOS) and Python versions to enable simple installation. Each release on GitHub is directly deployed to the Python package index (https://www.pypi.org; last access: 31 March 2022​​​​​​​) as well as conda-forge (conda-forge community, 2015). An extensive set of unit tests are performed automatically and continuously through GitHub actions.

3.5 Interoperability

To integrate GSTools into the Scientific Python Stack, we provide a set of interfaces to other packages. These include the above-mentioned packages ogs5py, meshio, PyVista and pyevtk (https://github.com/pyscience-projects/pyevtk, last access: 31 March 2022​​​​​​​) for mesh operations. Other packages for geostatistics are also supported, such as PyKrige (Sect. 3.3) and scikit-gstat (Mälicke, 2022), the latter having a focus on variography and can be used for more detailed variogram estimation. For both packages, interfaces are provided to convert covariance models of GSTools to or from their counterparts in the respective package. Another package worth mentioning is verde (Uieda, 2018), a Python library for processing and gridding spatial data. Some of the features provided there can be easily combined with capabilities of GSTools such as detrending data to preprocess inputs.

4 Workflows

Having explained the core features of GSTools, we now provide a couple of example applications covering the topic of kriging, variogram estimation, random field generation and coupling with other tools to achieve more elaborate workflows. The examples illustrate the abilities of GSTools and serve as a starting point for a user's project development. All shown code snippets are taken from the actual workflow scripts and are not self-contained.

Figure 22 Initialization of an OGS model with mesh generation.

Figure 23 Generating correlated log-normal SRFs adapted to the mesh settings of the numerical model for three connectivity structures following the Zinn and Harvey (2003) transformation.

4.1 Regression kriging vs. universal kriging: finding a north–south temperature trend

Kriging is a well-established interpolation method applied in many fields of natural science. We compare two options of incorporating auxiliary variables to calculate the kriging weights: (i) regression kriging (RK), where the trend of input data is estimated by regression and simple kriging is applied to the residuals, and (ii) universal kriging (UK), where the trend model is used as the internal drift in the kriging system. The methods differ in use of the covariance model. The linear RK does not incorporate spatial correlation information, while UK does through the drift function for calculating the kriging weights. Both methods are often considered to provide mathematically equal results, but we show that there are sensitive differences. The resources for this workflow are provided in Müller (2021).

As a data basis, we use measured temperature of the German weather service retrieved with the Python package wetterdienst (Gutzmann et al., 2021), which we examine for a linear north–south trend. We use the established spherical covariance model in its Yadrenko variant suitable for geographical coordinates. Variogram estimation and fitting results are shown in Fig. 17.

Figure 24 Tracer transport simulation results: spatially distributed concentration plumes after 15 d with transmissivity distributions in the background.

Figure 25 Workflow to generate an ensemble of transmissivity fields on a given mesh (a). A single realization in shown in the right plot (b).

Figure 26 Comparison the ensemble mean drawdown 〈h(r,t)〉 (a, d) with the effective head solution h_CG(r,t) (c, f) for two parameter sets. The vanishing absolute difference between both (b, e) shows that they perfectly agree.​​​​​​​

Figure 18 shows how to set up the UK estimator, including the drift function, and Fig. 19 shows the setup of the RK estimator. RK requires the preceding step of fitting the regression model for the trend of the Detrended kriging routine. The interpolation results are shown in Fig. 20, indicating that both methods provide equally good results.

Figure 27 Transmissivity of the Herten aquifer analogue and locations of virtual observations marked as black dots.

Figure 28 Variogram estimation and resulting experimental (dots) and fitted variogram γ (line) of the Herten aquifer analogue.

Figure 21 shows the estimated mean trends for both UK and RK, revealing completely contrary results. The RK result indicates an increase in mean temperature with increasing latitude, which seems reasonable given a raising terrain elevation from the Baltic Sea in the north towards the Alps in the south. The estimated mean of UK shows the opposite with temperature decreasing with latitude. A potential explanation here is the general temperature increase towards the Equator. While the UK mean fits better with the cross-section at 10^∘ longitude (Fig. 21), the RK mean fits the scatter diagram better, as expected.

4.2 Heterogeneous transport simulation: the impact of connectivity

The combination of ogs5py and GSTools makes it possible to quickly set up and run subsurface flow and transport simulations in a heterogeneous aquifer setting. The critical step is the generation of a spatially distributed hydraulic conductivity distribution, adapted to the numerical simulation grid. We further demonstrate GSTools' ability to generate different connectivity structures, and we discuss their impact on transport results. The resources for this workflow are provided in Müller and Zech (2021a).

A flow and transport model is initialized through an instance of the OGS class from ogs5py, with simple mesh generation (Fig. 22) and specification of model parameters and boundary conditions. Random fields are initialized through the SRF class (Fig. 23). By passing the subclass model.msh, mesh details are transferred for generating distributed values at particular mesh locations, even for unstructured grids. The subroutine transform.zinnharvey allows for generating Gaussian structures where the mean values of the field are not connected but the low- or high-conductivity areas are connected, using the transformation by Zinn and Harvey (2003). Note that the correlation lengths undergo rescaling (Gong et al., 2013). The concept of connectivity follows the paper by Zinn and Harvey (2003), where connectedness refers to connected paths of extreme or special values in the conductivity field.

Figure 29 One realization of the conditioned SRF (a) and absolute difference (b) between the “true” (Fig. 27) and conditioned transmissivity, showing increasing differences with distance from the conditioning area (rectangle).

Figure 30 Generation of an ensemble of 20 conditional realizations and transmissivity transects T(x) at y=4 m. The thick blue line is the true transmissivity (Fig. 27), and the shaded area shows the range of 1 standard deviation calculated from 20 realizations of conditioned fields. Black points indicate the observations.

Simulated tracer plumes in Fig. 24 show the particular effects of connectivity: the plume remains relatively compact for classical Gaussian fields, where mean values are connected. Transformed fields lead to more disrupted, dynamic plumes, which is mostly caused by trapping in areas of connected low conductivity and preferential flow in connected high-conductivity areas.

4.3 Characterizing mean drawdowns of a pumping test ensemble

Combining flow simulations in ogs5py with random fields of GSTools allows for performing Monte Carlo studies to identify ensemble mean behaviour. Zech et al. (2016) made use of this workflow to prove the applicability of an effective drawdown solution for pumping tests in random conductivity. We present a short form of their workflow, which is accessible in Müller and Zech (2021b).

The flow model is initialized through the OGS class with model parameters and boundary conditions creating the convergent flow setting of a pumping test. The mesh generation and time stepping can be specifically adapted to the non-uniform flow conditions. Ensembles of heterogeneous transmissivity fields are generated with the SRF class where reproducibility is controlled by the seed and where normal fields are converted in place with normalizer.LogNormal as shown in Fig. 25.

The implementation of the randomization method (Sect. 2.2.2) allows for the adaption of random fields to the non-uniform grid. The associated variance upscaling follows the coarse graining procedure for Gaussian variograms according to Eq. (20).

Calculated ensemble means can be compared to analytical solutions (Fig. 26), such as Theis' solution for homogeneous media or the effective drawdown solution by Zech et al. (2016), making use of their implementations in welltestpy and AnaFlow.

4.4 Geostatistical exercises with the Herten aquifer

We demonstrate how to estimate variograms and how to condition spatial random fields on observations using data from the Herten aquifer analogue (Bayer et al., 2011). The aquifer analogue was created from surveying multiple outcrop faces at a gravel pit, situated in the Rhine valley in southern Germany. The 2D information was interpolated to a 3D data set, including hydraulic, thermal and chemical information. The workflow files are provided in Schüler and Müller (2021).

We determine spatial correlations through variogram estimation using gstools.variogram. First, we identify the full transmissivity structure. The aquifer analogue data are given in a facies structure with one hydraulic conductivity value K per facies. We calculate transmissivity by integrating the hydraulic conductivity over the vertical axis $T(x,y)=\int K(x,y,z)\mathrm{d}z$. The structured transmissivity is shown in Fig. 27, which we consider as true distribution for the following exercises.

We select 13×13 virtual observations on a rectangular grid (Fig. 27), covering a subarea of about 42 m². These observations are used to determine the empirical variogram shown in Fig. 28. We fit an exponential covariance model to the data, which suits well with a coefficient of determination of R²=0.913.

We use the fitted exponential variogram model and ordinary kriging to create conditioned spatial random fields with CondSRF. Figure 29 shows one realization and the absolute difference to the true transmissivity (Fig. 27). Differences grow with increasing distance from observations. This trend can be even better seen in a transmissivity transect shown in Fig. 30. The standard deviation calculated from 20 realizations of conditioned SRFs shows that deviations from the reference field are significantly lower close to observation points.

5 Conclusions

The GSTools package provides a Python-based platform for geostatistical applications. It is similar to software packages like gstat for R or stand-alone packages like TPROGS (Carle, 1999), GSLIB (Deutsch and Journel, 1997) and S-GeMS (Remy, 2005). However, we believe that a comprehensive and ready-made geostatistical software package for Python has advantages, simply through the choice of the programming language, as it has a gentle learning curve, is often used as a glue language and is widely adopted by the scientific community. Salient features of GSTools are its random field generation and its versatile covariance model. It is furthermore integrated with other Python packages, like PyKrige (Murphy et al., 2021), ogs5py (Müller et al., 2020) or scikit-gstat (Mälicke, 2022), and provides interfaces to meshio (Schlömer et al., 2021) and PyVista (Sullivan and Kaszynski, 2019). GeoStat Examples (https://github.com/GeoStat-Examples, last access: 31 March 2022​​​​​​​) provides a number of applications, including the four presented workflows. They showcase the abilities of GSTools and can serve as a starting point for practitioners to develop their own solutions for the geostatistical problems they face.