Curvelet-Regularized SPDE Inversion on Piecewise-Planar Fractures with Trace-Graph Coupling

Technical Summary: Curvelet-Regularized SPDE Inversion on Piecewise-Planar Fractures with Trace-Graph Coupling

Problem Statement
Classical geostatistical methods, such as kriging and Gaussian Random Field (GRF) models, typically assume correlation is a function of Euclidean distance. This assumption fails in fractured media, where transport and flow propagate preferentially along connected fracture networks rather than through intact rock. Consequently, points close in Euclidean space may be weakly coupled, while distant points connected by a fracture corridor may be strongly coupled. Furthermore, fracture-controlled fields often exhibit anisotropic, curvilinear, and multiscale features (e.g., ridges, fronts, channels) that are poorly represented by isotropic smoothness priors alone. The paper addresses the challenge of reconstructing a latent scalar field (e.g., log-transmissivity, geochemical grade) on fracture surfaces given sparse point measurements and 3D fracture trace polylines.

Methodology
The proposed framework operates in two main stages: geometric inference and sparse-to-dense reconstruction.

1. Geometric Inference (3D Traces to Plane Charts):

   • Plane Fitting: 3D trace polylines are fitted to planes using PCA/SVD. Degenerate traces (near-collinear) are flagged and excluded or assigned post-hoc.
   • Clustering: Plane hypotheses are clustered using a combined metric of normal vector deviation and offset misfit (DBSCAN-style) to identify dominant fracture planes.
   • Local Coordinates: For each dominant plane, a local orthonormal basis $(u, v)$ is established. 3D samples are projected onto these local charts, and a 2D grid is constructed to discretize the latent field $z_p(u, v)$.

2. Reconstruction Models:

   • Baseline (Gaussian/SPDE): A standard Gaussian MAP estimator is formulated using a Sparse Precision Matrix derived from a Stochastic Partial Differential Equation (SPDE). The precision operator is defined as $Q = (\kappa^2 I + L_{uv})^2$, where $L_{uv}$ is the discrete 2D Laplacian. This provides a smooth, isotropic baseline.
   • Curvelet Regularization: To capture directional, fracture-aligned structures, the model adds an $\ell_1$ penalty on undecimated discrete curvelet (UDCT) coefficients. The objective function minimizes the data misfit, the SPDE quadratic prior, and the curvelet sparsity term:
     $$\min_z \frac{1}{2}\|W^{1/2}(Hz - y)\|_2^2 + \frac{1}{2}z^\top Q z + \lambda \|Cz\|_1$$
     where $C$ is the curvelet transform.
   • Trace-Graph Coupling (Version 2 Extension): To enforce "along-fracture" distance, a graph Laplacian $L_G$ is constructed on the trace network. A coupling term $z^\top P^\top L_G P z$ is added, where $P$ maps grid values to graph nodes. This penalizes differences between connected vertices along the fracture network, regardless of Euclidean distance.

3. Optimization:

   • The non-smooth curvelet term is solved using the Alternating Direction Method of Multipliers (ADMM).
   • The algorithm splits the problem into a sparse linear solve (for the quadratic SPDE block) and a coefficient-wise shrinkage step (for the curvelet block).
   • Convergence is guaranteed under standard convexity assumptions, with stopping criteria based on primal and dual residual norms.

Key Contributions

• Framework for Fractured Media: A reproducible pipeline that infers piecewise-planar fracture supports from 3D traces and reconstructs property fields from sparse data.
• Directional Sparsity Prior: The formal integration of curvelet transforms into SPDE-based geostatistics to specifically target curvilinear singularities and anisotropic features inherent to fractures.
• Connectivity-Driven Regularization: The definition of an "along-fracture" distance metric via graph Laplacians, allowing the model to enforce smoothness along fracture networks rather than in Euclidean space.
• Modular Architecture: The paper distinguishes between a "Version 1" (per-plane inference with Gaussian and curvelet priors) and a "Version 2" extension (explicit graph coupling and intersection constraints), providing a clear roadmap for increasing physical complexity.
• Reproducible Benchmarks: The authors provide synthetic benchmarks (S1: curvilinear anisotropy; S2: graph-driven connectivity) and a full suite of ablation studies to isolate the contributions of Euclidean smoothness, directional sparsity, and connectivity priors.

Results
Numerical experiments on synthetic data demonstrate:

• Directional Sparsity (S1): The curvelet-regularized model significantly outperforms ordinary kriging and the pure SPDE baseline in reconstructing ridge-like anomalies aligned with fracture traces, achieving lower RMSE and higher $R^2$ on held-out samples.
• Connectivity (S2): In scenarios where ground truth is driven by graph connectivity rather than Euclidean distance, the addition of the graph prior (Version 2) improves reconstruction accuracy. The full model (SPDE + Curvelets + Graph) yields the best performance, effectively capturing the non-Euclidean correlation structure.
• Efficiency: The ADMM solver converges efficiently, leveraging sparse linear algebra for the quadratic block and fast transforms ($O(M \log M)$) for the curvelet block.

Significance and Scope
The paper positions itself as a methods framework rather than a complete physics solver. Its primary significance lies in providing a statistically controlled interpolation/inversion layer that can be coupled to downstream flow or transport models.

• Modest Claims: The authors explicitly state they do not claim "physics completeness." The framework is designed to reconstruct fields conditioned on fracture geometry, not to solve the governing PDEs of fluid flow directly.
• Foundation for Future Work: The "Version 2" extension and the roadmap (Section 14) outline future integration with nonlinear physics and PDE-constrained inversions.
• Validation: The paper emphasizes a rigorous validation protocol, including sensitivity studies on hyperparameters ($\kappa, \lambda, \rho$), grid resolution, and sampling density, establishing a baseline for network-aware geostatistics in fractured media.