Autoregressive Physics-Informed Neural Operator Surrogates and Joint 4-D Seismic-Production Bayesian History Matching

1.1 Engineering application

• Production data from three wells cannot resolve inter-well permeability connectivity in the 100,188-cell Tau reservoir, motivating spatially distributed 4-D seismic observations.The synthetic model has about 2×10^5 unknown property values but only O(10^2) production measurements.Source: §1.1 · p. 2

1.2.1 From response surfaces to neural operators

• Classical emulators and CNN/PINN models are described as design-ensemble-limited or resolution-locked, whereas neural operators learn solution operators with discretisation-invariant, quasi-linear-cost evaluations.The paper places its FNO within the operator-learning family and uses a discretised PDE residual as a physics penalty.Source: §1.2.1 · pp. 2–3

1.2.3 Bayesian inverse problems and ensemble data assimilation

• The paper adopts a function-space Bayesian posterior and ensemble Kalman methods for derivative-free inversion, while identifying transport or score-based methods as a route beyond the Gaussianity assumption.Its DCT and VCAE charts address finite-ensemble rank deficiency and geological structure in latent parameter spaces.Source: §1.2.3 · p. 4

1.2.4 4-D seismic history matching

• Time-lapse seismic can provide areal information absent from sparse well data, but its value is case dependent and signals may be weak or ambiguous.Here, a closed-form petro-elastic transform of the surrogate state synthesises seismic inside the rollout, with expected information gain motivated by a linearised Fisher-information argument.Source: §1.2.4 · pp. 4–5

1.3 Research gap and approach

• The proposed framework combines a closed-loop autoregressive PINO surrogate, in-rollout seismic and Peaceman production observables, and a two-stage Bayesian inversion.The inversion uses aREKI followed by ES-MDA; DCT or VCAE charts reduce the petrophysical parameter space.Source: §1.3 · p. 5, Figure 1 · p. 6, §1.4 · p. 6
• The paper positions its contribution as an integrated framework for three-phase hydrocarbon history matching, not as the first combination of surrogates and seismic data.Its stated combination includes full-history autoregressive rollouts, consistent production and seismic observables, Peaceman well responses, and ensemble updates tied to a function-space posterior.Source: §1.3 · p. 5

1.2.2 Autoregressive surrogates and exposure bias

• Teacher-forced models face exposure bias in closed-loop rollouts; the paper trains its operator autoregressively with truncated backpropagation, a discrete residual, and noise injection.The paper presents these choices as extensions of prior sequential-PINO work, alongside well-response heads and joint seismic-production inversion.Source: §1.2.2 · pp. 3–4, §1.3 · p. 5

1.4 Contributions and outline

• The stated contributions span autoregressive PINO theory, Peaceman well modelling, Bayesian inversion, Gassmann seismic synthesis, reduced petrophysical charts, and a two-stage aREKI schedule.The paper also reports a Tau-model inverse-modelling demonstration with a blind forecast window and explicitly stated scope.Source: §1.4 · pp. 5–6

2.1.1 Continuum model

• The black-oil model describes water, oil, and gas flow with solution gas, using pressure and saturation fields as its primary state variables.Under the small-capillary approximation, the system reduces to an elliptic global pressure equation coupled to advection-dominated saturation transport.Source: §2.1.1 · p. 6

2.1.2 Constitutive relations

• Fault transmissibility multipliers scale inter-cell fluxes, with zero transmissibility representing a sealing fault; fault fields enter the surrogate as static input channels.The constitutive setup includes Stone II relative permeabilities and pressure-dependent fluid properties.Source: §2.1.2 · p. 7, §2.1.3 · p. 7

2.1.3 PVT, compressibility and the equilibration timescale

• Pressure equilibration is estimated at about 115 days, motivating report intervals that grow from one day to 183 days after pressure equilibration.The stated diffusivity hierarchy contrasts globally relaxing pressure with locally advecting saturation.Source: §2.1.3 · p. 7

2.1.4 Finite-volume discretisation and the discrete residual

• The finite-volume residual is reused as a PINO penalty without inversion, requiring one O(N) assembly per step and no symbolic PDE differentiation.The reference simulator drives the residual below 10^-6 with Newton–Krylov.Source: §2.1.4 · p. 7, §2.1.8 · p. 8

2.1.5 Fractional flow and front-speed: the source of saturation error

• Saturation-front position error can grow linearly in time after a pressure-gradient error, motivating greater water-field loss weight and increment penalties.The text connects this front-direction amplification to mildly declining long-horizon oil R2 while pressure R2 remains stable.Source: §2.1.5 · p. 7

2.1.6 The composite loss in detail

• The composite objective combines absolute and temporal-increment penalties, with larger Peaceman-head weight because well observations are sparse but informative.Pressure and water weights are 2, oil is 1, and Peaceman is 10; increment weights are one-tenth their corresponding absolute weights.Source: §2.1.6 · pp. 7–8

2.1.7 FNO complexity and parameter budget

• The four-layer FNO uses 32 latent channels and 16 retained Fourier modes per dimension, yielding about 8.4×10^6 parameters per spectral layer.The paper contrasts its O(N log N) spectral forward pass with the O(N^1.5–N^2) implicit solve and attributes the operator choice to resolution invariance.Source: §2.1.7 · p. 8, §2.1.8 · p. 8

2.1.9 Data extraction and the training tensor

• The training data comprise 200 simulation realisations over 20 report steps on the 46×99×22 grid, and the surrogate maps augmented controls, state, completion geometry, and time to the next state.For this undersaturated model, free-gas saturation is identically zero and the surrogate advances pressure, water saturation, and oil saturation.Source: §2.1.9 · p. 8, §2.1.10 · p. 8

2.2.1 The well-block pressure problem

• Peaceman modelling relates grid-block pressure to flowing bottom-hole pressure through an equivalent drainage radius for an unresolved well-block pressure profile.The paper uses the isotropic radius re = 0.2√(ΔxΔy) and well radius rw = 0.1 m.Source: §2.2 · p. 8, §2.2.1 · p. 9

2.2.2 Layer well index and phase rates

• Layer well indices and phase rates are aggregated across completions, and the assimilated producer-pressure observable is a transmissibility-weighted drainage-pressure proxy rather than flowing bottom-hole pressure.Producer wells have 20 completion intervals and the injector has 15; assimilated responses include rates, pressures, and cumulative water cut.Source: §2.2.2 · p. 9

2.2.3 Why a learned Peaceman head

• Five dedicated one-dimensional FNO heads predict Peaceman well responses because the saturation-dependent rate map is strongly nonlinear near water breakthrough.The per-well features couple global pressure and saturation predictions to multilayer completion physics, which the comparison study reports as decisive for rate fidelity.Source: §2.2.3 · pp. 9–10

2.3.2 The eleven-head composite operator

• The composite operator has 11 heads: three dynamic-state heads, three initial-state heads, and five Peaceman well-response heads.The initial-state gradient path is detached before the dynamic rollout graph, decoupling the two training paths.Source: Table 2 · p. 10, §2.3.2 · p. 11

2.3.1 The Fourier neural operator

• The FNO retains 16 Fourier modes per dimension; vertical resolution is complete, but lateral fault discontinuities and saturation fronts remain under-resolved.A higher-capacity per-cell Conv1d decoder sharpens local detail, while the projection-rate argument does not apply across saturation fronts.Source: §2.3.1 · pp. 10–11

2.3.3 Exposure bias: a total-variation analysis

• The exposure analysis uses Wasserstein-1 distance because arbitrarily small differences between deterministic flows can still yield mutually singular laws.When the rollout Lipschitz constant L exceeds one, the error bound can grow geometrically, weakening one-step empirical risk as a proxy for deployment loss.Source: §2.3.3 · p. 11

2.3.3 Exposure bias: a total-variation analysis; 2.3.4 Closed-loop (autoregressive) training

• Closed-loop training feeds surrogate predictions back into the rollout, addressing the mismatch between training on true states and deployment on predicted states.The proposed stability advantage is conditional: a lower closed-loop Lipschitz constant is empirical, and the trajectory must remain within the specified tube.Source: §2.3.3 · p. 11, §2.3.4 · p. 12

2.3.5 The PINO physics penalty in detail

• The PINO penalty adds the discrete black-oil residual to data fitting, but contraction under a trained residual defect requires control of both the residual and its derivative.The paper describes the resulting contraction condition as structural rather than certified because the relevant defect quantities are not estimated.Source: §2.3.5 · pp. 12–13

2.3.6 Covariate-shift (noise) injection; 2.3.7 K-step truncated BPTT; 2.3.8 Training algorithm; 2.3.9 The truncated-BPTT computational graph

• Noise injection trains on perturbed feedback states, while K=10 truncated BPTT detaches the graph at window boundaries and bounds backward memory.Under the stated contraction assumptions, the truncation bias decays geometrically; for ρ ≲ 0.5 and K=10, its factor is about 2 × 10^-3.Source: §2.3.6 · p. 13, §2.3.7 · pp. 13–14, §2.3.9 · p. 14

2.3.11 Forward training results

• 24.90 final loss and 23.57 best validation followed 700 training epochs, with rapid early pressure and initial-state fitting, then saturation learning and Peaceman refinement.The loss components converge within approximately 200 epochs, while the reported three-phase analysis describes continued refinement through epoch 700.Source: §2.3.10 · p. 14, §2.3.11 · p. 15, Figure 3 · p. 15

2.4 Closed-loop inference

• At inference, each step uses the model’s previous output, with zero teacher forcing; Gaussian state perturbations provide a stochastic variant for uncertainty rollouts.Memory recycling enables 190 ensemble members on one 192 GB GPU.Source: §2.4 · p. 15
• Across the 20-step rollout, pressure differences stay below 2500 kPa and saturation errors remain localised at the injector front.The errors grow only mildly over the rollout.Source: §2.4 · p. 16
• For Tau_1/Tau_2, the surrogate’s well-response R2 is >0.9999 for WOPR and 0.98 for WBHP_PROD.Figure 5 also reports pressure/saturation R2 > 0.999 and high-fidelity reproduction of the Gassmann-derived seismic delta.Source: Figure 4 · p. 16, Figure 5 · p. 16

2.5.1 The Cluster-Classify-Regress decomposition

• CCR clusters near-well feature and response pairs, classifies inputs to a cluster, then applies that cluster’s local expert.For Tau, its experts are regularised, early-stopped XGBoost trees.Source: §2.5.1 · p. 18

2.5.2 Why the operator wins

• FNO WOPR R2 is > 99.99% versus CCR 18.7%, and FNO WWPR R2 is 82.8% versus CCR −72.4%.FNO inference takes 0.579 s on one GB200, versus 42 s for CCR on 144 CPU cores.Source: §2.5.2 · p. 18, Table 4 · p. 18
• CCR wins on injector WBHP, dominated by the local injectivity index; the operator wins where the global field matters.CCR’s 10 producer features cannot capture the global pressure field or track the saturation front across the rollout.Source: §2.5.2 · p. 18

3 The Bayesian Inverse Problem

• The inverse problem is constructed rigorously as the foundation for the ensemble methods that follow.Source: §3 · p. 18

3.1 Measure-theoretic foundations

• The unknown is the pair (logK, φ) on a function space; the grid discretisation has dimension 2N = 200,376, and separability supports the measure-theoretic ensemble representation.The framework is stated to survive mesh refinement.Source: §3.1 · p. 18, §3.1.1 · pp. 18–19

3.1.2 The Gaussian prior and its Cameron–Martin space

• A Gaussian prior requires a positive, self-adjoint, trace-class covariance; its Cameron–Martin norm supplies the prior-precision regulariser, while prior draws can be rougher than that space.The framework also allows a non-Gaussian geostatistical prior represented by samples, with a Gaussian reference in the VCAE latent space.Source: §3.1.2 · p. 19

3.1.3 Forward operator and the data model

• The forward map composes the rollout, Peaceman heads, and, for joint assimilation, the Gassmann seismic map; the data model uses observation-noise covariance Γ.The posterior assumptions require lower boundedness, local upper boundedness, and local Lipschitz properties of the potential in the unknown and data.Source: §3.1.3 · p. 19, §3.1.4 · p. 19

3.1.5 Posterior as a Radon–Nikodym density

• Under the potential assumptions, the posterior is defined as a Radon–Nikodym density with respect to the prior, avoiding a nonexistent infinite-dimensional Lebesgue density.Source: §3.1.5 · p. 20

3.1.6 Specialisation to a Lebesgue density

• On the finite grid, the posterior has a Lebesgue density; the function-space formulation instead uses the prior-relative density and supports mesh-independent ensemble solving.The grid MAP regulariser is the prior precision, consistent with the function-space Onsager–Machlup functional.Source: §3.1.6 · p. 20

3.1.7 Well-posedness in the Hellinger metric

• For bounded data, the posterior is locally Lipschitz in Hellinger distance as a function of the data, and its square-integrable moments are locally Lipschitz as well.Source: §3.1.7 · p. 21
• The surrogate posterior retains existence and well-posedness when its potential satisfies the stated assumptions, but closeness to the full-physics posterior is a separate requirement.The posterior’s local Hellinger stability also bounds errors in square-integrable posterior expectations when the forward-map error is controlled.Source: §3.1.7 · p. 21

3.1.8 Ensemble (particle) approximation

• Ensemble Kalman methods approximate the posterior with particles and an empirical Gaussian projection; with Ne ≪ d, rank deficiency motivates reparametrisation and localisation.The prior ensemble uses 190 samples from a pool of 1400 geostatistical realisations, which can preserve higher-order geological structure.Source: §3.1.8 · p. 21, §3.1.9 · p. 22

3.2.1 What 4-D seismic measures

• ∆Vp maps swept and depleted volume, providing spatial information about inter-well connectivity that point-well production data cannot see.The closed-form petro-elastic map transforms the surrogate state into time-lapse seismic velocity.Source: §3.2.1 · p. 22, §3.2 · p. 22

3.2.2 Gassmann fluid substitution

• The Gassmann channel derives saturated bulk modulus and velocity from pressure, oil saturation, porosity, dry-frame properties, and fluid properties.The dry frame is pressure-sensitive, while the shear modulus is fluid-independent under Gassmann’s assumption.Source: §3.2.2 · pp. 22–23
• The differentiable, backend-agnostic map is evaluated at baseline and monitor acquisition steps, with inactive cells replaced by the active-cell mean.This implementation supports placement inside an automatic-differentiation loop and dispatch across GPUs.Source: §3.2.2 · p. 23, §3.2.4 · p. 23

3.2.3 The 4-D observable and its compression

• The reported runs use the raw 4,554-value areal map; an alternative DCT encoding retains 5% of modes while preserving the innovation norm and reducing dimension by about 20×.The DCT truncation removes high-wavenumber modes while retaining the swept-volume signature.Source: §3.2.3 · p. 23

3.2.4 Noise balancing in the joint problem

• Optional seismic-noise rescaling balances modality leverage but departs from the stated likelihood, so the resulting posterior is not the exact Bayesian posterior for Γ.The paper reports a run with this option and notes that seismic information gain should be shown to persist without it.Source: §3.2.4 · p. 23

3.3.1 Why reduce dimension at all

• With d = 2N ≈ 2 × 10^5 parameters and 190 members, the sample covariance has rank at most 189, restricting updates to the ensemble-anomaly span.The section presents dimension reduction and covariance localisation as remedies; the reported run used grid-space inversion with sufficient rank.Source: §3.3.1 · p. 23

3.3.2 DCT parametrisation (linear chart)

• The DCT chart retains about 225 of 4,554 areal modes, reducing effective dimension to about 10^2 and yielding a full-rank reduced covariance.Its orthonormal linear basis captures smooth-field variance but cannot represent sharp multimodal facies boundaries without Gibbs ringing.Source: §3.3.2 · p. 24

3.3.3 VCAE parametrisation (nonlinear chart)

• The VCAE chart learns nonlinear latent coordinates from the prior ensemble and decodes them to geologically plausible fields, including connected geobodies and sharp facies contacts.Its latent prior is standard Gaussian, and the chart can represent genuinely multimodal posteriors that a single Gaussian or linear DCT basis cannot.Source: §3.3.3 · p. 24
• Both charts restrict EKI to a space the ensemble can sample, restoring full rank to the reduced cross-covariance.DCT is a fixed orthonormal frame; VCAE is a learned chart adapted to the prior manifold.Source: §3.3.3 · p. 24

3.4.1 Ensemble Kalman Inversion; 3.4.2 Continuous-time limit and the gradient-flow structure

• In the continuous-time limit, the EKI ensemble mean follows preconditioned gradient descent on the misfit in the covariance metric, confined to the ensemble-anomaly span.In the linear-Gaussian case, the flow converges to the MAP/Tikhonov minimiser and the particle spread reproduces posterior covariance.Source: §3.4.2 · pp. 24–25

3.4.3 Adaptive regularisation: aREKI

• aREKI adaptively chooses inflation toward the discrepancy principle’s target misfit, while limiting ensemble collapse and respecting the ES-MDA likelihood budget.At convergence, the stated calibration is E[2J] ≈ Nop, or one unit of χ^2 per observation.Source: §3.4.3 · p. 25

3.4.4 ES-MDA as exact tempered assimilation

• ES-MDA composes multiple tempered likelihood updates into the full likelihood in the linear-Gaussian case, using gentler steps to handle forward-map nonlinearity.The deterministic schedule αk = Na avoids the variance-criterion degeneracy that can affect aREKI on a preconditioned ensemble.Source: §3.4.4 · p. 25

3.4.5 Covariance localisation

• Covariance localisation tapers spurious long-range sample correlations; it was disabled for the reported run and is recommended for Ne < 100 or finer grids.The compactly supported taper zeroes correlations beyond 2r0.Source: §3.4.5 · p. 25

3.4.6 The two-stage production→joint schedule

• The two-stage schedule first conditions on production with aREKI, restores about 2% of ensemble variance, then assimilates production and seismic jointly with ES-MDA.The joint phase uses ES-MDA because the Phase-1-conditioned ensemble has collapsed misfit variance, where aREKI’s variance criterion degenerates.Source: §3.4.6 · p. 26

3.4.7 Posterior selection

• 0.876 (joint) versus 2.556 (production) on the normalized misfit scores selects P2, with a gap of 1.68.The posterior score uses prior-RMSE-normalized misfit in each modality, with P2 preferred in a tie.Source: §3.4.7 · p. 26

3.4.8 Non-monotone convergence and Jensen bias

• The ensemble-mean misfit can rise transiently under nonlinear EKI even as the best-member RMSE improves, so posterior statistics are extracted from the final ensemble.The paper attributes the transient behavior to an ensemble nonlinearity bias with indefinite sign.Source: §3.4.8 · p. 26

3.4.9 Why joint helps: a linearised information argument

• In the linearised argument, seismic information adds constraints across production’s null space, but the covariance-dominance inequality is not a proof for the nonlinear problem or seismic model error.Production sensitivities are concentrated in a low-dimensional subspace, while the areal seismic sensitivities include the swept inter-well region.Source: §3.4.9 · pp. 26–27

4.1.1 Workflow stages

• The ISIM workflow trains an autoregressive PINO surrogate, computes production and seismic observables, then applies aREKI to production followed by ES-MDA to both data types.The stages also select a posterior and extract percentile fields and best-fit members.Source: §4.1.1 · p. 27

4.1.2 PyTorch DistributedDataParallel

• The four-GPU forward job uses DDP gradient averaging, mixed precision, and K=10 truncated backpropagation; the inverse stage batches 190 surrogate rollouts.Near-ideal weak scaling is expected but was not measured; the reported 19.9 s/epoch is for four GPUs.Source: §4.1.3 · p. 28, §4.1.2 · p. 28, §4.1.4 · p. 28

4.2.1 Data split and run configuration

• The experiment uses 190 inversion members and withholds the final five of 20 report steps for blind forecasting, but its synthetic seismic data share the inversion’s petro-elastic model.Forward accuracy is assessed on only five held-out realizations, limiting the strength of generalization claims.Source: §4.2.1 · pp. 28–29

4.2.2 Joint vs. production-only: posterior comparison

• Joint assimilation recovers the high-permeability corridor and fault block, while its posterior water-rate ensemble is tighter than production-only assimilation.The joint posterior structure matches the true model; Figure 11 also shows tighter concentration over the forecast window.Source: §4.2.2 · p. 29, Figure 10 · p. 32, Figure 11 · p. 33

4.2.3 Producer bottom-hole pressure: prior-to-posterior collapse

• The joint posterior BHP ensemble tracks observations through the forecast window, but its P10–P90 band narrows within the 25% production-noise envelope.The paper cautions that a band tighter than the observation noise may indicate over-collapse and should be considered alongside spread diagnostics.Source: §4.2.3 · p. 29, Figure 12 · p. 33

4.2.4 Cost reduction

• The ensemble-mean RMSE drops by nearly two orders of magnitude, and MAP/best models lie at the low-RMSE tail across all four responses.The joint seismic constraint sharpens the water-rate distribution most strongly.Source: §4.2.4 · p. 35, Figure 14 · p. 35, §4.2.5 · p. 35

4.2.6 Quantitative dominance and summary

• The joint score is 0.876, below the prior-normalized reference of 1, while the production-only score is 2.556; P2 is selected.Table 5 compares production-only P1 with joint seismic–production P2; lower scores are better.Source: §4.2.6 · p. 36, Table 5 · p. 36

4.2.7 Posterior uncertainty quantification

• Posterior uncertainty is narrower over the seismically imaged corridor and wider at the aquifer periphery, which neither data modality observes.The reported deliverable is a posterior distribution summarized by permeability percentiles, not a single model.Source: §4.2.7 · p. 36

5.1 Limitations and scope

• Full-physics validation of selected posterior members is outstanding, so some apparent match may reflect surrogate-error structure.The study calls for running selected MAP and best members through the full-physics simulator against observations.Source: §5.1 · p. 37
• The Tau experiment is synthetic end to end: seismic data use the inversion petro-elastic code and omit model error, repeatability noise, and registration or depth-conversion error.Source: §5.1 · p. 37
• The two-stage schedule gives production data total weight about 2, tightening the P2 posterior in row(Jp); its absolute spread is not a calibrated credibility statement.The passage recommends a partial production budget in Phase 2, or treating Phase 1 as a preconditioner, so the full schedule has a single budget.Source: §5.1 · p. 37
• A diagonal seismic error covariance overstates information because the spatially correlated ΔVp map has effective degrees of freedom O(10^2), not 4554.Source: §5.1 · p. 37
• The operator has about 2 × 10^8 parameters, but training uses 190 realisations and the held-out test set only five, which cannot exclude memorisation.Source: §5.1 · p. 37
• The study does not report spread diagnostics despite a posterior predictive band narrower than observation noise, and the assumed Jacobian contraction lacks direct measurement and ablation.The Peaceman heads use a two-well feature layout, so transfer to another well count or grid requires retraining.Source: §5.1 · p. 37

5.2 Conclusions

• R^2 > 0.999 on black-oil state and well responses, with a 0.579 s rollout; joint matching sharpens the posterior and recovers inter-well connectivity on the synthetic experiment.The blind forecast window is tracked, but full-physics verification of the posterior remains outstanding.Source: §5.2 · pp. 37–38
• The closed-form, GPU-batched Gassmann channel is fused into the rollout at negligible marginal cost, making the surrogate a 4-D seismic forward model.Source: §5.2 · pp. 37–38

• 4-D means time-lapse seismic; MAP means maximum a posteriori; aREKI means adaptive regularised ensemble Kalman inversion.Source: Nomenclature · p. 38
• PINO means physics-informed neural operator, FNO Fourier neural operator, DCT discrete cosine transform, and VCAE variational convolutional auto-encoder.Source: Nomenclature · p. 38
• EKI means ensemble Kalman inversion, ES-MDA ensemble smoother with multiple data assimilation, and TBPTT truncated backpropagation through time.Source: Nomenclature · p. 38

• Tau model data cannot be released and the software is unavailable, but the specified workflow can be independently reimplemented with other reservoir models and simulators.Source: Data and code availability · p. 38