Theoretical Variogram Models: Spherical, Exponential & Gaussian

An empirical variogram is a scatter of noisy points; kriging needs a smooth, mathematically valid function evaluated at any distance. Fitting a theoretical variogram model bridges that gap by replacing the raw lag estimates with a continuous curve whose parameters, nugget, sill, and range, encode the spatial continuity of your variable. Choosing the right family and fitting it correctly is the single most consequential modelling decision in geostatistics, because those parameters flow directly into the interpolation weights and prediction variance. This page sits within Variogram Modeling & Semivariance Analysis and covers the admissible model families, their equations, and a fully annotated scikit-gstat fitting workflow.

Prerequisites #

• Python 3.10+
• numpy>=1.24, scipy>=1.10, scikit-gstat>=1.0, geopandas>=0.14, matplotlib>=3.7
• A completed empirical variogram estimation with sensible lag bins
• Sample coordinates in a projected (metric) CRS so that ranges are expressed in metres
• At least 30 to 50 well-distributed samples; more if you expect a nugget or long range

Mathematical Core #

An empirical variogram estimates the semivariance $\hat{\gamma}(h)$ at a set of lag distances $h$. A theoretical model $\gamma(h; \boldsymbol{\theta})$ is a parametric function chosen so that the fitted curve is admissible, meaning it produces a valid covariance structure for kriging. Every model is built from three interpretable parameters:

• Nugget $c_0$: the apparent semivariance at $h \to 0^+$, capturing measurement error and micro-scale variation below the shortest sample spacing.
• Partial sill $c$: the additional structured variance contributed by spatial correlation, so the total sill is $c_0 + c$.
• Range $a$: the distance over which spatial correlation is expressed; beyond it, pairs are effectively independent.

The admissibility constraint #

A model is only usable if it is conditionally negative-definite. Formally, for any set of $n$ locations and any real weights $\lambda_i$ satisfying $\sum_i \lambda_i = 0$,

$$\sum_{i=1}^{n} \sum_{j=1}^{n} \lambda_i \lambda_j \, \gamma(\|\mathbf{x}_i - \mathbf{x}_j\|) \le 0 .$$

Equivalently, the associated covariance function $C(h) = (c_0 + c) - \gamma(h)$ must be positive-definite. This is why you cannot fit an arbitrary polynomial or spline to the empirical points: a curve that fits the dots beautifully but violates this condition can produce a kriging system with a negative prediction variance, which is physically meaningless. Only a small library of proven families is admissible, and the ones below are the workhorses.

Spherical model #

$$\gamma(h) = c_0 + c\left[\frac{3h}{2a} - \frac{1}{2}\left(\frac{h}{a}\right)^{3}\right] \quad \text{for } 0 < h \le a, \qquad \gamma(h) = c_0 + c \ \text{ for } h > a .$$

The spherical model rises almost linearly near the origin, then bends over and reaches the sill exactly at the range $a$. Its clean, finite range makes it the pragmatic default for many earth-science variables (ore grade, soil properties, rainfall).

Exponential model #

$$\gamma(h) = c_0 + c\left[1 - \exp\!\left(-\frac{h}{a}\right)\right] .$$

The exponential model approaches the sill asymptotically and never reaches it. It rises linearly from the origin, indicating a continuous but non-smooth field. Its practical, or effective range is the distance at which it attains 95 percent of the sill, approximately $3a$.

Gaussian model #

$$\gamma(h) = c_0 + c\left[1 - \exp\!\left(-\frac{h^{2}}{a^{2}}\right)\right] .$$

The Gaussian model rises parabolically (flat) near the origin, implying an extremely smooth, mean-square-differentiable phenomenon such as a potentiometric surface. Its effective range is roughly $\sqrt{3}\,a$. A pure Gaussian model with no nugget produces an ill-conditioned kriging matrix, so always include a small nugget in practice.

Matern model #

$$\gamma(h) = c_0 + c\left[1 - \frac{1}{2^{\nu-1}\Gamma(\nu)}\left(\frac{h}{a}\right)^{\nu} K_\nu\!\left(\frac{h}{a}\right)\right],$$

where $K_\nu$ is the modified Bessel function of the second kind and $\Gamma$ is the gamma function. The smoothness parameter $\nu$ generalises the family: $\nu = 0.5$ reproduces the exponential model, and $\nu \to \infty$ approaches the Gaussian. Fitting $\nu$ lets the data choose its own near-origin smoothness rather than forcing an exponential-or-Gaussian dichotomy.

Comparing the shapes #

The near-origin behaviour is the diagnostic that separates the families: Gaussian is flat (parabolic), exponential is linear, and spherical is linear but reaches a hard sill. Match that shape to the physical smoothness of your variable, then let the fit confirm it.

Annotated Implementation #

scikit-gstat bundles empirical estimation and theoretical fitting into a single Variogram object. Selecting model='spherical' (or 'exponential', 'gaussian', 'matern') triggers a least-squares fit; the estimated parameters are then available on .parameters.

import numpy as np
import geopandas as gpd
from skgstat import Variogram

# ── 1. Load projected samples (metres) ──────────────────────────────────────
gdf = gpd.read_file("soil_samples.gpkg").to_crs("EPSG:32633")  # UTM 33N
coords = np.column_stack([gdf.geometry.x, gdf.geometry.y])      # shape (n, 2)
values = gdf["zinc_ppm"].to_numpy(dtype=float)                  # target variable

# ── 2. Fit a spherical model to the empirical variogram ─────────────────────
V = Variogram(
    coordinates=coords,
    values=values,
    model="spherical",     # theoretical family
    estimator="matheron",  # classical semivariance estimator
    bin_func="even",       # equal-width lag bins
    n_lags=15,             # number of lag classes
    maxlag="median",       # cap lags at the median pairwise distance
    fit_method="trf",      # trust-region least squares
    fit_sigma="linear",    # down-weight distant, noisy lags
    normalize=False,
)

# ── 3. Read the fitted parameters ───────────────────────────────────────────
# .parameters returns [range, sill, nugget] for the spherical/exp/gauss models
rng, sill, nugget = V.parameters
print(f"nugget c0 = {nugget:.3f}")
print(f"sill      = {sill:.3f}   (partial sill = {sill - nugget:.3f})")
print(f"range a   = {rng:.1f} m")
print(f"fit RMSE  = {V.rmse:.4f}")

A few contract details worth pinning down. The Variogram constructor performs both the empirical binning and the model fit on instantiation; there is no separate .fit() call required, but you can trigger a refit after changing an attribute:

# Reconfigure and refit in place
V.model = "exponential"   # switching the family triggers an automatic refit
V.n_lags = 20             # changing binning also refits
params = V.parameters     # [effective_range, sill, nugget] for exponential

# Evaluate the fitted model at arbitrary distances (for kriging or plotting)
h = np.linspace(0, V.maxlag, 200)
gamma_h = V.fitted_model(h)     # callable returning theoretical semivariance

Note that for the exponential and Gaussian models scikit-gstat reports the effective range in .parameters[0], not the raw scale parameter $a$ from the equations above. This is a deliberate convenience so that the reported range is comparable across families; keep it in mind when you translate parameters into another library’s kriging call.

Fitting and comparing every admissible family #

Rather than eyeballing one model, fit the full admissible set and rank them. Model selection should penalise complexity so that the extra parameter of the Matern family has to earn its place.

import numpy as np
from skgstat import Variogram

def fit_and_score(coords, values, families=("spherical", "exponential",
                                            "gaussian", "matern")):
    """Fit each admissible model and rank by RMSE and AIC.

    Returns a list of dicts sorted by AIC (lower is better).
    """
    results = []
    for fam in families:
        V = Variogram(coords, values, model=fam, n_lags=15,
                      maxlag="median", fit_method="trf", fit_sigma="linear")

        # Residual sum of squares between empirical and fitted semivariance
        emp = V.experimental                      # empirical semivariances
        bins = V.bins                             # lag bin centres
        fit = V.fitted_model(bins)
        rss = np.nansum((emp - fit) ** 2)
        m = np.sum(~np.isnan(emp))                # number of usable lags
        k = len(V.parameters)                     # free parameters (+1 for matern)

        # Small-sample-corrected AIC for least-squares fits
        aic = m * np.log(rss / m) + 2 * k
        results.append({"model": fam, "rmse": V.rmse, "aic": aic,
                        "params": V.parameters})

    return sorted(results, key=lambda r: r["aic"])

ranking = fit_and_score(coords, values)
for r in ranking:
    print(f"{r['model']:<12} RMSE={r['rmse']:.4f}  AIC={r['aic']:.2f}")

best = ranking[0]
print(f"\nSelected model: {best['model']}  ->  {best['params']}")

The AIC above uses the least-squares form $\text{AIC} = m \ln(\text{RSS}/m) + 2k$, where $m$ is the number of populated lag bins and $k$ the number of fitted parameters. Because the Matern model carries an extra smoothness parameter, its RSS must improve enough to overcome the $+2$ penalty before it is preferred. When two models tie on AIC, prefer the one whose near-origin behaviour matches the physics of your variable.

Diagnostic and Parameter Configuration #

The empirical binning drives everything downstream, so the lag configuration deserves as much care as the model choice; the details are covered in empirical variogram estimation. Within the fitting step, three knobs matter most:

Weighting scheme (fit_sigma). Ordinary least squares (fit_sigma=None) treats every lag equally, which lets sparse, noisy far lags distort the near-origin fit that kriging depends on. Set fit_sigma='linear' or 'exp' to weight lags by reliability, or supply per-lag pair counts explicitly. Weighted least squares that weights each lag by $N(h)/\gamma(h)^2$ is the classical Cressie criterion and is a robust default.

# Explicit pair-count weighting (Cressie-style emphasis on near lags)
V = Variogram(coords, values, model="spherical", n_lags=15,
              fit_method="trf", fit_sigma="linear")

# Inspect how many pairs populate each lag — thin far lags are unreliable
for h, n in zip(V.bins, V.bin_count):
    flag = "  <-- sparse" if n < 30 else ""
    print(f"lag {h:8.1f} m : {n:4d} pairs{flag}")

Nugget handling. By default scikit-gstat fits the nugget as a free parameter. If domain knowledge says the process is continuous with negligible measurement error, force a zero nugget with use_nugget=False to avoid absorbing structured variance into $c_0$. Conversely, a fitted nugget approaching the full sill signals a pure nugget effect: no spatial structure was resolved, and kriging will collapse toward the global mean.

Maximum lag. Fitting to lags beyond half the study extent uses pairs supported by too few observations and inflates the apparent range. maxlag='median' (median pairwise distance) or a fixed fraction of the diagonal keeps the fit anchored to well-estimated lags.

Output Interpretation #

Read the fitted parameters against the empirical scatter and the physics of the variable:

• Nugget-to-sill ratio. $c_0 / (c_0 + c) < 0.25$ indicates strong spatial structure; $0.25$ to $0.75$ moderate; $> 0.75$ weak structure dominated by noise or micro-scale variation. This ratio, not the absolute nugget, tells you how much kriging can improve on the global mean.
• Range in context. A range far shorter than your sample spacing means most prediction locations fall beyond the correlation distance, so kriging degrades to the mean; a range spanning the whole study area suggests an unmodelled trend that should be handled with universal kriging.
• Effective vs nominal range. For exponential and Gaussian models, remember the reported effective range is the 95-percent-of-sill distance. Do not confuse it with the scale parameter when hand-computing covariances.
• RMSE and AIC together. A low RMSE with a physically implausible shape (for example a Gaussian fit to obviously rough data) is a warning sign, not a success. Cross-check the winning model’s near-origin behaviour.

Good signs: a monotone rise to a stable sill, a nugget below a quarter of the sill, and a range comfortably larger than the mean sample spacing but smaller than half the domain. Warning signs: a variogram that keeps climbing without levelling (trend or non-stationarity), a sill that dips at long lags (a hole effect or periodicity), or a fitted range that pins to the maxlag boundary.

Production Considerations #

Computational cost. Empirical variogram estimation is $O(n^2)$ in the number of pairs; the theoretical fit itself is cheap (a handful of parameters via scipy.optimize). For $n > 10{,}000$ points, subsample or estimate the variogram on a spatially stratified sample rather than all $\binom{n}{2}$ pairs. scikit-gstat supports maxlag to prune the far, expensive lags that add little information.

Reproducibility. Pin scikit-gstat, numpy, and scipy versions. The fitted parameters depend on the optimiser (fit_method='trf' vs 'lm'), the lag binning, and the weighting scheme; record all three alongside the CRS EPSG code so a stored model is fully reproducible.

Serialisation for kriging. Persist the three parameters (nugget, sill, range) plus the family name and CRS as a small JSON or dictionary. Downstream kriging engines (pykrige, gstools) accept these directly; do not re-fit on every pipeline run.

import json

model_card = {
    "family": V.model.__name__ if callable(V.model) else str(V.model),
    "nugget": float(V.parameters[2]),
    "sill": float(V.parameters[1]),
    "range_m": float(V.parameters[0]),
    "crs": "EPSG:32633",
    "estimator": "matheron",
}
with open("variogram_model.json", "w") as fh:
    json.dump(model_card, fh, indent=2)

Numerical stability. A pure Gaussian model with a zero nugget produces a nearly singular kriging covariance matrix and unstable weights. Always retain a small nugget (even 1 to 3 percent of the sill) for Gaussian and Matern fits with high smoothness to keep the linear system well-conditioned.

Troubleshooting #

Next Steps #

With an admissible model chosen, drill into the mechanics of fitting spherical, exponential and Gaussian variogram models and the fine detail of estimating nugget, sill and range parameters. If your variogram behaves differently along different directions, the fitted range becomes directional, which is the subject of anisotropy and directional variograms. The finished model then feeds the interpolation weights in the kriging, interpolation and surface generation techniques workflow.

Related

• Fitting Spherical, Exponential & Gaussian Variogram Models — step-by-step least-squares fits for each family
• Estimating Nugget, Sill & Range Parameters — interpreting and constraining the three core parameters
• Empirical Variogram Estimation — binning sample pairs into the empirical curve you fit here
• Anisotropy & Directional Variograms — direction-dependent ranges for anisotropic fields
• Kriging, Interpolation & Surface Generation Techniques — where the fitted model becomes prediction weights

← Back to Variogram Modeling & Semivariance Analysis