Isotropic and directional variograms

§3.1 and §3.2 worked under a quiet assumption that is almost never strictly true for real geological data: that the variogram $\gamma$ depends only on the SCALAR lag distance $h = |\mathbf{h}|$, not on the direction of the lag vector. A pair of samples 100 m apart in the east-west direction was treated as carrying the same information about spatial correlation as a pair 100 m apart in the north-south direction. This is the isotropy assumption, and it is the simplest possible spatial-correlation model. It is also the one geology breaks first.

Real geological fields almost always have a preferred direction. Fluvial sandstone bodies are long along the paleocurrent direction and short across it. Bedded sediments are continuous along the bedding plane and decorrelate quickly perpendicular to it. Mineral veins follow structural trends. Pollution plumes spread along groundwater flow. The list goes on: the depositional, tectonic, hydrodynamic, and chemical processes that produced the rock leave a directional imprint, and the variogram you would measure along the long axis of that imprint differs — sometimes dramatically — from the one you would measure across it. §3.3 lifts the isotropy assumption and develops the machinery for measuring direction-dependent spatial structure from data.

The story is short and the algebra is small. The §3.2 Matheron estimator already does almost all the work; you only need to add a directional filter on the pairs you accept into each bin. Two widgets stage the story: the first GENERATES a 2-D anisotropic field with tunable major range, minor range, and orientation, and computes the four directional variograms so you can watch them separate as you crank the anisotropy ratio; the second takes the four observed ranges and FITS a 2-D anisotropy ellipse so you can read off the three numbers (major axis, minor axis, orientation) that fully describe geometric anisotropy in 2-D. By the end you have a defensible experimental description of direction-dependent spatial structure, and §3.4 will codify that description into the anisotropy-ellipsoid model that Part 4 fits and Part 5 kriging consumes.

Why isotropy fails for real geology

The isotropy assumption is mathematically clean. It is the geological reality that fails, and it fails in a structured way that makes the failure detectable from data. Three families of process produce direction-dependent spatial structure in nearly every reservoir-characterisation, mining, hydrogeology, or environmental dataset you are likely to encounter.

• Sedimentary fabric. Deposition lays down rock in layers, with grain orientation, sorting, and lithification varying systematically along the bedding plane and perpendicular to it. A horizontal sand layer in a reservoir is laterally continuous over tens to hundreds of metres but only metres thick vertically. The variogram along any horizontal direction reaches its sill at much larger $h$ than the variogram along the vertical. Even WITHIN the horizontal plane, if the depositional system had a directional flow (a river, a delta lobe, a tidal current), the sediment fabric will be longer along the paleocurrent direction than across it.
• Structural fabric. Faults, fractures, folds, and joint sets impose linear or curvilinear directions on the rock mass, and properties that are sensitive to those features — fracture density, mineralisation grade, hydraulic conductivity — vary slowly along the fabric direction and fast perpendicular to it. A vein-hosted gold deposit is the textbook case: the grade variogram along the vein strike has range tens to hundreds of metres; perpendicular to the vein it has range of a metre or less.
• Process-driven fabric. Even where the rock itself is approximately isotropic, the PROCESSES that operate on it can impose directional structure. Groundwater contamination spreads along the regional flow direction. Vegetation strips orient along moisture or aspect gradients. Acid mine drainage migrates along surface drainage. In each case the variogram of the affected variable is anisotropic with the major axis along the process direction.

The practical consequence is that in any real geostatistical study, you should EXPECT anisotropy and let the data tell you whether it is mild enough to ignore. The isotropic γ(h) of §3.1–§3.2 is then a special case — what you get when the four directional variograms happen to overlap — rather than the default. §3.3 is how you check.

The directional variogram γ(h; φ) and how to compute it

The directional variogram is the population analogue of the §3.1 definition restricted to a single direction:

where $\mathbf{h}$ is now treated as a VECTOR with both a magnitude $h = |\mathbf{h}|$ and a direction $\phi = \arg(\mathbf{h})$. The isotropic case is just the special case in which $\gamma(\mathbf{h})$ depends only on $|\mathbf{h}|$. In the general case, $\gamma(\mathbf{h})$ is a function of the full vector. For 2-D in-plane work the common practical parameterisation is

where $\hat{\mathbf{e}}_\phi = (\cos\phi, \sin\phi)$ is the unit vector along direction $\phi$. The directional variogram is the function $\phi \mapsto \gamma(h; \phi)$ at each scalar lag $h$ — a 2-D surface plotted against $(h, \phi)$. The variogram is also direction-reversal symmetric: $\gamma(h; \phi) = \gamma(h; \phi + \pi)$, because the variance of the increment does not depend on which end of the pair is labelled "first". So $\phi$ effectively lives on $[0, \pi)$ — only the unsigned axis matters, not the arrow direction.

From data, you ESTIMATE $\gamma(h; \phi)$ with a directional version of the Matheron 1962 estimator from §3.2. The change from the isotropic version is a single filter on the pair set: a pair $(i, j)$ contributes to the bin centred on direction $\phi_0$ and lag $h_k$ if

• its separation magnitude $|\mathbf{h}_{ij}|$ falls inside the lag bin $[k\Delta, (k+1)\Delta]$ (the §3.2 rule, unchanged), AND
• the angular distance between its separation azimuth $\phi_{ij}$ and the target direction $\phi_0$ is at most a chosen angular tolerance $\Delta\phi$ (the new rule).

Inside the qualifying pair set the estimator is unchanged — squared-increment average with the $1/2$ factor:

The angular tolerance $\Delta\phi$ is a user choice with the same trade-off as the lag bin width $\Delta$: tight cones give cleaner direction-by-direction resolution but smaller $N(h; \phi_0)$ per bin (so noisier $\hat{\gamma}$); wide cones smooth the directional dependence by averaging across a wider azimuth range. The default for 2-D screening is $\Delta\phi = \pm 22.5^{circ}$ — four directions (0°, 45°, 90°, 135°) at this tolerance tile the half-circle $[0, \pi)$ without overlap and without gaps. Once you have identified a preferred direction, tighter tolerances of $\pm 15^{circ}$ or $\pm 10^{circ}$ give a sharper view of the principal-axis variograms (Goovaerts 1997 §4.4, Isaaks & Srivastava 1989 ch. 7, Deutsch & Journel 1998 program gamv).

One refinement matters at long lag: as $h$ grows, even a fixed angular tolerance $\Delta\phi$ corresponds to an ever-WIDENING physical cone (the lateral offset grows linearly with $h$). At very large $h$ the cone can become so wide that the directional filter is effectively meaningless — a pair at $h = 500$ m and $\Delta\phi = 22.5^{circ}$ can be up to 200 m off the target azimuth, which is itself a large spatial offset. The standard fix is the BANDWIDTH cap: a pair also has to fall within a chosen lateral offset $b$ of the target direction. Practically, $b$ is set to roughly the minimum sample spacing or to a fraction of the expected major-axis range. The widget below uses a tolerance-only filter (no explicit bandwidth) for clarity, but real GSLIB gamv runs always set both (Deutsch & Journel 1998 §II.2).

Directional variograms in action

The first §3.3 widget. A 2-D anisotropic Gaussian random field with tunable major range $a_{\max}$, minor range $a_{\min}$, and major-axis direction $\theta$ is generated via Cholesky against an anisotropic Spherical covariance. 120 samples are drawn from the field on the unit square, and the four directional empirical variograms (0°, 45°, 90°, 135°, each with $\Delta\phi = \pm 22.5^{circ}$) are computed and overlaid on a single set of axes. As you slide the anisotropy parameters you can watch the four curves merge or separate in real time.

Three things to do with this widget. First, set both ranges equal ($a_{\max} = a_{\min} = 0.30$): the four directional variograms collapse onto each other and you have recovered the isotropic case of §3.1–§3.2 — a useful sanity check that the directional binning is correctly implemented. Second, set $a_{\max} = 0.50$ and $a_{\min} = 0.10$ with $\theta = 0^{circ}$ (major axis east-west): the 0° curve has the longest range, the 90° curve has the shortest, and 45° and 135° fall between them at roughly equal positions because they make a 45° angle with the principal axes. Third, hold the ranges fixed and slide $\theta$ through 90°: the FAR variogram becomes the NEAR variogram and vice versa — the assignment of "which curve is the major direction" rotates with $\theta$. The four-direction sweep is a SCREEN; the actual major-axis direction will not in general land exactly on one of the four compass-aligned bins, which is why the second widget below fits an ellipse rather than just reading off the longest of the four curves.

The dashed curves overlaid on each directional empirical variogram are the TRUE directional variogram of the generating process — what you would get from infinite samples. The empirical points wiggle around these true curves with the same $1/\sqrt{N(h; \phi)}$ noise scaling as in §3.2. The directional filter cuts the pair count per bin by roughly a factor of four (one direction takes 1/4 of the half-circle), so bins that had $N(h) = 80$ in the isotropic case might have $N(h; \phi) = 20$ in the directional case — borderline by the rule-of-30. This is the cost of resolving direction: you need more samples or a coarser bin width to keep each direction-by-bin combination above the noise floor. In practice you often see directional variograms that are clearly separated in their LONG-LAG behaviour (where the structural signal is strong relative to the noise) but jagged at SHORT-LAG, where small-sample noise dominates each directional bin.

Two faces of anisotropy: geometric and zonal

When the four (or more) directional variograms separate, they can do so in two qualitatively different ways. The distinction matters for what model you build downstream.

Geometric anisotropy. All directional variograms have the SAME sill but different ranges. The four curves rise from zero at different slopes and reach a common asymptote. Geometric anisotropy is the GEOLOGICAL DEFAULT — it is what you get from a stationary field whose correlation is stretched anisotropically without changing the underlying variance. The variogram surface $\gamma(\mathbf{h})$ traced out by the level set $\gamma = c$ for some constant $c <$ sill is an ELLIPSE (in 2-D) or an ELLIPSOID (in 3-D). The two semi-axes of the ellipse are the major range $a_{\max}$ and the minor range $a_{\min}$; the orientation of the major axis is $\theta$. These three numbers fully describe 2-D geometric anisotropy. The standard models are simply isotropic models computed at a REDUCED distance $h' = h \sqrt{(\cos\phi'/a_{\max})^2 + (\sin\phi'/a_{\min})^2}$ with $\phi' = \phi - \theta$ (Goovaerts 1997 eq 4.20, Chilès & Delfiner 2012 §2.5.2).

Zonal anisotropy. Different directional variograms have DIFFERENT sills (in addition, usually, to different ranges). The four curves rise at different rates AND reach different asymptotes. Zonal anisotropy cannot be captured by stretching a single isotropic structure — it requires NESTED structures, one of which has its sill confined to a single direction. The standard model treats the field as the sum of an isotropic component plus a directional component whose variance only contributes along a specific axis: see §4.4 on nested structures. Physically, zonal anisotropy is usually a SIGN of non-stationarity rather than a true property of a stationary anisotropic field. The classic case is a 3-D reservoir whose vertical variogram reaches a sill at 1.0 (the true total variance) but whose horizontal variograms only reach 0.6 — because each horizontal slice is confined to a single stratigraphic layer, and the between-layer variance is not represented in the horizontal pair set (Goovaerts 1997 §4.4.1, Isaaks & Srivastava 1989 ch. 16). Real "zonal" anisotropy almost always disappears once you condition on the stratigraphic context that was masking it.

For 2-D screening with the first widget above, the field is generated with strict GEOMETRIC anisotropy — all four directional variograms reach the same sill $\sigma^2 = 1$, only the ranges differ. This is the case you should expect to see most of the time in real data inside a single facies / single stationarity domain. The second widget below assumes geometric anisotropy and fits the corresponding ellipse.

The anisotropy ellipse — three numbers describe 2-D anisotropy

The second §3.3 widget closes the loop. Given the four observed directional ranges from the empirical variograms (one per bin direction), it fits a 2-D anisotropy ellipse and reports the three numbers that fully describe geometric anisotropy: the major-axis length $a_{\max}$, the minor-axis length $a_{\min}$, and the major-axis orientation $\theta$. The anisotropy ratio $a_{\max} / a_{\min}$ falls out as the derived fourth number.

Adjust the four range sliders to whatever values you read off the previous widget's curves, or use one of the preset buttons (isotropic; mild 1.5:1 anisotropy; strong 3:1 along the north-south axis; strong 4:1 along the NE-SW axis at $\theta \approx 45^{circ}$). The polar rose plot on the left shows the four observed ranges as coloured dots at their compass directions, with the fitted ellipse traced in green. The major axis is drawn as an amber dashed line through the centre at angle $\theta$; the minor axis is drawn as a pink dashed line at $\theta + 90^{circ}$. The readout pills on the right call out the three anisotropy numbers separately, with the anisotropy ratio highlighted. The widget grid-searches over $(\theta, a_{\max}, a_{\min})$ at 2° angular resolution and 0.01 axis-length resolution, which is more than fine enough to fit four observed ranges with three free parameters.

The pedagogical move is the same as the §3.1 anatomy widget: the three numbers ($\theta$, $a_{\max}$, $a_{\min}$) are not abstract fit parameters — they ARE the directional spatial structure of the field. The major axis points along the depositional, structural, or process direction; the major-axis length is how far you can interpolate confidently in that direction; the minor-axis length is how far you can interpolate perpendicular to it. The ratio of the two is how strongly directional the field is. In §4.3 these three numbers will feed directly into the search-ellipse and the anisotropic-Spherical / Exponential / Gaussian models that Part 5 kriging plugs into.

What anisotropy tells you about the rock

The anisotropy ellipse is more than a fitting artefact — it carries direct geological / process information that you can read off the three numbers. A few examples illustrate the connection.

• Fluvial sandstone reservoir. Anisotropy ratio typically 3:1 to 5:1, major axis along the paleocurrent direction. The major-axis range is the average length of a sand body along its long axis; the minor-axis range is the body width across the channel. Kriging with this ellipsoid produces elongated prediction zones aligned with the inferred paleocurrent — a defensible reservoir-model output that honours the depositional process.
• Tidal-channel reservoir. Anisotropy ratio often 5:1 to 10:1, major axis along the channel sinuosity direction (which can be a sweeping average direction if the channels meander). For very high anisotropy ratios, the two-point variogram model starts to fail to capture the connectivity of the channel network — this is one of the motivations for multipoint statistics in Part 9. The variogram still gives you the AXES of the anisotropy correctly; what it loses is the CONTINUITY of the long axis.
• Stratified iron-ore body. Anisotropy ratio in 2-D plan view often near isotropic (1:1 to 2:1, since within a horizontal stratigraphic layer the mineralisation is roughly uniform); 3-D variograms reveal the strong vertical anisotropy with horizontal range tens of times the vertical. This is a case where 2-D analysis underestimates the true anisotropy because the dominant axis is out of the analysis plane.
• Vein-hosted gold deposit. Extreme anisotropy ratio, often 10:1 or higher, major axis along the structural fabric (vein strike). The minor-axis range can be smaller than the closest sample spacing — meaning that across the vein you are essentially in the nugget regime even though along the vein you have well-developed spatial structure. This is the kind of dataset for which directional variogram modelling is not optional; an isotropic model would average the vein-parallel and vein-perpendicular variograms together and produce silently miscalibrated reserves estimates.
• Groundwater contamination plume. Anisotropy ratio set by the ratio of longitudinal to transverse dispersion in the aquifer, often 5:1 to 50:1, major axis along the regional flow direction (which can be inferred from the head map). The anisotropy ellipse of the contaminant concentration is a direct fingerprint of the flow-and-transport process.

In each case the three numbers $(\theta, a_{\max}, a_{\min})$ have a direct physical meaning. Reporting an anisotropy ellipse without naming the process that produced it is a missed opportunity. Reporting it WITH the process is one of the most useful pieces of information geostatistics extracts from a dataset.

Where §3.3 ends and §3.4 picks up

Be precise about scope. §3.3 has shown you (i) why the isotropy assumption almost always fails for real geology, (ii) how to estimate the directional variogram $\hat{\gamma}(h; \phi)$ from data by adding an angular filter to the §3.2 binning, (iii) the distinction between geometric and zonal anisotropy, and (iv) how to compress a set of directional ranges into a 2-D ellipse described by three numbers. What §3.3 has NOT done is fit a permissible MODEL to that ellipse — the four observed ranges are still empirical, the ellipse is a graphical aid, and the directional variogram is still presented as a set of binned data points rather than a continuous function that the kriging system can consume.

§3.4 takes up the modelling problem. It introduces the FULL anisotropy parameterisation that Part 4.3 will fit and Part 5 kriging will consume: the 2-D and 3-D anisotropy ellipsoids written in terms of rotation matrices and diagonal range matrices, the conventions for naming the three (or six) angles in 3-D, the search-ellipse that determines which samples enter a kriging neighbourhood, and the diagnostic plots (rose diagram, variogram map) that practitioners use to identify anisotropy at the EDA stage before any formal fitting. §3.5 then introduces indicator variograms (categorical / threshold-based), §3.6 introduces robust estimators (Cressie-Hawkins and friends), and Part 4 closes Parts 3 and 4 by fitting permissible models — including anisotropic ones — to the experimental variograms produced here.

Try it

• In the directional-variogram-explorer, set $a_{\max} = a_{\min} = 0.30$ (isotropic). The four directional variograms should overlap closely. They will not be IDENTICAL because each direction has finite $N(h; \phi)$ and so each is a separate noisy realisation of the same underlying isotropic γ. How much wiggle around the gold theoretical curve is normal? (Hint: about $\sqrt{2/N}$ as a coefficient of variation, with $N \approx 20$ per bin — so ~30% wiggle.)
• In the directional-variogram-explorer, set $a_{\max} = 0.50$, $a_{\min} = 0.12$, $\theta = 0^{circ}$ (major axis east-west). Look at the four empirical variograms. Which has the longest practical range? Which has the shortest? Are the 45° and 135° curves roughly equal? Why? (Hint: 45° and 135° each make the same angle with both principal axes, so they sample the same intermediate reduced distance.)
• Same widget, same settings. Now click "New realisation" four or five times. The empirical variograms wiggle, but the SEPARATION between the major and minor directions persists across realisations. This is the test of "real anisotropy" vs "noise": real anisotropy is reproducible across realisations; noise is not.
• In the directional-variogram-explorer, hold $a_{\max} = 0.50$ and $a_{\min} = 0.12$ and slide $\theta$ from 0° to 90°. Watch which of the four curves takes the longest range as $\theta$ rotates: at $\theta = 0^{circ}$ it is the 0° curve; at $\theta = 45^{circ}$ it is the 45° curve; at $\theta = 90^{circ}$ it is the 90° curve. The major direction is whatever direction $\theta$ points the major axis to.
• In the anisotropy-ellipse widget, click the "Strong NE-SW" preset. Read off the three fitted numbers. The ratio should be ~4:1 with $\theta$ near 45°. Now move the 0° slider up to 0.30 to break the consistent ellipse — what happens to the fit? (Hint: the RSS rises and the ellipse no longer passes cleanly through all four observed dots. The fit is still the BEST 2-D ellipse, but the four observations are no longer mutually consistent.)
• In the anisotropy-ellipse widget, set all four observed ranges to 0.28 (isotropic). The fitted ellipse becomes a CIRCLE with $a_{\max} = a_{\min} = 0.28$. Note that $\theta$ is then UNDEFINED — any orientation gives the same circle. The grid search will pick some arbitrary $\theta$; in practice you should report "isotropic" rather than $\theta$. The anisotropy ratio is the practical signal: ratio close to 1 ⇒ isotropic.
• Without coding: a colleague brings you an experimental variogram analysis of a fluvial sandstone reservoir with horizontal directional variograms in four directions, all reaching the same sill of 0.04 (porosity variance) but with practical ranges 35 m at 0°, 22 m at 45°, 14 m at 90°, and 21 m at 135°. What anisotropy ratio do they have? Where is the major direction? Is this geometric or zonal anisotropy? (Hint: ratio = 35/14 = 2.5:1; major direction near 0° (east-west) since that range is longest; geometric, since all four reach the same sill.)
• Without coding: a 3-D variogram analysis of an ore-grade dataset reports a vertical sill of 1.2 (the total variance) and four horizontal directional sills of 0.7, 0.8, 0.7, 0.75 with horizontal ranges of 80 m, 110 m, 75 m, 95 m. What kind of anisotropy is this? What does it suggest about the geology? (Hint: the horizontal sills are CONSISTENTLY below the vertical sill — this is ZONAL anisotropy. The vertical direction sees the full between-stratum variance; horizontal slices only see the within-stratum variance. The geology is stratigraphically layered. The ratio of within- to between-strata variance is 0.75/1.2 ≈ 0.6, so about 40% of the variance is BETWEEN layers.)

Pause and reflect: the directional variogram is the place where your knowledge of the depositional or structural process meets the data. If you ALREADY know the major direction from the geology (the channel orientation, the bedding strike, the fault trend), the four-direction sweep is a CONFIRMATION. If you do not, the data tells you. Which way around is more common in your domain — do you usually walk in knowing the major direction and use the variograms to estimate the ranges, or do you walk in blind and use the variograms to discover the direction itself?

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

You have the directional variogram $\gamma(h; \phi)$ as the population generalisation of $\gamma(h)$ that admits direction-dependent spatial structure. You have its empirical estimator $\hat{\gamma}(h; \phi_0)$: bin pairs by lag distance AND by azimuth (cone tolerance $\Delta\phi$, default $\pm 22.5^{circ}$ in 2-D), with an optional bandwidth cap for the long-lag pairs. You know the default 2-D four-direction sweep (0°, 45°, 90°, 135°) and that it is a SCREEN, not a final fit — the major and minor axes typically fall between the four compass directions, and the second widget's ellipse fit recovers their precise orientation. You can distinguish geometric anisotropy (same sill, different ranges — the geological default) from zonal anisotropy (different sills, usually a non-stationarity diagnostic) and you have a vocabulary for reporting which kind your data shows.

You have the three-number description of 2-D geometric anisotropy: major-axis range $a_{\max}$, minor-axis range $a_{\min}$, major-axis direction $\theta$, plus the derived anisotropy ratio $a_{\max}/a_{\min}$. You can read all three off the anisotropy ellipse, and you can connect each one to a depositional / structural / process feature of the rock. You have a working sense of what anisotropy ratios are typical in real geology — mild 1.5:1 for nearshore sediments, moderate 3:1 to 5:1 for fluvial channels and shoreface bars, strong 5:1 to 10:1 for tidal channels and bedded reservoirs, extreme 10:1 or higher for vein-hosted ores and structurally controlled mineralisation.

§3.4 picks up the formal anisotropy ellipsoid: rotation matrices in 2-D and 3-D, the GSLIB angle conventions, the search ellipse for kriging, and the variogram map / rose diagram that you would use at the EDA stage to spot anisotropy before any formal binning. §3.5 swaps the continuous $Z$ for an indicator $I(\mathbf{s}; z_c) = \mathbf{1}{Z(\mathbf{s}) \le z_c}$ and runs the same directional binning on the indicator — the bridge to indicator kriging in Part 8. §3.6 hardens the squared-difference estimator against outliers with the Cressie-Hawkins 1980 robust form. Part 4 then fits permissible anisotropic models to the empirical variograms you produced here. Part 5 kriging finally consumes them. The directional-binning machinery developed in this section is the foundation; everything downstream is an elaboration.

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 directional / anisotropic generalisation of $\gamma$ in the regionalised-variable framework. Defines the reduced-distance approach to anisotropic models.)
• 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 model. The definitive treatment of when the two-point variogram framework breaks for anisotropic data.)
• Goovaerts, P. (1997). Geostatistics for Natural Resources Evaluation. Oxford University Press. (§4.4 on anisotropic variograms — the practitioner-level treatment that this section's organisation most closely follows. The four-direction sweep convention, the $\pm 22.5^{circ}$ tolerance default, and the reduced-distance parameterisation $h' = \sqrt{(\cos\phi'/a_{\max})^2 + (\sin\phi'/a_{\min})^2}$ are all developed here.)
• Isaaks, E.H., Srivastava, R.M. (1989). An Introduction to Applied Geostatistics. Oxford University Press. (Chapter 7 on spatial continuity, especially §7.4 on anisotropy. The most readable practitioner-level treatment of directional binning and the angular-tolerance choice; the rose-diagram visualisation used in this section's widget originates from the GSLIB/Stanford tradition this book represents.)
• Deutsch, C.V., Journel, A.G. (1998). GSLIB: Geostatistical Software Library and User's Guide (2nd ed.). Oxford University Press. (The gamv program in GSLIB implements directional variograms with both an angular tolerance $\Delta\phi$ and an explicit bandwidth cap. The §II.2 documentation lays out the directional-cone formula used widely in industry practice.)
• 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 rotation-matrix parameterisation that §3.4 will pick up. The modern comprehensive reference.)
• Pyrcz, M.J., Deutsch, C.V. (2014). Geostatistical Reservoir Modeling (2nd ed.). Oxford University Press. (Reservoir-engineering treatment of anisotropy that emphasises the connection between the variogram ellipse and depositional process in fluvial, deltaic, and shoreface settings. The anisotropy-vs-process examples in this section's "connection to geology" subsection are drawn from this reference.)
• Wackernagel, H. (2003). Multivariate Geostatistics: An Introduction with Applications (3rd ed.). Springer. (Treats geometric and zonal anisotropy and the decomposition of zonal structures into nested components. Useful complement to Goovaerts and Chilès & Delfiner on why zonal anisotropy is usually a non-stationarity diagnostic.)