The stochastic geostatistical inversion approach is widely used in subsurface inverse problems to estimate unknown parameter fields and corresponding uncertainty from noisy observations. However, the approach requires a large number of forward model runs to determine the Jacobian or sensitivity matrix, thus the computational and storage costs become prohibitive when the number of unknowns, m, and the number of observations, n increase. To overcome this challenge in large‐scale geostatistical inversion, the Principal Component Geostatistical Approach (PCGA) has recently been developed as a “matrix‐free” geostatistical inversion strategy that avoids the direct evaluation of the Jacobian matrix through the principal components (low‐rank approximation) of the prior covariance and the drift matrix with a finite difference approximation. As a result, the proposed method requires about K runs of the forward problem in each iteration independently of m and n, where K is the number of principal components and can be much less than m and n for large‐scale inverse problems. Furthermore, the PCGA is easily adaptable to different forward simulation models and various data types for which the adjoint‐state method may not be implemented suitably. In this paper, we apply the PCGA to representative subsurface inverse problems to illustrate its efficiency and scalability. The low‐rank approximation of the large‐dimensional dense prior covariance matrix is computed through a randomized eigen decomposition. A hydraulic tomography problem in which the number of observations is typically large is investigated first to validate the accuracy of the PCGA compared with the conventional geostatistical approach. Then the method is applied to a large‐scale hydraulic tomography with 3 million unknowns and it is shown that underlying subsurface structures are characterized successfully through an inversion that involves an affordable number of forward simulation runs. Lastly, we present a joint inversion of head and tracer test data using MODFLOW and MT3DMS as coupled black‐box forward simulation solvers. These applications demonstrate the advantages of the PCGA, i.e., the scalability to high‐dimensional inverse problems and the ability to utilize multiple forward models as black boxes.

The quasi‐linear geostatistical approach is for weakly nonlinear underdetermined inverse problems, such as Hydraulic Tomography and Electrical Resistivity Tomography. It provides best estimates as well as measures for uncertainty quantification. However, for its textbook implementation, the approach involves iterations, to reach an optimum, and requires the determination of the Jacobian matrix, i.e., the derivative of the observation function with respect to the unknown. Although there are elegant methods for the determination of the Jacobian, the cost is high when the number of unknowns, m, and the number of observations, n, is high. It is also wasteful to compute the Jacobian for points away from the optimum. Irrespective of the issue of computing derivatives, the computational cost of implementing the method is generally of the order of m²n, though there are methods to reduce the computational cost. In this work, we present an implementation that utilizes a matrix free in terms of the Jacobian matrix Gauss‐Newton method and improves the scalability of the geostatistical inverse problem. For each iteration, it is required to perform K runs of the forward problem, where K is not just much smaller than m but can be smaller that n. The computational and storage cost of implementation of the inverse procedure scales roughly linearly with m instead of m² as in the textbook approach. For problems of very large m, this implementation constitutes a dramatic reduction in computational cost compared to the textbook approach. Results illustrate the validity of the approach and provide insight in the conditions under which this method perform best.

The Surface Water and Ocean Topography (SWOT) satellite has the potential to transform global hydrologic science by offering simultaneous and synoptic estimates of river discharge and other hydraulic variables. Discharge is estimated from SWOT observations of water surface elevation, width, and slope. A first assessment using just the highest quality SWOT measurements, over the first 15 months (March 2023–July 2024) of the mission evaluated at 65 gauged reaches shows results consistent with pre‐launch expectations. SWOT estimates track discharge dynamics without relying on any gauge information: median correlation is 0.73, with a correlation interquartile range of 0.51–0.89. SWOT estimates capture discharge magnitude correctly in some cases but are biased (median bias is 50%) in others. There are already a total of 11,274 ungauged global locations with highest quality SWOT measurements where SWOT discharge is expected to accurately track discharge variations: this value will increase as SWOT data record length grows, algorithms are refined and SWOT measurements are reprocessed. This first look indicates that SWOT discharge is performing as expected for SWOT data that achieve performance requirements, providing observed information on discharge variations in ungauged basins globally.