Anisotropy: range, sill, geometric vs zonal

§3.3 introduced the directional variogram $\hat\gamma(h; \phi)$ and ran the standard four-direction sweep that lets you SEE direction-dependent spatial structure in a dataset. The §3.3 ellipse widget then closed the loop by fitting ($\theta, a_{\max}, a_{\min}$) to four observed ranges — a hint that 2-D anisotropy is parsimoniously described by just three numbers. §3.4 takes the next step. It develops the formal anisotropy framework that downstream chapters depend on, makes the geometric / zonal distinction precise, walks through the full anisotropy-extraction workflow with both widgets supporting the reader at each step, and lays out the diagnostic pitfalls that produce silently wrong anisotropy conclusions on real datasets.

The §3.4 message has three parts. First, geometric anisotropy is a CHANGE OF COORDINATES on an isotropic model: rotate the frame so the major axis lines up with the first coordinate, then divide by the (major, minor) ranges, and the variogram of the anisotropic field equals an isotropic variogram evaluated in that reduced frame. Three numbers ($\theta, a_{\max}, a_{\min}$) suffice in 2-D; five numbers do it in 3-D. Second, zonal anisotropy is genuinely different: it requires NESTED structures whose individual sills add up to different directional totals, and in real datasets it usually signals non-stationarity rather than a property of one stationary field. Third, extracting anisotropy from a finite dataset is a workflow with well-known traps — spurious-zonal effects from finite-domain sampling, the masquerade of a missed trend as anisotropy, ellipse fits that are poorly constrained when sample density is sparse in oblique directions. The first §3.4 widget makes the geometric / zonal distinction tangible by showing two synthetic fields side-by-side. The second §3.4 widget walks the full workflow from synthetic field to fitted ($\theta, a_{\max}, a_{\min}$) and shows the reader exactly how well the fit recovers the latent truth.

Geometric anisotropy as a reduced-distance change of coordinates

The clean statement of geometric anisotropy is this: an anisotropic variogram is an ISOTROPIC variogram evaluated on a rotated, stretched lag. Concretely, let $\gamma_0(\cdot)$ be an isotropic variogram model (Spherical, Exponential, Gaussian — anything from Part 4) with sill $\sigma^2$ and (isotropic) range 1. Then the anisotropic variogram is

where $R(\theta)$ is the 2-D rotation that takes the major axis to the first coordinate axis, and $S = \mathrm{diag}(a_{\max}, a_{\min})$ is the range matrix. The combination $R^\top S^{-1}$ does two things in sequence: rotate the lag vector $\mathbf{h}$ so the major axis is horizontal, then divide its first and second components by the major and minor ranges respectively. Compute the Euclidean norm of the result, and that scalar is the REDUCED DISTANCE $h'$. The anisotropic variogram at lag $\mathbf{h}$ equals the isotropic variogram at the scalar reduced distance $h'$.

The closed-form 2-D reduced distance is

where $h = |\mathbf{h}|$ and $\phi = \arg(\mathbf{h})$. The 1-D directional RANGE at azimuth $\phi$ — the value of $h$ at which the isotropic $\gamma_0$ reaches its sill in that direction — is then the reciprocal of the bracketed term:

This is the ANISOTROPY ELLIPSE. As $\phi$ sweeps from 0 to $\pi$, $a(\phi)$ traces out one half of an ellipse with semi-axes $a_{\max}$ and $a_{\min}$ at angle $\theta$ (the full $\phi \in [0, 2\pi)$ traces the ellipse twice because the variogram is direction-reversal symmetric). The three numbers ($\theta, a_{\max}, a_{\min}$) describe this ellipse completely (Goovaerts 1997 §4.4; Chilès & Delfiner 2012 §2.5.2; Isaaks & Srivastava 1989 ch. 7).

The SILL of the anisotropic variogram is unchanged from the isotropic $\gamma_0$. At every direction $\phi$, as $h \to \infty$, $h' \to \infty$ as well (because $a(\phi)$ is finite), so $\gamma_0(h') \to \sigma^2$. Geometric anisotropy stretches the CORRELATION structure but leaves the underlying TOTAL VARIANCE alone. This is the diagnostic that distinguishes geometric from zonal anisotropy at the empirical level: a set of directional variograms that all rise to the SAME sill but at different rates is geometric; a set of directional variograms that rise to DIFFERENT sills is zonal.

In 3-D the reduced-distance form generalises directly:

where $R$ is now a 3-D rotation matrix and $(a_1, a_2, a_3)$ are the three principal-axis ranges. A 3-D rotation has three free angles (typically described in GSLIB and SGeMS as the azimuth, the dip, and a final rake / roll angle; conventions vary by software so always check the package documentation). With three ranges plus the rotation, 3-D geometric anisotropy has six free parameters in the most general case; a common simplification is to assume one of the three axes is vertical (no dip / no rake), reducing the rotation to a single azimuth angle and giving four free parameters total. This is the standard parameterisation used in stratigraphic reservoir modelling and most mining geostatistics (Deutsch & Journel 1998 §II.1.4).

Zonal anisotropy as nested structures and a non-stationarity flag

Zonal anisotropy is qualitatively different. It cannot be reduced to a single isotropic model evaluated in a stretched frame, because no rescaling can change the SILL of an isotropic model. The standard way to write a zonal-anisotropic variogram is as a NESTED sum of structures:

where $\gamma_{\rm iso}$ is an isotropic (or geometrically-anisotropic) structure with sill $\sigma^2_{\rm common}$ that contributes in every direction, and $\gamma_{\rm zonal}$ is a structure whose sill $\sigma^2_{\rm zonal}$ contributes ONLY along a single principal direction (the variability is "striped" along one axis). The reduced-distance argument of $\gamma_{\rm zonal}$ depends only on the along-axis component of the lag; perpendicular displacements see zero contribution from this term. Net effect: the variogram sill along the zonal direction is $\sigma^2_{\rm common} + \sigma^2_{\rm zonal}$; perpendicular sills are $\sigma^2_{\rm common}$ alone (Goovaerts 1997 §4.4.1; Wackernagel 2003 ch. 5; Chilès & Delfiner 2012 §2.5.3).

Where does zonal anisotropy come from physically? Almost always from a system that has been INADEQUATELY DECLUSTERED INTO STATIONARITY DOMAINS. The textbook case is a 3-D reservoir whose vertical variogram reaches a sill of, say, $\sigma^2_{\rm tot} = 1.0$ (the full between-and-within layer variance) but whose horizontal variograms only reach $\sigma^2_{\rm within} = 0.6$ — because each horizontal slice is contained within a single stratigraphic layer and the BETWEEN-layer variance is not sampled by any horizontal pair. Reporting this as zonal anisotropy is technically correct, but the better description is that the reservoir was incorrectly treated as a single stationarity domain: the layers should have been declustered first (Isaaks & Srivastava 1989 ch. 16). Once you condition on the stratigraphic context, the residual horizontal and vertical variograms typically reach the same sill — apparent zonal anisotropy is gone.

For this reason, when you see zonal anisotropy in an EDA report, the diagnostic question is not "what nested-structure model do I fit?" but "what stratigraphic / facies / process boundary is my analysis domain crossing?" Sometimes the answer is "no boundary, this is genuinely a zonal field with an along-axis stripe structure" — striped sand-dune fields and fault-controlled mineralisation are real examples. More often the answer is "I have not yet partitioned my domain correctly". Either way, zonal anisotropy is a flag to investigate rather than a parameter to fit blindly.

In the §3.4 widgets below, the first compares a strictly geometric anisotropic field with a strictly zonal one constructed by adding an along-major-axis nested structure. The two fields share the same major direction and ranges; only the sill structure differs. Toggling between them makes the diagnostic concrete: the SILL is what tells you which kind of anisotropy you are dealing with.

Geometric vs zonal — side by side

The first §3.4 widget. Two anisotropic Gaussian fields are generated on the unit square, both with major direction $\theta = 30^{circ}$ and major / minor ranges $a_{\max} = 0.40$, $a_{\min} = 0.12$. The GEOMETRIC field uses an anisotropic Spherical model with a uniform sill of 1.0 in every direction. The ZONAL field is the sum of an isotropic Spherical structure (sill 1.0, range = $a_{\min}$) and a "striped" Spherical structure (sill 1.0, range = $a_{\max}$) whose argument depends only on the along-major component of the lag — so the zonal field has a sill of 1.0 across-axis but 2.0 along-axis. The reader toggles between the two fields and watches the four directional variograms rearrange themselves.

Three observations to make with this widget. First, the FIELD HEATMAPS look qualitatively similar — both show stretched correlation along the $\theta = 30^{circ}$ direction — but the zonal field has noticeably more vertical "amplitude" along the major axis because of the extra striped sill. The eye does not distinguish geometric from zonal at the visual level reliably; the variogram is what reveals the distinction. Second, the GEOMETRIC variogram panel shows four directional variograms rising to roughly the same asymptote (the sill ≈ 1.0) at different rates. The sill ratio across the four directions is close to 1:1 (small deviations come from finite-sample noise — try clicking around to convince yourself the ratio stays near 1). Third, the ZONAL variogram panel shows the major-direction variograms (the 0° and 45° curves, both close to the $\theta = 30^{circ}$ axis) rising to a sill near 2.0, while the perpendicular variograms (90° and 135°) saturate near 1.0. The sill ratio across directions is now clearly above 1 — typically ~2:1 — which is the empirical signature of zonal anisotropy.

The diagnostic message is the sill ratio readout. If max(sill_φ) / min(sill_φ) is within a few percent of 1, the data is consistent with geometric anisotropy (with the small deviation attributable to sampling noise). If the ratio is clearly above 1 — say above 1.25 — the data is showing zonal behaviour, and the next move is to ask WHY the sill differs across directions before fitting a model. Is the analysis domain crossing a stratigraphic boundary? Is one direction systematically longer than another so it samples a wider stationarity context? Has a trend been residualised? Answers to these questions decide whether the right next step is a nested-structure fit (rare) or a re-partitioning of the analysis domain into separate stationarity sub-domains (common).

The full anisotropy-extraction workflow

Given a dataset suspected of geometric anisotropy, the standard practitioner workflow is:

• Compute $\hat\gamma(h; \phi_k)$ in $k = 4$ to $k = 8$ directions. In 2-D, four directions at 0°, 45°, 90°, 135° with ±22.5° tolerance is the minimum (the §3.3 widget). Six directions at 30° spacing with ±15° tolerance is sharper if sample density allows. Eight at 22.5° spacing with ±11.25° is sharper still but quickly runs into per-bin sample-size problems. Pick the highest resolution your $N(h; \phi)$ counts can sustain by the rule-of-30 (§3.2).
• Identify the candidate major-axis direction. The direction whose variogram has the LONGEST practical range — the lag at which it first crosses ~95% of its sill — is the candidate. Note that this is just a candidate; the true major axis will rarely fall EXACTLY on one of the binned directions.
• Check the perpendicular direction. The variogram in the direction PERPENDICULAR to the candidate should have the SHORTEST practical range. If it does not, something has gone wrong: either the anisotropy is not geometric (rotate-and-stretch model fails), or the candidate major direction is wrong, or the sample density in oblique directions is too sparse to constrain the fit.
• Fit a 2-D ellipse to the {(φ_k, a_k)} pairs. Three parameters ($\theta, a_{\max}, a_{\min}$) fitted to 4–8 observations is well-posed. A grid search over $\theta \in [0^{circ}, 180^{circ})$ and $(a_{\max}, a_{\min}) \in [0, R_{\max}]$ at fine resolution recovers the fit in <100 ms — the second §3.4 widget below does exactly this. More sophisticated estimators (weighted least squares with the variogram-cloud variances as weights; Cressie 1985 WLS) refine the precision but the grid search is the conceptual default.
• Report the four numbers: $\theta, a_{\max}, a_{\min}$, and the derived anisotropy ratio $a_{\max} / a_{\min}$. The ratio is the headline; the absolute ranges feed into Part 5 kriging neighbourhood sizes; the orientation $\theta$ feeds into the kriging search ellipse.
• Interpret the result geologically. Does $\theta$ align with the paleocurrent, bedding strike, fault trend, or process direction inferred from independent geology? Does the anisotropy ratio fit the expected range for the depositional / structural setting (see the table in the next subsection)? An anisotropy result that aligns with independent geology is far more defensible than one that does not.

The second §3.4 widget walks this entire workflow on a synthetic dataset where the latent ($\theta, a_{\max}, a_{\min}$) are known exactly. The reader sets the latent direction and ratio with sliders and then sees the recovered fit alongside the truth — a controlled experiment in how well an empirical ellipse fit recovers the latent geometry from a finite sample.

Four observations. First, with the default parameters ($\theta = 60^{circ}$, ratio = 4:1) and 180 samples, the fit typically recovers $\theta$ to within 5°–15° of the truth and the ratio to within 0.5 of the truth. Click "New realisation" several times: the fit jitters but stays in the right neighbourhood. Second, push the ratio slider down toward 1:1 (isotropic). The fit recovers ratio ≈ 1 reliably, but $\theta$ becomes meaningless — there is no real major axis when the field is isotropic, and the grid search just picks the angle that minimises a noisy residual. Always pair the ratio with the orientation; if the ratio is near 1, report "isotropic" and ignore $\theta$. Third, push the ratio up toward 6:1 with $\theta = 60^{circ}$. The recovery is now extremely clean — strong anisotropy with adequate sample density is the easy case. Fourth, try $\theta$ near 90° vs $\theta$ near 45°. The fit at the cardinal direction (0°, 90°) is sometimes slightly cleaner than at the oblique direction because more pair-counts fall into the cardinal directional bins; this is a small effect with 6-direction binning but worth knowing about.

The workflow as a whole is the gateway from RAW DATA to MODEL PARAMETERS. By the end of this widget, the reader has three numbers that the Part 4 model-fitter and the Part 5 kriging algorithm can consume. The §3.4 toolkit is what bridges the experimental variogram of Parts 3.1–3.3 with the model-based prediction of Parts 4 and 5.

Typical anisotropy ratios across rock-types

An anisotropy ratio that aligns with the expected range for the depositional / structural setting is much more credible than a ratio that does not. The table below is a rough guide; specific values vary widely with scale, sampling design, and how the practitioner defines practical range. Sources: Goovaerts (1997) §4.4 industry-default examples; Pyrcz & Deutsch (2014) §3 reservoir-modelling case studies; Isaaks & Srivastava (1989) chapters 7 & 16 mining-grade datasets; Chilès & Delfiner (2012) §2.5 multi-domain reference.

• Sands and dunes (2-D plan-view): 1.5:1 to 3:1 along the prevailing wind / current direction. The grain-size and sorting fabric is mild, so the variogram ratio rarely exceeds 3:1 in the horizontal plane. 3-D vertical / horizontal ratios for the same systems can be 5:1 or higher.
• Fluvial channels (plan-view): 3:1 to 5:1 along the paleocurrent direction; can reach 10:1 for low-sinuosity straight channels. The major-axis range is the average length of a single sand body; the minor-axis range is the width across the channel.
• Tidal channels and sinuous fluvial: 5:1 to 10:1, sometimes higher. The two-point variogram begins to break down at very high ratios because it cannot represent the CONNECTIVITY of the meander network — this is one of the motivations for multipoint statistics (Part 9). Still, the AXES of the anisotropy come out cleanly even at 10:1.
• Shoreface bars and barrier-island deposits: 3:1 to 5:1 along the depositional strike (parallel to the paleoshoreline). The major-axis range corresponds to bar continuity along strike; the minor-axis range is the dip-direction width.
• Bedded sedimentary reservoirs in 3-D: Horizontal / vertical ratios of 10:1 to 100:1. The horizontal-plane 2-D anisotropy is often mild (within-layer transport is roughly isotropic) but the cross-bedding ratio is extreme. This is the case where 2-D plan-view analysis dramatically UNDERESTIMATES the true 3-D anisotropy.
• Vein-hosted ore deposits: 10:1 to 50:1, with the major axis along the structural fabric (vein strike). The minor-axis range can be smaller than the closest sample spacing — across the vein the analyst is in the nugget regime even though along the vein the variogram is well developed.
• Stratiform iron / copper ore bodies: 2:1 to 5:1 in 2-D plan view (the within-layer mineralisation is roughly uniform in plan view); 3-D ratios with the vertical can be 10:1 or much higher.
• Groundwater contamination plumes: 5:1 to 50:1 along the regional flow direction, with the ratio set by the ratio of longitudinal to transverse dispersion in the aquifer. The flow direction is independently inferable from the head map — a strong consistency check.
• Vegetation, soil, and surface-process variables: Mild 1.5:1 to 3:1, with the major axis along the dominant gradient (slope, aspect, drainage). Often weaker than geological anisotropy, partly because the underlying processes are themselves more spatially distributed.

Two principles to extract from the table. (i) Anisotropy ratios above ~5:1 are common in geology but rare in environmental / surface-process datasets. If your environmental study reports a 20:1 ratio, the result is suspect — check for residual trends, finite-domain effects, or sample-density bias before believing it. (ii) The HEADLINE NUMBER for any geostatistical study is the anisotropy ratio combined with the orientation $\theta$, with the orientation tied to a named geological / process feature. "Anisotropy ratio 4:1 along the inferred N30E paleocurrent direction" is a defensible result. "Anisotropy ratio 4:1" without a direction or a story is not.

Three pitfalls that quietly corrupt anisotropy analyses

Three failure modes account for most bogus anisotropy reports in practice. Each has a diagnostic and a fix.

• Insufficient pair counts in oblique directions. With sparse data, the angular cone for a 45° or 135° direction can have far fewer pairs than the 0° or 90° cones, simply because the sample geometry is not random in azimuth (most sample grids are rectangular). The oblique-direction empirical variograms are then NOISY and can swing the ellipse fit toward a spurious orientation. Diagnostic: always report $N(h; \phi)$ alongside $\hat\gamma(h; \phi)$; if oblique-direction counts are much smaller than cardinal-direction counts, treat the oblique-direction ranges as uncertain. Fix: widen the angular tolerance (back from ±15° to ±22.5°) to gain pairs, or report the fit with an uncertainty bar instead of a single point estimate.
• Nonstationarity masquerading as anisotropy. A field with a deterministic trend — say, a linear gradient $Z(x, y) = ax + by + \varepsilon(x, y)$ — produces directional variograms that look anisotropic even when the residual $\varepsilon$ is isotropic. The variogram in the direction of the trend grows without bound; perpendicular to it the variogram is bounded; the ratio looks like extreme anisotropy. Diagnostic: compute the mean of the field in moving windows along each direction; if the mean changes with location, you have a trend, not just anisotropy. Fix: residualise the trend (linear regression, polynomial fit, kriging-with-trend in Part 5.3) and compute the variograms on the residuals. The residual variograms often reveal a much milder underlying anisotropy than the raw ones.
• The spurious-zonal effect. Even with a perfectly stationary field, the apparent SILL of a directional variogram can be biased upward when the analysis domain is much LONGER in that direction than in others — the long-direction pairs sample a wider stratigraphic context, so they see more total variance than the short-direction pairs. This produces an artefactual zonal anisotropy that vanishes if you reanalyse a square sub-region of the domain. Diagnostic: always report the domain dimensions along each principal direction; if one dimension is >2× another, the spurious-zonal effect is in play. Fix: cap the maximum lag in each direction at half the SHORT-direction dimension of the domain — the standard half-domain rule from §3.2 generalised to anisotropic settings (Goovaerts 1997 §4.3.5; Isaaks & Srivastava 1989 ch. 7).

The discipline that distinguishes good anisotropy analysis from bad analysis is paranoia about each of these three traps. Every reported anisotropy result should arrive with explicit answers to: how many pairs went into each direction? has the trend been residualised? is the domain anisotropy ratio consistent with the data anisotropy ratio? Anisotropy results that pass all three checks are far more credible than those that do not.

Anisotropy as a fingerprint of geological process

An anisotropy ellipse with no geological interpretation is just three numbers. An anisotropy ellipse tied to a depositional or structural process is one of the most powerful pieces of information geostatistics extracts from a dataset. Three quick examples cement the point.

Fluvial sandstone reservoir. A horizontal-plane analysis recovers $\theta = 25^{circ}$ E of N, $a_{\max} = 350$ m, $a_{\min} = 90$ m. Anisotropy ratio is 3.9:1, with the major axis pointing along the inferred paleocurrent direction from independent core data (paleocurrent indicators, ripple cross-bedding orientations). The major-axis range = 350 m is the average length of a single sand body along its long axis; the minor-axis range = 90 m is the across-channel width. Kriging with this ellipse produces elongated prediction zones in the N25E direction that honour the depositional fabric and produce reservoir models materially different from an isotropic baseline.

Vein-hosted gold deposit. A horizontal-plane analysis recovers $\theta = 75^{circ}$, $a_{\max} = 80$ m, $a_{\min} = 4$ m. Anisotropy ratio is 20:1, with the major axis aligned with the local vein-strike from structural mapping. The minor-axis range is smaller than the closest drillhole spacing — across the vein, the grade is effectively unpredictable from neighbouring samples. Kriging without this anisotropy produces a smeared grade map that underestimates the high-grade narrow stripe and overestimates the low-grade halo; the resulting resource estimate is silently biased low for the recoverable ounces.

Groundwater contamination plume. A 2-D plume of dissolved TCE in a sandy aquifer shows $\theta = 105^{circ}$, $a_{\max} = 250$ m, $a_{\min} = 30$ m. Anisotropy ratio is 8.3:1, with the major axis pointing along the regional groundwater-flow direction inferred from the local head map. The anisotropy ratio matches the literature value for longitudinal/transverse dispersion in clean sands (5:1 to 15:1). This independent consistency strengthens the geostatistical model and lets it be used for plume-front predictions with greater confidence than an isotropic baseline would warrant.

In each case the workflow is the same: extract the anisotropy ellipse from the geostatistics, compare to the independent geological / hydrogeological inference, and report only the result that survives both. This is the discipline that separates anisotropy analysis as a credible science from anisotropy analysis as a parameter-fitting exercise.

Where §3.4 ends and §3.5–§3.6 + Part 4 pick up

§3.4 has given you the framework. You can write an anisotropic variogram in reduced-distance form; you can distinguish geometric from zonal anisotropy by the sill ratio across directions; you can run the full ellipse-extraction workflow on a finite dataset and report ($\theta, a_{\max}, a_{\min}$) with awareness of the three classic pitfalls; you have a rough sense of what anisotropy ratios to expect across rock-types. What §3.4 has NOT done is fit a permissible MODEL to your anisotropic empirical variogram. The §3.4 widgets show empirical points and a fitted ellipse; the ellipse is a graphical aid, not a continuous variogram function that the kriging system can consume.

§3.5 swaps the continuous $Z(\mathbf{s})$ for an indicator $I(\mathbf{s}; z_c) = \mathbf{1}{Z(\mathbf{s}) \le z_c}$ and runs the same directional / anisotropic binning on the indicator — the bridge to indicator kriging and categorical-facies modelling in Part 8. §3.6 develops the robust Cressie-Hawkins estimator that hardens the squared-difference $\hat\gamma$ against outliers; the same anisotropic-binning machinery applies. Part 4 then closes the experimental-and-model loop: §4.1 develops permissible isotropic models (Spherical, Exponential, Gaussian, power, Matérn); §4.3 fits anisotropic versions of those models to the directional empirical variograms produced here; §4.4 develops nested structures including the formal zonal-anisotropic models that this section's second widget hints at. Part 5 kriging finally consumes the ($\theta, a_{\max}, a_{\min}$, sill, nugget) parameter vector to produce optimal predictions.

Try it

• In the anisotropy-types widget, click "Geometric". Read off the four directional sills from the readout pills. They should all be close to 1.0 (within ~10% of each other). Now click "Zonal". The 0° and 45° sills should rise to roughly 2.0 (these directions are closest to the major axis $\theta = 30^{circ}$), while 90° and 135° stay near 1.0. Note that the SILL RATIO max:min pill is the diagnostic: ~1:1 in geometric mode, ~2:1 in zonal mode.
• Same widget, zonal mode. The major-axis direction is θ = 30°. Of the four canonical directions (0°, 45°, 90°, 135°), which one is CLOSEST to the major axis? Which is FURTHEST? Use the field heatmap to confirm visually. (Hint: 30° is between 0° and 45°, so both 0° and 45° see partial along-major variance; 120° = θ + 90° is the perpendicular direction, closest to 135°; 90° is also somewhat perpendicular.)
• In the anisotropy-fit-from-data widget, set latent θ = 60° and anisotropy ratio = 4. The recovered θ_fit should be within ~10° of 60° and the ratio within ~0.5 of 4. Note the |θ_fit − θ_latent| pill. Click "New realisation" five times in a row and watch how the fit jitters. The orientation jitter is the natural uncertainty of an ellipse fit with 180 samples.
• Same widget. Lower the anisotropy ratio to 1.0 (isotropic). Watch what happens to the fitted θ. It will JUMP between widely different values across realisations — the field has no preferred direction, so the fit picks essentially arbitrary orientations. The takeaway: anisotropy ratio close to 1 means the orientation is meaningless; report "isotropic" rather than θ.
• Same widget. Raise the anisotropy ratio to 6.0 with latent θ = 60°. The fit becomes very clean — θ_fit stays close to 60° across realisations and the ratio stays close to 6. Strong anisotropy with adequate sample density is the easy case; this is what successful anisotropy analyses look like.
• Without coding: you analyse a horizontal-plane variogram in four directions on a 600 m × 200 m dataset (the long dimension is east-west). You find the east-west sill = 1.4 and the north-south sill = 1.0, both reaching their sills at h = 80 m. Is this geometric or zonal anisotropy? What other explanation should you check FIRST? (Hint: sills differ → zonal; but the domain is 3× longer in east-west than north-south, so the SPURIOUS-ZONAL pitfall is in play. Re-analyse a 200 m × 200 m square sub-region before claiming zonal anisotropy.)
• Without coding: a fluvial-channel reservoir study reports anisotropy ratio 12:1 with the major axis at $\theta = 75^{circ}$ E of N, while independent paleocurrent indicators point to N20E. What's most likely going on? (Hint: the geostatistical θ and the geological paleocurrent disagree by ~85° — likely a missed trend in the data or a misidentification of the major direction by an under-sampled oblique direction. Re-run the workflow on residualised data and check the oblique-direction pair counts before reporting the ratio.)
• Without coding: an environmental-soil study reports an anisotropy ratio of 25:1 for soil moisture across a 1 ha field. Is this credible? (Hint: 25:1 is at the high end even for vein-hosted ore; in environmental settings ratios above ~5:1 are suspect. Check for surface-drainage or aspect-slope trends; residualise; recompute. The mild residual anisotropy is the real story.)

Pause and reflect: the anisotropy ellipse compresses an enormous amount of spatial information into three numbers. Those three numbers can be defensible or they can be artefacts of pair-counting, trend, and domain shape. The discipline of §3.4 is to test which.

What you now know — and what §3.5 onward builds on it

You can write the anisotropic variogram in reduced-distance form: $\gamma(\mathbf{h}) = \gamma_0(|R(\theta)^\top S^{-1} \mathbf{h}|)$, where $R(\theta)$ is the rotation that aligns the major axis with the first coordinate and $S$ is the diagonal range matrix. Three numbers ($\theta, a_{\max}, a_{\min}$) describe 2-D geometric anisotropy; five-to-six numbers do it in 3-D. You can distinguish geometric anisotropy (same sill across directions) from zonal anisotropy (different sills), and you know that zonal anisotropy in real data is almost always a non-stationarity diagnostic — a flag to repartition the analysis domain rather than to fit a more complex model.

You can run the full workflow from data to ellipse: compute $\hat\gamma(h; \phi_k)$ in 4–8 directions, extract a practical range per direction, fit a 2-D ellipse to the ($\phi_k, a_k$) pairs by grid search, and report ($\theta, a_{\max}, a_{\min}$) plus the derived ratio. You have a rough sense of what ratios to expect across rock-types (sands 1.5–3:1, fluvial 3–5:1, tidal 5–10:1, vein-hosted ores 10–50:1, environmental fields rarely above 5:1) and you know that the headline number is the ratio combined with an orientation that points at a named geological / process feature.

You can spot the three classic pitfalls: insufficient pair counts in oblique directions, nonstationarity masquerading as anisotropy, and the spurious-zonal effect from finite-domain sampling. Each has a diagnostic and a fix; the discipline of paranoia about each is what separates defensible anisotropy analysis from parameter-fitting theatre.

§3.5 reuses this directional / anisotropic binning machinery on the indicator $I(\mathbf{s}; z_c)$ — the bridge to categorical-facies modelling and indicator kriging in Part 8. §3.6 hardens the squared-difference estimator against outliers with the Cressie-Hawkins 1980 robust form. Part 4 fits permissible anisotropic models to the empirical variograms you produced here. Part 5 kriging finally consumes ($\theta, a_{\max}, a_{\min}$, sill, nugget) to produce optimal predictions. The anisotropy framework developed in this section is the conceptual hinge between the experimental variograms of Part 3 and the model-based prediction of Parts 4 and 5.

Mark section complete →

References

• Matheron, G. (1971). The Theory of Regionalized Variables and Its Applications. Les Cahiers du Centre de Morphologie Mathématique, No. 5. Fontainebleau, France: École Nationale Supérieure des Mines de Paris. (The canonical reference for the reduced-distance / anisotropic generalisation of $\gamma$ in the regionalised-variable framework. The rotation-and-stretch parameterisation that this section formalises traces to Matheron's original development.)
• Cressie, N. (1993). Statistics for Spatial Data (revised ed.). Wiley. (§2.5 "Anisotropy" — the modern mathematical-statistics reference for the distinction between geometric and zonal anisotropy and the conditions under which each is a defensible stationary-field model. The definitive treatment of the limits of the two-point variogram framework for highly anisotropic data.)
• Goovaerts, P. (1997). Geostatistics for Natural Resources Evaluation. Oxford University Press. (§4.4 on anisotropic variograms — the practitioner reference whose organisation §3.4 most closely follows. The reduced-distance parameterisation $h' = h\sqrt{(\cos(\phi-\theta)/a_{\max})^2 + (\sin(\phi-\theta)/a_{\min})^2}$, the geometric / zonal taxonomy, and the typical-ratio table all draw on this source.)
• Isaaks, E.H., Srivastava, R.M. (1989). An Introduction to Applied Geostatistics. Oxford University Press. (Chapters 7 and 16 on spatial continuity and 3-D anisotropy — the most readable practitioner-level treatment. Chapter 16's discussion of horizontal-vs-vertical anisotropy in bedded reservoirs is the source of this section's zonal-anisotropy / stratigraphic interpretation.)
• Deutsch, C.V., Journel, A.G. (1998). GSLIB: Geostatistical Software Library and User's Guide (2nd ed.). Oxford University Press. (§II.1.4 documents the 3-D anisotropy rotation convention used by the GSLIB programs vmodel, kt3d, sgsim, and others. The convention for the three rotation angles (azimuth, dip, rake) is industry-standard in mining geostatistics.)
• Chilès, J.-P., Delfiner, P. (2012). Geostatistics: Modeling Spatial Uncertainty (2nd ed.). Wiley. (§2.5.2 on anisotropic variograms, including the formal characterisation of geometric and zonal anisotropy and the conditions under which each is a permissible stationary model. The modern comprehensive reference.)
• Wackernagel, H. (2003). Multivariate Geostatistics: An Introduction with Applications (3rd ed.). Springer. (Chapter 5 on the linear model of regionalisation, which provides the formal framework for the nested-structure decomposition of zonal anisotropy. The cleanest mathematical treatment of why zonal anisotropy requires multiple nested structures.)
• Pyrcz, M.J., Deutsch, C.V. (2014). Geostatistical Reservoir Modeling (2nd ed.). Oxford University Press. (Reservoir-engineering treatment that emphasises the connection between the anisotropy ellipse and depositional process in fluvial, deltaic, shoreface, and carbonate settings. The anisotropy-vs-rock-type table in this section draws on case studies from this reference.)