Additivity: Definition, Semivariogram Models, and Geostatistics

In petroleum geostatistics, additivity is the mathematical property that allows two or more valid (admissible) semivariogram models to be summed with positive coefficients to produce a new model that is itself valid. Because the resulting nested model remains positive semi-definite, it can be used directly in kriging and stochastic simulation without introducing non-physical negative estimation variances. Additivity is the theoretical foundation for fitting complex, multi-scale spatial variability in reservoir characterization models, where a single basic model is rarely sufficient to capture all of the structural features seen in the experimental semivariogram computed from well data.

Key Takeaways

• Any linear combination of admissible semivariogram models with strictly positive coefficients is itself admissible, preserving the positive semi-definiteness required for kriging.
• The experimental semivariogram gamma(h) = 0.5 * E[(Z(x+h) - Z(x))^2] quantifies how spatial correlation decays with separation distance h; the fitted model must honour additivity to guarantee valid kriging variances.
• Nested models combine a nugget component (representing micro-scale variability or measurement error) with one or more structured components, each with its own sill and range.
• The four basic admissible models most commonly used in petroleum applications are the spherical, exponential, Gaussian, and power models; any positive-weight sum of these is also admissible.
• Additivity extends to multivariate settings (co-regionalization) and to anisotropic models, enabling reservoir geologists and petrophysicists to represent directional variability in permeability and porosity simultaneously.

Definition and Theoretical Background

The semivariogram (also called the variogram) is the central tool of geostatistics. For a spatial random function Z(x), the semivariogram at lag h is defined as:

gamma(h) = 0.5 * E[ (Z(x + h) - Z(x))^2 ]

where E denotes statistical expectation and x is a location vector in two or three dimensions. In practice, gamma(h) is estimated from a finite set of data pairs separated by distance h (within some tolerance band) and then fitted with a mathematical model. The fitted model must be conditionally negative semi-definite, which is equivalent to its associated covariance function C(h) = C(0) - gamma(h) being positive semi-definite. This condition, rooted in Bochner's theorem of harmonic analysis, ensures that the covariance matrix assembled during kriging is always positive definite, making the kriging system solvable and the resulting estimation variance non-negative.

The additivity principle states that if gamma_1(h) and gamma_2(h) are both admissible semivariogram models, then:

gamma(h) = w_1 * gamma_1(h) + w_2 * gamma_2(h), with w_1 > 0 and w_2 > 0

is also admissible. This can be extended to any finite number of components. The proof follows directly from the fact that a positive linear combination of positive semi-definite functions is positive semi-definite. In physical terms, each component in a nested model represents a distinct spatial scale of variability, and the total variability at any lag is the sum of contributions from all scales.

How It Works

When a geostatistician computes an experimental semivariogram from wireline log data or core measurements across a set of wells, the resulting scatter of gamma estimates rarely conforms to a single basic model. A typical reservoir dataset might show a rapid initial rise of gamma with lag (indicating short-range heterogeneity from diagenetic pore-fill or thin-bed lamination), followed by a more gradual approach to a second, higher sill (reflecting long-range depositional architecture such as facies belts or sequence boundaries). Fitting a single spherical model to this pattern would either overestimate short-range variability or underestimate the total variance, producing a biased interpolation.

The nested model solution uses the additivity principle to decompose the experimental semivariogram into a nugget component plus two or more structured components. A practical example for a fluvial sandstone reservoir might be:

• Nugget component: gamma_0 = C_0 * nugget(h), with sill C_0 = 0.05 representing measurement noise and sub-core-scale variability in porosity.
• Short-range spherical component: gamma_1 = C_1 * Sph(h; a_1), with sill C_1 = 0.25 and range a_1 = 50 m (164 ft), capturing diagenetic heterogeneity within individual flow units.
• Long-range spherical component: gamma_2 = C_2 * Sph(h; a_2), with sill C_2 = 0.70 and range a_2 = 800 m (2,625 ft), capturing facies-belt-scale variability tied to sequence stratigraphy.

The total sill of the nested model equals C_0 + C_1 + C_2 = 1.00, which equals the total variance of the porosity dataset (after standardisation). Each component weight is positive, satisfying the additivity condition. When this model is fed into ordinary kriging or sequential Gaussian simulation (SGS), the algorithm correctly honours both the local heterogeneity at scales below 50 m and the regional trend at scales up to 800 m. Without the nested structure, either the small-scale or the large-scale variability would be smeared or lost, degrading the quality of the reservoir model.

The fitting procedure itself is typically iterative: the modeller makes an initial visual estimate of the number of structures, their sills, and their ranges from the experimental semivariogram, then uses weighted least squares or maximum likelihood to refine the parameters. Software packages such as GSLIB, Petrel (Schlumberger/SLB), and RMS (Roxar/Emerson) all implement nested model fitting and enforce the additivity constraint automatically by restricting all component weights to positive values. In Petrel, the variogram modelling panel displays each nested component as a colour-coded curve, with the sum shown as a bold line against the experimental points, giving the modeller immediate visual feedback on the quality of fit.

Basic Admissible Models Used in Nested Structures

Four mathematical models dominate petroleum geostatistics practice. Each is admissible on its own and can therefore participate in any nested combination.

• Spherical model: Rises linearly near the origin, reaches the sill exactly at the range a, and stays flat beyond. It is the most commonly used model in reservoir modelling because its finite range aligns well with the concept of a correlation length tied to a depositional body. The mathematical form is Sph(h; a) = 1.5*(h/a) - 0.5*(h/a)^3 for h <= a, and = 1 for h > a.
• Exponential model: Rises more steeply near the origin than the spherical model and approaches the sill asymptotically. The practical range (distance at which 95% of the sill is reached) is approximately 3a. It is favoured for highly heterogeneous media such as fractured carbonates or tight gas sands where short-range variability dominates.
• Gaussian model: Has a parabolic behaviour near the origin (indicating a very smooth, continuously differentiable spatial variable) and an inflection point near 0.58a before approaching the sill asymptotically. It is used for laterally continuous, gradationally varying properties such as acoustic impedance or net-to-gross in deltaic systems.
• Power model: Has no sill and models non-stationary variability where gamma(h) grows indefinitely as a power of h. It is used sparingly, typically for regional-scale trend components in basin-wide models, and must always be paired with bounded components in a nested structure to prevent infinite kriging variances at large lags.
• Nugget model: A pure discontinuity at the origin: gamma(0) = 0 and gamma(h) = 1 for all h > 0. It represents variability at scales smaller than the smallest lag sampled (micro-scale heterogeneity, measurement error, or both). A pure nugget effect means there is no spatial correlation at any sampled scale.

Anisotropy and Directional Nested Models

Most reservoir properties are anisotropic: they vary more rapidly in one direction than another. Horizontal permeability in a braided fluvial system, for instance, may be correlated over hundreds of metres along the palaeocurrent direction but only tens of metres perpendicular to it. Vertical permeability is typically correlated over only a few metres due to layering. Additivity handles anisotropy by assigning a separate geometric anisotropy ratio (the ratio of the range along the major axis to the range along the minor axis) to each nested component.

In the 3D case, each component gamma_i(h) is evaluated using a transformed lag distance that accounts for both the azimuthal anisotropy in the horizontal plane and the vertical-to-horizontal range ratio. The total nested model gamma(h) = sum_i [w_i * gamma_i(h)] inherits the positive semi-definiteness of its components regardless of the anisotropy parameters, because anisotropy is introduced through a linear transformation of the coordinate system and linear transformations preserve positive semi-definiteness. This means a geostatistician can assign different anisotropy ratios and orientations to different nested components, for example one component with a north-northeast major axis for a depositional trend and a second component with a roughly isotropic geometry for diagenetic overprinting, and the additivity guarantee still holds.

In practice, estimating anisotropy parameters requires computing directional experimental semivariograms in multiple azimuths and at different dip angles. In a sparse well dataset (5 to 20 wells covering hundreds of square kilometres), horizontal anisotropy is often constrained by analogue outcrop data, seismic attribute maps, or depositional process models rather than from the well data alone. The vertical semivariogram, however, can usually be computed reliably from high-resolution wireline logs with sample intervals of 0.15 m (0.5 ft) or finer.

Co-Regionalization: Multivariate Additivity

When two or more reservoir properties (for example, porosity and acoustic impedance derived from seismic inversion) are modelled jointly, the cross-semivariogram between variables must also be admissible. The linear model of co-regionalization (LMC) extends additivity to the multivariate case: the cross-semivariogram between variables u and v at scale i is gamma_uv_i(h) = b_uv_i * gamma_i(h), where gamma_i(h) is the same basic admissible model used for each variable at that scale, and the matrix of co-regionalization coefficients [b_uv_i] must be positive semi-definite. This ensures that the joint covariance matrix across all variables and all locations is also positive semi-definite.

The LMC framework is widely used in co-kriging and co-simulation workflows where a secondary variable with dense spatial coverage (such as a seismic-derived acoustic impedance cube) is used to improve the spatial interpolation of a primary variable known only at well locations (such as porosity measured from cores or wireline logs). Because the admissibility of each component is guaranteed by additivity, the co-simulation honours the correlation structure between the two variables at every scale simultaneously.