`variogram'

Variogram models are coded as the sum of one or more simple models (and optionally an anisotropy structure). A simple variogram model is denoted by

$c \, Mod(a)$

with $c$ the vertical (variance) scaling factor, $Mod$ the model type, and $a$ the range (horizontal, distance scaling factor) of this simple model. If $Mod$ is a transitive model (i.e. after some distance, or asymptotically as $h \rightarrow \infty$ the simple model reaches a certain maximum) then $c$ is the partial sill of that model and the simple covariogram that corresponds to this variogram model is:

$c(1 - Mod(a))$

Figure 4.1: Example variogram 8 Nug() + 12 Sph(10). Semivariance representation (a) and covariance representation (b)

(a) (b)

Fig. 4.1 gives an example of the variogram model 8 Nug() + 12 Sph(10) as semivariance and covariance representation. Unit models available in gstat are listed in table 4.1.

Table 4.1: Simple variogram models in gstat: the building blocks for a variogram model

model	syntax	$\gamma(h)$	$h$ range
Nugget	1 Nug(0)	0	$h = 0$
		1	$h > 0$
Spherical	1 Sph(a)	$\frac{3h}{2a}-\frac{1}{2}(\frac{h}{a})^3$	$0 \le h \le a$
		$1$	$h > a$
Exponential	1 Exp(a)	$1 - \exp(\frac{-h}{a})$	$h \ge 0$
Linear	1 Lin(0)	$h$	$h \ge 0$
Linear-with-sill4.2	1 Lin(a)	$\frac{h}{a}$	$0 \le h \le a$
		$1$	$h > a$
Circular	1 Cir(a)	$\frac{2 h}{\pi a}\sqrt{1 - (\frac{h}{a})^2} + \frac{2}{\pi} \arcsin\frac{h}{a}$	$0 \le h \le a$
		$1$	$h > a$
Pentaspherical	1 Pen(a)	$\frac{15h}{8a}- \frac{5}{4}(\frac{h}{a})^3 +\frac{3}{8}(\frac{h}{a})^5$	$0 \le h \le a$
		$1$	$h > a$
Gaussian	1 Gau(a)	$\gamma(h) = 1 - \exp(-(\frac{h}{a})^2)$	$h \ge 0$
Bessel4.3	1 Bes(a)	$1 - \frac{h}{a}{K_1}(\frac{h}{a})$	$h \ge 0$
Logarithmic	1 Log(a)	$0$	$h = 0$
		$\log(h + a)$	$h > 0$
Power	1 Pow(a)	$h ^ a$	$h \ge 0, 0 < a \le 2$
Periodic	1 Per(a)	$1 - \cos(\frac{2\pi h}{a})$	$h \ge 0$

Note that the Exp(), Gau() and Bes() models reach their sill asymptotically (as $h \rightarrow \infty$). The `effective range', is the distance where the variogram reaches 95% of its maximum, and this is 3 a for Exp(a), $\sqrt{3}$ a for Gau(a) and 4 a for Bes(a). The logarithmic and power model are unbounded (and are therefore not suitable for covariance modelling or simple kriging).

Pseudo cross variograms may have a non-zero value for $h = 0$. An intercept can be defined as a constant added to the variogram model, e.g.

1.5 + 0.5 Nug() + 2.2 Sph(20)

or, equivalently

1.5 Int() + 0.5 Nug() + 2.2 Sph(20)

and can be fitted only when a sample variogram estimate is available at zero distance.

Geometric anisotropy

Geometric anisotropy can be modelled for each individual simple model by addition of two or five anisotropy parameters after the range, e.g.

$c \, Mod(a,p,s)$

for 2-d anisotropy, or, for 3-d anisotropy:

$c \, Mod(a,p,q,r,s,t)$

Using anisotropy, the variogram model range parameter ($a$) is the maximum range, that of the major direction of continuity (direction of spatial correlation at longest distances). The range in the direction perpendicular to the major direction is the minor range. The anisotropy ratio is the ratio between the minor range and the major range (a value between 0 and 1).

A two-dimensional range ellipse is defined by $(a,p,s)$, a three-dimensional ellipsoid is defined by $(a,p,q,r,s,t)$. In the two-dimensional case $p$ is the angle for the principal direction of continuity (measured in degrees, clockwise from positive Y, north), and $s$ is the anisotropy ratio. So, the range in the major direction ($p$) is $a$, and the range in the minor direction ($p+90$) is $as$.

Figure 4.3: anisotropy ellipse

In three dimensions $p$ is the angle for the principal direction of continuity (measured in degrees, clockwise from Y), $q$ is the dip angle for the principal direction of continuity (measured in positive degrees up from horizontal), $r$ is the third rotation angle to rotate the two minor directions around the principal direction defined by $p$ and $q$. A positive angle acts clockwise while looking in the principal direction. Anisotropy ratios $s$ and $t$ are the ratios between the major range and each of the two minor ranges. (Note that $c \, Mod(a,p,s)$ is equivalent to $c \, Mod(a,p,0,0,s,1)$.)

Zonal anisotropy

Zonal anisotropy (anisotropy in the sill) can be obtained by defining geometric anisotropy with large anisotropy ratios. For instance, when the spatial working domain is not too large (the largest distance in the area considered does not exceed, let's say, 100), then the model

1 Sph(2e5, 90, 1e-5)

will have nearly zero values (and thus be practically absent) in direction 90 (east-west), whereas it will reach the sill (1) at distance 2 in direction 0 (north-south).