Field-level applications involving reservoir flow one-way-coupled with geomechanics—such as CO \(_2\) or hydrogen storage, and hydraulic fracturing—require accurate calculation of a self-equilibrating initial stress field, which requires balancing inputs such as pore pressure, intact rock density, and in-situ stresses. Moreover, we must honor specific measurements, which inevitably complicate the computation of this stress field due to inherent uncertainties. To address this compatibility problem, we minimize a well-defined objective function to ensure accurate and timely results. We employ a superposition principle and an inversion approach to calculate a consistent stress field that fits the available measurements and satisfies equilibrium conditions. To enforce compatibility, we minimize an objective function that allows us to adjust the input degrees of freedom (DOF) represented by piecewise linear polygonal functions of depth, thereby mimicking well-log data. We utilize a local Gauss-Newton (GN) optimizer that halves the number of iterations compared to analogous gradient-based optimizers. Additionally, we interpolate the resulting predictions with Duchon splines (a type of radial basis function) to estimate the sensitivity matrix, rather than using the traditional QR decomposition technique. We present extensive numerical results, including analytical nonlinear least-squares problems for benchmarking and calibration, a 2D plane-strain elasticity case, and full 3D synthetic field problems involving both homogeneous and heterogeneous poroelastic scenarios. These are one-way coupled with our in-house reservoir simulator. To assess our implementation and verify its consistency, we create inverse problems with known solutions. We employ this approach to tackle large synthetic reservoirs, demonstrating its scalability on parallel computers. We provide tabulated results, including the number of iterations, objective function values, and floating-point operations (FLOPs), to evaluate the performance of the inversion process. The results indicate that our approach produces consistent stress fields that conform to the predefined measurements while maintaining minimal computational cost in terms of iterations and simulation jobs. The included numerical results suggest that approximating the sensitivity matrix of the optimization process using Duchon splines outperforms the QR method, encouraging its use in inverse nonlinear least-squares approaches for the geomechanics initialization problem. We demonstrate scalability to large problems with hundreds of unknowns and robustness to noisy observations, suggesting the approach’s suitability for realistic field inversion problems.
