Sequential Gaussian Simulation (SGS)
Learning objectives
- State the SGS algorithm: random visit order + local kriging + Normal sampling
- Recognise SGS as the scalable workhorse for large grids (O(N k³) cost)
- Distinguish the search neighbourhood k from the global covariance
- Document the trade-off: scalability vs approximation of the global covariance
- Apply SGS to generate ensemble realisations for resource estimation
Cholesky/LU (§7.2) gives exact unconditional simulation but doesn't scale beyond ~1000 grid nodes. SEQUENTIAL GAUSSIAN SIMULATION (SGS; Isaaks 1990, Deutsch & Journel 1998) achieves SCALABILITY via a SEQUENTIAL approach: visit grid nodes one at a time, simulate each from local kriging conditional on data + already-simulated nodes. The result is approximate but scales to millions of nodes.
The SGS algorithm
- INITIALISE: known data set = hard data.
- RANDOMISE: random visit order of unsampled grid nodes.
- FOR each node in visit order:
- a. Find k NEAREST known points (data + already-simulated) within search ellipsoid.
- b. KRIGE: compute mean prediction and kriging variance using these k points.
- c. SAMPLE: draw value from N(prediction, kriging variance).
- d. ADD: include simulated value in the known data set.
- END FOR. The resulting grid is one SGS realisation.
Different random seeds → different realisations. Each realisation honours data + approximates the variogram. The ensemble characterises spatial uncertainty.
Why SGS is scalable
For each of N nodes, SGS does ONE local kriging with k points (typically k = 20-50). Cost per node: O(k³) for matrix factorisation. Total cost: O(N × k³). With k constant, this is LINEAR in N — vastly cheaper than Cholesky's O(N³).
| Method | Time | Memory | Scalable? |
|---|---|---|---|
| Cholesky | O(N³) | O(N²) | N ≤ ~1000 |
| SGS | O(N k³) | O(N + k²) | N up to 10^7+ |
The scalability-vs-exactness trade-off
SGS's local-neighbourhood approach is APPROXIMATE: it considers only k nearest known points, not the full data set. The global covariance is reproduced APPROXIMATELY, not exactly. Larger k = better covariance reproduction but slower simulation. Standard practice: k = 20-50.
Modern best practice: VALIDATE the simulation by comparing the realised variogram (from many simulated grids) to the assumed variogram model. Should match within Monte Carlo error.
Why random visit order matters
If you visit nodes in REGULAR ORDER (e.g., left-to-right), early-visited nodes have FEWER previously-simulated neighbours; later nodes have MORE. This induces SYSTEMATIC BIAS in the realised covariance.
RANDOM visit order distributes this asymmetry uniformly across the grid. Modern SGS implementations always randomise.
Modern variants of SGS
- p-field simulation: post-process to enforce histogram matching.
- Direct simulation: bypass the Gaussian assumption via empirical conditional CDFs.
- Multi-Gaussian SGS for facies: simulate Gaussian field, threshold to get facies.
- Sequential indicator simulation (SISIM): indicator-variogram-based facies sim (§8.3).
- Multipoint statistics (MPS): training-image-based simulation (§9.2).
For continuous variables (porosity, grade), SGS is the modern default.
Try it
- Default range = 2, k = 10. The widget runs full SGS: visits 44 unsampled nodes (50 grid - 6 data) in random order, krigs and samples each. Final blue curve passes EXACTLY through green data points.
- Re-simulate (different seed). Same data + variogram, but different realisation. The data are honoured identically; the unsampled regions vary — exactly the uncertainty quantification we want.
- Crank range to 5.0 (long correlation). Realisations are smoother. Each kriged point is strongly influenced by neighbours.
- Drop range to 0.5. Realisations are much rougher; nodes are nearly independent. SGS still completes quickly.
- Drop k to 3. Sparse search neighbourhood. Each kriging uses fewer points; the variogram is approximated more poorly. Realisations may show artifacts. Increase k to 25-30 for better reproduction; standard production uses k = 20-50.
A reservoir model simulates permeability on a 100x100x10 grid (100,000 nodes). With SGS (k=30) it takes ~10 minutes. The same problem with Cholesky would take ~10^15 operations — infeasible. Why?
What you now know
SGS is the scalable workhorse for conditional simulation. Random visit + local kriging + Normal sampling per node. O(N k³) — linear in N. Approximate global covariance but accurate for k ≥ 20. Modern best practice: validate by re-computing realised variogram. §7.4 next: conditioning on hard data — the algorithmic mechanism that ensures realisations honour observations exactly.
References
- Isaaks, E.H. (1990). "The application of Monte Carlo methods to the analysis of spatially correlated data." PhD thesis, Stanford University. (Original SGS.)
- Gómez-Hernández, J.J., Journel, A.G. (1993). "Joint sequential simulation of multi-Gaussian fields." Geostatistics Tróia '92, 85–94.
- Deutsch, C.V., Journel, A.G. (1998). GSLIB, 2nd ed. Oxford. (GSLIB SGSIM program.)
- Goovaerts, P. (1997). Geostatistics for Natural Resources Evaluation, Chapter 8. Oxford.
- Pyrcz, M.J., Deutsch, C.V. (2014). Geostatistical Reservoir Modeling, 2nd ed. Oxford.