Majken C. Looms

SIPPI: A Matlab toolbox for sampling the solution to inverse problems with complex prior information

From a probabilistic point-of-view, the solution to an inverse problem can be seen as a combination of independent states of information quantified by probability density functions. Typically, these states of information are provided by a set of observed data and some a priori information on the solution. The combined states of information (i.e. the solution to the inverse problem) is a probability density function typically referred to as the a posteriori probability density function. We present a generic toolbox for Matlab and Gnu Octave called SIPPI that implements a number of methods for solving such probabilistically formulated inverse problems by sampling the a posteriori probability density function. In order to describe the a priori probability density function, we consider both simple Gaussian models and more complex (and realistic) a priori models based on higher order statistics. These a priori models can be used with both linear and non-linear inverse problems. For linear inverse Gaussian problems we make use of least-squares and kriging-based methods to describe the a posteriori probability density function directly. For general non-linear (i.e. non-Gaussian) inverse problems, we make use of the extended Metropolis algorithm to sample the a posteriori probability density function. Together with the extended Metropolis algorithm, we use sequential Gibbs sampling that allow computationally efficient sampling of complex a priori models. The toolbox can be applied to any inverse problem as long as a way of solving the forward problem is provided. Here we demonstrate the methods and algorithms available in SIPPI. An application of SIPPI, to a tomographic cross borehole inverse problems, is presented in a second part of this paper.

Monitoring of saline tracer movement with vertically distributed self-potential measurements at the HOBE agricultural test site, Voulund, Denmark

The self-potential (SP) method is sensitive to water fluxes in saturated and partially saturated porous media, such as those associated with rainwater infiltration and groundwater recharge. We present a field-based study at the Voulund agricultural test site, Denmark, that is, to the best of our knowledge, the first to focus on the vertical self-potential distribution prior to and during a saline tracer test. A coupled hydrogeophysical modeling framework is used to simulate the SP response to precipitation and saline tracer infiltration. A layered hydrological model is first obtained by inverting dielectric and matric potential data. The resulting model that compares favorably with electrical resistance tomography models is subsequently used to predict the SP response. The electrokinetic contribution (caused by water fluxes in a charged porous soil) is modeled by an effective excess charge approach that considers both water saturation and pore water salinity. Our results suggest that the effective excess charge evolution prior to the tracer injection is better described by a recent flux-averaged model based on soil water retention functions than by a previously proposed volume-averaging model. This is the first time that raw (i.e., without post-processing or data-correction) vertically distributed SP measurements have been explained by a physically based model. The electrokinetic contribution cannot alone reproduce the experimental SP data during the tracer test and an electro-diffusive contribution (caused by concentration gradients) is needed. The predicted amplitude of this contribution is too small to perfectly explain the data, but the shape is in accordance with the field data. This discrepancy is attributed to imperfect descriptions of electro-diffusive phenomena in partially saturated soils, unaccounted soil heterogeneity, and discrepancies between the measured and predicted electrical conductivities in the tracer infiltration area. This study opens the way for detailed long-term field-based investigations of the SP method in vadose zone hydrology.

SIPPI: A Matlab toolbox for sampling the solution to inverse problems with complex prior information: Part 2—Application to crosshole GPR tomography

Abstract We present an application of the SIPPI Matlab toolbox, to obtain a sample from the a posteriori probability density function for the classical tomographic inversion problem. We consider a number of different forward models, linear and non-linear, such as ray based forward models that rely on the high frequency approximation of the wave-equation and ‘fat’ ray based forward models relying on finite frequency theory. In order to sample the a posteriori probability density function we make use of both least squares based inversion, for linear Gaussian inverse problems, and the extended Metropolis sampler, for non-linear non-Gaussian inverse problems. To illustrate the applicability of the SIPPI toolbox to a tomographic field data set we use a cross-borehole traveltime data set from Arrenaes, Denmark. Both the computer code and the data are released in the public domain using open source and open data licenses. The code has been developed to facilitate inversion of 2D and 3D travel time tomographic data using a wide range of possible a priori models and choices of forward models.

Geostatistical inference using crosshole ground-penetrating radar

High-resolution tomographic images obtained from crosshole geophysical measurements have the potential to provide valu- able information about the geostatistical properties of unsaturat- ed-zonehydrologic-statevariablessuchasmoisturecontent.Un- der drained or quasi-steady-state conditions, the moisture con- tentwillreflectthevariationofthephysicalpropertiesofthesub- surface, which determine the flow patterns in the unsaturated zone.Deterministicleast-squaresinversionofcrossholeground- penetrating-radarGPRtraveltimesresultinsmooth,minimum- variance estimates of the subsurface radar wave velocity struc- ture,whichmaydiminishtheutilityoftheseimagesforgeostatis- tical inference. We have used a linearized stochastic inversion technique to infer the geostatistical properties of the subsurface radar wave velocity distribution using crosshole GPR travel- times directly. Expanding on a previous study, we have deter- minedthatitispossibletoobtainestimatesofglobalvarianceand mean velocity values of the subsurface as well as the correlation lengths describing the subsurface velocity structures. Accurate estimation of the global variance is crucial if stochastic realiza- tions of the subsurface are used to evaluate the uncertainty of the inversion estimate. We have explored the full potential of the geostatisticalinferencemethodusingseveralsyntheticmodelsof varying correlation structures and have tested the influence of different assumptions concerning the choice of covariance func- tionanddatanoiselevel.Inaddition,wehavetestedthemethod- ology on traveltime data collected at a field site in Denmark. There, inferred correlation structures indicate that structural dif- ferences exist between two areas located approximately 10 m apart,anobservationconfirmedbyaGPRreflectionprofile.Fur- thermore, the inferred values of the subsurface global variance and the mean velocity have been corroborated with moisture- content measurements, obtained gravimetrically from samples collectedatthefieldsite.

Modeling cosmic ray neutron field measurements

The cosmic-ray neutron method was developed for intermediate-scale soil moisture detection, but may potentially be used for other hydrological applications. The neutron signal of different hydrogen pools is poorly understood and separating them is difficult based on neutron measurements alone. Including neutron transport modeling may accommodate this shortcoming. However, measured and modeled neutrons are not directly comparable. Neither the scale nor energy ranges are equivalent, and the exact neutron energy sensitivity of the detectors is unknown. Here, a methodology to enable comparability of the measured and modeled neutrons is presented. The usual cosmic-ray soil moisture detector measures moderated neutrons by means of a proportional counter surrounded by plastic, making it sensitive to epithermal neutrons. However, that configuration allows for some thermal neutrons to be measured. The thermal contribution can be removed by surrounding the plastic with a layer of cadmium, which absorbs neutrons with energies below 0.5 eV. Likewise, cadmium-shielding of a bare detector allows for estimating the epithermal contribution. First, the cadmium difference method is used to determine the fraction of thermal and epithermal neutrons measured by the bare and plastic-shielded detectors, respectively. The cadmium difference method results in linear correction models for measurements by the two detectors, and has the greatest impact on the neutron intensity measured by the moderated detector at the ground surface. Next, conversion factors are obtained relating measured and modeled neutron intensities. Finally, the methodology is tested by modeling the neutron profiles at an agricultural field site and satisfactory agreement to measurements is found. This article is protected by copyright. All rights reserved.

Monitoring unsaturated flow and transport using cross-borehole geophysical methods

Recent research has shown that cross-borehole georadar and electrical resistivity tomography (ERT) can provide data on soil moisture content and conductivity variations in the vadose zone at a more appropriate spatial scale than traditional techniques. A field site has been established in Denmark on a 20-30 m layer of unsaturated melt water sand and gravel deposits. Two identical field setups have been established each having four ERT and four georadar boreholes. The boreholes are drilled to a depth of 12 m, and form a cross consisting of two lines. Along each line the outer two boreholes (7 m apart) are equipped with ERT instrumented PVC-tubes (electrodes every 50 cm) while the inner two boreholes (5 m apart) have access tubes for georadar. Two different tracer infiltration experiments have been performed; (i) natural infiltration and (ii) forced infiltration experiments. In the natural infiltration tracer experiment, 200 L tracer was distributed on a 10 m by 10 m area, simulating a 2 mm rain event. In the forced infiltration tracer experiment, tracer was applied during the first 4 hours of the experiment at a constant rate of 250 l/hr (5.1 mm/hr) followed by irrigation of clean water at the same rate to accelerate flow during the successive 10 days. When using the two geophysical monitoring techniques simultaneously, it is possible to monitor variations in fluid conductivity resulting from tracer infiltration by combining the water content images from cross-borehole georadar with the bulk conductivity images achieved from the cross-borehole ERT. In both experiments, water content and conductivity data were collected prior to and after the applied tracer. The two sets of experiments provide data describing two vastly different flow conditions. In the natural infiltration experiment the upper boundary is determined by the naturally occurring precipitation and actual evapotranspiration. In this case little temporal variation in soil moisture content is observed and the tracer migration mainly takes place in autumn and winter where a surplus net-precipitation exists. In contrast, flow and transport velocities are much higher in the forced infiltration experiment and the upper boundary condition for this experiment is much more well-defined. Both experiments are interpreted by numerical hydraulic models which allow estimation of large-scale unsaturated hydraulic and transport parameters by appropriate inverse modeling.