Graph theory for analyzing pair-wise data: application to geophysical model parameters estimated from interferometric synthetic aperture radar data at Okmok volcano,Alaska

Graph theory is useful for analyzing time-dependent model parameters estimated from interferometric synthetic aperture radar (InSAR) data in the temporal domain. Plotting acquisition dates (epochs) as vertices and pair-wise interferometric combinations as edges defines an incidence graph. The edge-vertex incidence matrix and the normalized edge Laplacian matrix are factors in the covariance matrix for the pair-wise data. Using empirical measures of residual scatter in the pair-wise observations, we estimate the relative variance at each epoch by inverting the covariance of the pair-wise data. We evaluate the rank deficiency of the corresponding least-squares problem via the edge-vertex incidence matrix. We implement our method in a MATLAB software package called GraphTreeTA available on GitHub (https://github.com/feigl/gipht). We apply temporal adjustment to the data set described in Lu et al. (Geophys Res Solid Earth 110, 2005) at Okmok volcano, Alaska, which erupted most recently in 1997 and 2008. The data set contains 44 differential volumetric changes and uncertainties estimated from interferograms between 1997 and 2004. Estimates show that approximately half of the magma volume lost during the 1997 eruption was recovered by the summer of 2003. Between June 2002 and September 2003, the estimated rate of volumetric increase is \((6.2 \, \pm \, 0.6) \times 10^6~\mathrm{m}^3/\mathrm{year} \). Our preferred model provides a reasonable fit that is compatible with viscoelastic relaxation in the five years following the 1997 eruption. Although we demonstrate the approach using volumetric rates of change, our formulation in terms of incidence graphs applies to any quantity derived from pair-wise differences, such as range change, range gradient, or atmospheric delay.     

Parameter estimation in nonlinear environmental problems

Popular parameter estimation methods, including least squares, maximum likelihood, and maximum a posteriori (MAP), solve an optimization problem to obtain a central value (or best estimate) followed by an approximate evaluation of the spread (or covariance matrix). A different approach is the Monte Carlo (MC) method, and particularly Markov chain Monte Carlo (MCMC) methods, which allow sampling from the posterior distribution of the parameters. Though available for years, MC methods have only recently drawn wide attention as practical ways for solving challenging high-dimensional parameter estimation problems. They have a broader scope of applications than conventional methods and can be used to derive the full posterior pdf but can be computationally very intensive. This paper compares a number of different methods and presents improvements using as case study a nonlinear DNAPL source dissolution and solute transport model. This depth-integrated semi-analytical model approximates dissolution from the DNAPL source zone using nonlinear empirical equations with partially known parameters. It then calculates the DNAPL plume concentration in the aquifer by solving the advection-dispersion equation with a flux boundary. The comparison is among the classical MAP and some versions of computer-intensive Monte Carlo methods, including the Metropolis–Hastings (MH) method and the adaptive direction sampling (ADS) method.     

On finite volume method implementation of poro‐elasto‐plasticity soil model

Accurate prediction of the interactions between the nonlinear soil skeleton and the pore fluid under loading plays a vital role in many geotechnical applications. It is therefore important to develop a numerical method that can effectively capture this nonlinear soil‐pore fluid coupling effect. This paper presents the implementation of a new finite volume method code of poro‐elasto‐plasticity soil model. The model is formulated on the basis of Biot's consolidation theory and combined with a perfect plasticity Mohr‐Coulomb constitutive relation. The governing equation system is discretized in a segregated manner, namely, those conventional linear and uncoupled terms are treated implicitly, while those nonlinear and coupled terms are treated explicitly by using any available values from previous time or iteration step. The implicit–explicit discretization leads to a linearized and decoupled algebraic system, which is solved using the fixed‐point iteration method. Upon the convergence of the iterative method, fully nonlinear coupled solutions are obtained. Also explored in this paper is the special way of treating traction boundary in finite volume method compared with FEM. Finally, three numerical test cases are simulated to verify the implementation procedure. It is shown in the simulation results that the implemented solver is capable of and efficient at predicting reasonable soil responses with pore pressure coupling under different loading situations. Copyright © 2015 John Wiley & Sons, Ltd.     

Introducing Transport “Surprises” in the Classroom: The Visible Fracture

Heterogeneity in aquifer properties, and the influence of transport processes other than advection and dispersion, often produce transport “surprises” in that measurements become difficult to reconcile with predictions from the traditional advection-dispersion equation (ADE) that students are introduced to early in their hydrogeology training. Students should be aware of and prepared for the reasons why the ADE (as commonly applied) may not always “work” in the sense of generating valid predictions. Though the predictive limitations of the ADE have been frequently discussed in the hydrogeologic literature, our experience is that students are not appropriately skeptical of transport predictions. For this reason, we believe it is imperative that future generations of hydrogeologists are introduced to transport surprises early in their formative education. We present a laboratory activity, centered around a “Visible Fracture,” which was presented in a laboratory class period of 75 min. The Visible Fracture consists of two sheets of Plexiglas surrounding a heterogenous fluid-filled aperture. Heated fluid is injected into the fracture, and thermochromic liquid crystal (TLC) materials are used to visualize the temperature both in the flowing fluid and in the surrounding Plexiglas “host rock.” Visualization of the plume shows the complex shapes that can be produced due to macroscopic heterogeneity. Tracer particles within the fracture allow students to examine heterogeneous local advective velocities, and to observe retardation of the fluid temperature plume. Student self-reported knowledge surveys indicate greater conceptual understanding of transport non-idealities after experiencing this activity.     

A Potential-Based Inversion of Unconfined Steady-State Hydraulic Tomography

The importance of estimating spatially variable aquifer parameters such as transmissivity is widely recognized for studies in resource evaluation and contaminant transport. A useful approach for mapping such parameters is inverse modeling of data from series of pumping tests, that is, via hydraulic tomography. This inversion of field hydraulic tomographic data requires development of numerical forward models that can accurately represent test conditions while maintaining computational efficiency. One issue this presents is specification of boundary and initial conditions, whose location, type, and value may be poorly constrained. To circumvent this issue when modeling unconfined steady-state pumping tests, we present a strategy that analyzes field data using a potential difference method and that uses dipole pumping tests as the aquifer stimulation. By using our potential difference approach, which is similar to modeling drawdown in confined settings, we remove the need for specifying poorly known boundary condition values and natural source/sink terms within the problem domain. Dipole pumping tests are complementary to this strategy in that they can be more realistically modeled than single-well tests due to their conservative nature, quick achievement of steady state, and the insensitivity of near-field response to far-field boundary conditions. After developing the mathematical theory, our approach is first validated through a synthetic example. We then apply our method to the inversion of data from a field campaign at the Boise Hydrogeophysical Research Site. Results from inversion of nine pumping tests show expected geologic features, and uncertainty bounds indicate that hydraulic conductivity is well constrained within the central site area.     

Hydraulic Tomography: Continuity and Discontinuity of High‐K and Low‐K Zones

Hydraulic tomography is an emerging field and modeling method that provides a continuous hydraulic conductivity (K) distribution for an investigated region. Characterization approaches that rely on interpolation between one‐dimensional (1D) profiles have limited ability to accurately identify high‐K channels, juxtapositions of lenses with high K contrast, and breaches in layers or channels between such profiles. However, locating these features is especially important for groundwater flow and transport modeling, and for design and operation of in situ remediation in complex hydrogeologic environments. We use transient hydraulic tomography to estimate 3D K in a volume of 15‐m diameter by 20‐m saturated thickness in a highly heterogeneous unconfined alluvial (clay to sand‐and‐gravel) aquifer with a K range of approximately seven orders of magnitude at an active industrial site in Assemini, Sardinia, Italy. A modified Levenberg‐Marquardt algorithm was used for geostatistical inversion to deal with the nonlinear nature of the highly heterogeneous system. The imaging results are validated with pumping tests not used in the tomographic inversion. These tests were conducted from three of five clusters of continuous multichannel tubing (CMTs) installed for observation in the tomographic testing. Locations of high‐K continuity and discontinuity, juxtaposition of very high‐K and very low‐K lenses, and low‐K “plugs” are evident in regions of the investigated volume where they likely would not have been identified with interpolation from 1D profiles at the positions of the pumping well and five CMT clusters. Quality assessment methods identified a suspect high‐K feature between the tested volume and a lateral boundary of the model.     

Fitting Data Under Omnidirectional Noise: A Probabilistic Method for Inferring Petrophysical and Hydrologic Relations

While geophysicists recognize that all measurements of media properties are subject to noise, methods for fitting petrophysical relations commonly employ regression-based formulas, which assume no error in one of the properties being related (the “independent" variable). To derive a more rigorous method for fitting such relations, we take a probabilistic viewpoint on the problem of fitting petrophysical relations to sample correlation data. Unlike prior approaches, we take into account the fact that noise is present in both measured properties. Under basic assumptions, we derive a new objective function for such problems which is proportional to the data likelihood and which can be used for both model parameter optimization and model selection. We present several numerical experiments outlining the utility of our method, and compare results of our method against results of other commonly used methods, such as kriging, regression, and total distance minimization. The results of our applications using the maximum likelihood technique appear visually accurate, and we also provide quantitative comparisons of performance that suggest the method produces more desirable results.     

Stochastic Cost Optimization of Multistrategy DNAPL Site Remediation

This paper investigates numerical optimization of dense nonaqueous phase liquid (DNAPL) site remediation design considering effects of prediction and measurement uncertainty. Results are presented for a hypothetical problem involving remediation using thermal source reduction (TSR) and bioremediation with electron donor (ED) injection. Pump-and-treat is utilized as a backup measure if compliance criteria are not met. Remediation system design variables are optimized to minimize expected net present value (ENPV) cost. Adaptive criteria are assumed for real-time control of TSR and ED duration. Source zone dissolved concentration data enabled more reliable and lower cost operation of TSR than soil concentration data, but using both soil and dissolved data improved results sufficiently to more than offset the additional cost. Decisions to terminate remediation and monitoring or to initiate pump-and-treat are complicated by measurement noise. Simultaneous optimization of monitoring frequency, averaging period, and lookback periods to confirm decisions, in addition to remediation design variables, reduced ENPV cost. Results indicate that remediation design under conditions of uncertainty is affected by subtle interactions and tradeoffs between design variables, compliance rules, site characteristics, and uncertainty in model predictions and monitoring data. Optimized designs yielded cost savings of up to approximately 50% compared with a nonoptimized design based on common engineering practices. Significant improvements in accuracy and reductions in cost were achieved by recalibrating the model to data collected during remediation and re-optimizing design variables. Repeating this process periodically is advisable to minimize total costs and maximize reliability.