An algebraic formulation of geophysical inverse problems

Many geophysical inverse problems derive from governing partial differential equations with unknown coefficients. Alternatively, inverse problems often arise from integral equations associated with a Green's function solution to a governing differential equation. In their discrete form such equations reduce to systems of polynomial equations, known as algebraic equations. Using techniques from computational algebra one can address questions of the existence of solutions to such equations as well as the uniqueness of the solutions. The techniques are enumerative and exhaustive, requiring a finite number of computer operations. For example, calculating a bound to the total number of solutions reduces to computing the dimension of a linear vector space. The solution set itself may be constructed through the solution of an eigenvalue problem. The techniques are applied to a set of synthetic magnetotelluric values generated by conductivity variations within a layer. We find that the estimation of the conductivity and the electric field in the subsurface, based upon single-frequency magnetotelluric field values, is equivalent to a linear inverse problem. The techniques are also illustrated by an application to a magnetotelluric data set gathered at Battle Mountain, Nevada. Surface observations of the electric ( E _y ) and magnetic ( H _x ) fields are used to construct a model of subsurface electrical structure. Using techniques for algebraic equations it is shown that solutions exist, and that the set of solutions is finite. The total number of solutions is bounded above at 134 217 728. A numerical solution of the algebraic equations generates a conductivity structure in accordance with the current geological model for the area.     

Numerical inversion of deformation caused by pressure sources: application to Mount Etna (Italy)

The interpretation of geodetic data in volcanic areas is usually based on analytical deformation models. Although numerical finite element (FE) modelling allows realistic features such as topography and crustal heterogeneities to be included, the technique is not computationally convenient for solving inverse problems using classical methods. In this paper, we develop a general tool to perform inversions of geodetic data by means of 3-D FE models. The forward model is a library of numerical displacement solutions, where each entry of the library is the surface displacement due to a single stress component applied to an element of the grid. The final solution is a weighted combination of the six stress components applied to a single element-source. The pre-computed forward models are implemented in a global search algorithm, followed by an appraisal of the sampled solutions. After providing extended testing, we apply the method to model the 1993–1997 inflation phase at Mt Etna, documented by GPS and EDM measurements. We consider four different forward libraries, computed in models characterized by homogeneous/heterogeneous medium and flat/topographic free surface. Our results suggest that the elastic heterogeneities of the medium can significantly alter the position of the inferred source, while the topography has minor effect.     

Dynamic non-planar crack rupture by a finite volume method

Modelling dynamic rupture for complex geometrical fault structures is performed through a finite volume method. After transformations for building up the partial differential system following explicit conservative law, we design an unstructured bi-dimensional time-domain numerical formulation of the crack problem. As a result, arbitrary non-planar faults can be explicitly represented without extra computational cost. On these complex surfaces, boundary conditions are set on stress fluxes and not on stress values. Prescribed rupture velocity gives accurate solutions with respect to analytical ones depending on the mesh refinement, while solutions for spontaneous propagation are analysed through numerical means. An example of non-planar spontaneous fault growth in heterogeneous media demonstrates the good behaviour of the proposed algorithm as well as specific difficulties of such numerical modelling.     

Multi-scale Numerical Simulation and 3D Modeling for Deep Mineral Exploration in the Jiaojia Gold District,China

In this study, deposit- and district-scale three-dimensional (3D) fault-and-intrusion structure models were constructed, based on which a numerical simulation was implemented in the Jiaojia gold district, China. The numerical simulation of the models shows the basic metallogenic path and trap of the gold deposits using mineral system theory. The objective of this study was to delineate the uncertainty of the geometry or buffer zones of the ore-forming and ore-controlling fault-and-intrusion domains in 3D environment representing the exploration criteria extraction and the gold potential targeting in the study area. The fast Lagrangian analysis of continua in three dimensions was used as the platform to define the stress deformation fracture ore storage and the hydrothermal seepage channel zone based on the gold deposit features and metallogenic model in the study area. The validity of the numerical simulation was verified by comparing it with robust 3D geological models of the large Xincheng gold deposit. The potential targeting zones are analyzed for uncertainty and then evaluated by Boolean operation in a 3D geological model using the computer-aided design platform. The research results are summarized as follows. (1) In the pre-mineralization period, the Jiaodong fault’s left lateral movement created the Jiaojia network faults and formed a fracture zone with NW- to NNW-trending dips of 20° to 40°. (2) During the mineralization period, hydrothermal flow was associated with the intrusion geometry and features. However, it was constrained by the Jiaojia fault, which blocked the vadose flow into the upper wall rock and made the hydrothermal route close to the fault in the footwall fracture zones. (3) Three gold potential targets were identified by the numerical simulation results in the study area: the NW-trending Sizhuang gold deposit, the NW-trending zone of Jiaojia gold deposit, and the NE-trending zone of the Xincheng gold deposit. (4) The numerical simulation results show the fault-and-intrusion metallogenic domain and the hydrothermal alteration zones, which reflect the main ore-controlling and ore-forming factors of mineralization. The information obtained through the numerical simulation discussed here can be used to define exploration criteria in the study area.

     

Stabilization of Bottom Hole Temperature With Finite Circulation Time and Fluid Flow

As an alternative to finite element or finite difference modelling, analytical solutions are derived by the method of Laplace transformation and numerical results obtained for several models of the bottom hole temperature stabilization. Included in the models are the finite circulation time, the thermal property contrast between the borehole mud and the surrounding formation, and the presence of radial or lateral fluid flow in the formation, all of which are found to have significant effects on the dissipation of the thermal disturbance induced by drilling.
The mud circulation is considered to have the effect of either maintaining the borehole mud at a constant temperature or supplying a constant amount of heat per unit length per unit time to the borehole. For small circulation times, the former reduces to the 'zero circulation' model in which the mud circulation creates an instantaneous temperature anomaly at the hole bottom; for small borehole radii, the latter reduces to the line source model and the traditional 'Homer plot'.
For typical drilling operations in which the bottom hole temperatures are measured several hours to several tens of hours after the hole is shut in, the new models generally predict higher equilibrium formation temperatures than does the Horner plot. However, predictions from the various models converge if the BHTs are taken after the hole has been shut in for a period which is greater than about five times the circulation period.     

Ground deformation in a viscoelastic medium composed of a layer overlying a half-space: a comparison between point and extended sources

We obtain and compare analytical and numerical solutions for ground displacement caused by an overpressurized magma chamber placed in a linear viscoelastic medium composed of a layer over a half-space. Different parameters such as size, depth and shape of the chamber, crustal rheology and topography are considered and discussed. Numerical solutions for an axisymmetric extended source are computed using a finite element method (FEM). Analytical solutions for a point source are obtained using the dislocation theory and the propagator matrix technique. In both cases, the elastic solutions are used together with the correspondence principle of linear viscoelasticity to obtain the solution in the Laplace transform domain. Viscoelastic solutions in the time domain are derived inverting the Laplace transform using the Prony series method. The differences between the results allow us to constrain the applicability of the point source and the flat surface hypothesis, which are usually implicitly assumed when analytical solutions are derived. The effect of the topography is also considered. The results obtained show that neglecting the topographic effects may, in some cases, introduce an error greater than that implicit in the point-source hypothesis. Therefore, for an adequate modelling and interpretation of the time-dependent displacements, topography must be considered.     

Reduced complexity strategies for modelling urban floodplain inundation

Significant advances in flood inundation modelling have been made in the last decade through the use of a new generation of 2D hydraulic numerical models. These offer the potential to predict the local pattern and timing of flood depth and velocity, enabling informed flood risk zoning and improved emergency planning. With the availability of high resolution DEMs derived from airborne lidar, these models can theoretically now be routinely parameterized to represent considerable topographic complexity, even in urban areas where the potential exists to represent flows at the scale of individual buildings. Currently, however, computational constraints on conventional finite element and volume codes typically require model discretization at scales well below those achievable with lidar and are thus unable to make optimal use of this emerging data stream.In this paper we review two strategies that attempt to address this mismatch between model and data resolution in an effort to improve urban flood forecasts. The first of these strives for a solution by simplifying the mathematical formulation of the numerical model by using a computationally efficient 2D raster storage cell approach coupled to a 1D channel model. This parsimonious model structure enables simulations over large model domains offering the opportunity to employ a topographic discretization strategy which explicitly represents the built environment. The second approach seeks to further reduce the computational overhead of this raster method by employing a subgrid parameterization to represent the effect of buildings and micro-relief on flow pathways and floodplain storage. This multi-scale methodology enables highly efficient model applications at coarse spatial resolutions while retaining information about the complex geometry of the built environment.These two strategies are evaluated through numerical experiments designed to reconstruct a flood in the small town of Linton in southern England, which occurred in response to a 1 in 250 year rainfall event in October 2001. Results from both approaches are encouraging, with the spatial pattern of inundation and flood wave propagation matching observations well. Both show significant advantages over a coarse resolution model without subgrid parameterisation, particularly in terms of their ability to reproduce both hydrograph and inundation depth measurements simultaneously, without need for recalibration. The subgrid parameterization is shown to achieve this without contributing significant computational complexity and reduces model run-times by an order of magnitude.     

B-polarization induction in two thin half-sheets coupled to the mantle by a conducting crust—II. Solution by the mode-matching method, the fields above the conductors, and example results

The mode-matching method is used to obtain an exact analytical solution to the problem of B -polarization induction in two adjacent thin half-sheets, lying on a conducting layer that is terminated by a perfect conductor at finite depth. These components of the model represent, respectively, the Earth's conducting surface layers, crust, and mantle. In dimensionless variables, the model has three independent parameters, these being the two thin-sheet conductances and the layer thickness. The mode-matching solution obtained in this paper is shown to be identical lo that derived via the Wiener-Hopf method in a companion paper (Dawson 1996), and so provides additional verification of that solution. As was shown in the companion paper, the solution for the present model contains, as special limiting cases, those for three models considered earlier by various authors. The second part of the present paper addresses the solutions for the electric fields in the non-conducting half-space above the conductors, which represents the atmosphere. In the final part, sample numerical calculations are presented to illustrate the solution.     

A hybrid solution for wave propagation problems in regular media with bounded irregular inclusions

Summary. This paper describes a method that enables one to calculate the effects of localized heterogeneities on the wavefield in an otherwise regular medium. It does so by connecting a finite element solution for a heterogeneous inclusion to any type of solution for the regular medium, e.g. a reflectivity solution for a layered medium or an analytical solution for a simple half-space. Once the Green's functions for the regular medium are determined, the method reduces to a coupled set of algebraic equations for the wavefield, with the incident field and/or body forces as known variables. An efficient numerical scheme is derived for the solution of these equations.     

Effect of the overconsolidation ratio of soils in surface settlements due to tunneling

Construction of urban tunnels requires the control of surface subsidence to minimize any disturbance to nearby buildings and services. Past study of surface subsidence has been limited to mainly empiri...     

Mechanical and numerical models for sea ice dynamics on small-meso scale

On small-meso scale,the sea ice dynamic characteristics are quite different from that on large scale.To model the sea ice dynamics on small-meso scale,a new elastic-viscous-plastic(EVP) constitutive model and a hybrid Lagrangian-Eulerian (HLE) numerical method are developed based on continuum theory.While a modified discrete element model(DEM) is introduced to model the ice cover at discrete state.With the EVP constitutive model,the numerical simulation for ice ridging in an idealized rectangular basin is carried out and the results are comparable with the analytical solution of jam theory.Adopting the HLE numerical model,the sea ice dynamic process is simulated in a vortex wind field.The furthering application of DEM is discussed in details for modeling the discrete distribution of sea ice.With this study ,the mechanical and numerical models for sea ice dynamics can be improved with high precision and computational efficiency.     

Adaptive unstructured grid finite element simulation of two-dimensional magnetotelluric fields for arbitrary surface and seafloor topography

We present an adaptive unstructured triangular grid finite element approach for effectively simulating plane-wave diffusive electromagnetic fields in 2-D conductivity structures.
The most striking advantage of irregular grids is their potential to incorporate arbitrary geometries including surface and seafloor topography. Adaptive mesh refinement strategies using an a posteriori error estimator yield most efficient numerical solutions since meshes are only refined where required.
We demonstrate the robustness of this approach by comparison with analytical solutions and previously published numerical simulations. Maximum errors may systematically be reduced to, for example, 0.8 per cent for the apparent resistivity and 0.2° in the phase.
An additional accuracy study of the thickness of the air layer in E-polarization suggests to keep a minimum thickness depending on lateral conductivity contrasts within the earth.
Furthermore, we point out the new quality and flexibility of our simulation technique by addressing two marine magnetotelluric applications. In the first case, we discuss topographic effects associated with a synthetic sinusoidal sea bottom model and in the second case, we show a close-to-reality scenario using real bathymetry data from the East Pacific Rise at 17°S.     

Gravitational viscoelastic relaxation of eccentrically nested spheres

We present a semi-analytical solution to the 2-D forward modelling of viscoelastic relaxation in a heterogeneous model consisting of eccentrically nested spheres. Several numerical methods for 2-D and 3-D viscoelastic relaxation modelling have been applied recently, including finite-element and spectral-finite-difference schemes. The present semi-analytical approach provides a model response against which more general numerical algorithms can be validated. The eccentrically nested sphere solution has been tested by comparing it with the analytical solutions for viscoelastic relaxation in a homogeneous sphere and in two concentrically nested spheres, and good agreement was obtained.     

Data model for system conceptualization in groundwater studies

Development of a conceptualization of a hydrogeologic system serves as the basis of groundwater modeling. While existing groundwater data models are designed to store groundwater system information, none is designed to capture its conceptual view. This study addresses this need by presenting a new object-oriented Conceptualization Groundwater Data Model that represents a groundwater system as a series of aquifer layers with defined aquifer properties and water boundary conditions. A case study is presented that develops the conceptual view of the groundwater system beneath Konza Prairie. This single conceptualization is used to support groundwater models across existing technologies of finite difference, finite element, and analytical element methods. While the models each employ different mathematics, data input files, and formats, all models are founded on the same conceptualization process that is represented using this new data model. The case study illustrates the data model's promise as an effective mechanism for groundwater system conceptualization and data storage, and utility for various groundwater computational models. This conceptualization of a groundwater data model suggests a new focus on incorporating system conceptualization into data model design.     

Salt tectonics driven by differential sediment loading: stability analysis and finite-element experiments

At many continental margins, differential sediment loading on an underlying salt layer drives salt deformation and has a significant impact on the structural evolution of the basin. We use 2‐D finite‐element modelling to investigate systems in which a linear viscous salt layer underlies a frictional‐plastic overburden of laterally varying thickness. In these systems, differential pressure induces the flow of viscous salt, and the overburden experiences updip deviatoric tension and downdip compression. A thin‐sheet analytical stability criterion for the system is derived and is used to predict conditions under which the sedimentary overburden will be unstable and fail, and to estimate the initial velocities of the system. The analytical predictions are in acceptable agreement with initial velocity patterns of the numerical models. In addition to initial stability analyses, the numerical model is used to investigate the subsequent finite deformation. As the systems evolve, overburden extension and salt diapirism occur in the landward section and contractional structures develop in the seaward section. The system evolution depends on the relative widths of the salt basin and the length scale of the overburden thickness variation. In narrow salt basins, overburden deformation is localised and characterised by high strain rates, which cause the system to reach a gravitational equilibrium and salt movement to cease earlier than for wide salt basins. Sedimentation enhances salt evacuation by maintaining a differential pressure in the salt. Continued sedimentary filling of landward extensional basins suppresses landward salt diapirism. Sediment progradation leads to seaward propagation of the landward extensional structures and depocentres. At slow sediment progradation rates, the viscous flow can be faster than the sediment progradation, leading to efficient salt evacuation and salt weld formation beneath the landward section. Fast sediment progradation suppresses the viscous flow, leaving salt pillows beneath the prograding wedge.     

Semi-analytical solution for viscous Stokes flow in two eccentrically nested spheres

Applying the infinite Prandtl number approximation, a semi-analytical solution for computing 2-D axisymmetric viscous Stokes flow in a model consisting of two eccentrically nested spheres of different viscosities is derived. Since numerical codes based on spectral or finite techniques for modelling mantle flow in a spherical geometry in the presence of lateral viscosity variation are becoming more and more popular, reliable examples for testing and validating such codes are extremely useful. The eccentrically nested sphere solution was used to test a numerical algorithm based on a mixed spherical-harmonic finite-element formulation of the Stokes problem, and good agreement was obtained.     

Application of the boundary element method to a terrestrial heat flow problem

b
In this paper, we apply the boundary element method (BEM) to the 2-D steady state heat flow problem of what would be the perturbation to the regional temperature gradient, and hence heat flow density, determined from temperatures measured in a borehole that passes close to, but does not penetrate, a body of anomalous thermal properties. This type of problem with an infinity boundary is particularly well suited to the BEM.
The results have been compared with those obtained from analytical solutions for bodies of simple shape; it is found that for the worst case of a close approach to a boundary of small radius of curvature, a numerical modelling error of less than 1 per cent can still be obtained provided the length of each element is less than the shortest distance between the calculation point and the object.     

Inversion for multiple parameter classes

Many geophysical data, such as the frequencies of the free oscillations of the Earth, depend on more than one type of model parameter. For inverse problems depending on multiple parameter classes, an iterative solution procedure is introduced in which each parameter class can be treated in the same way. This approach has considerable advantages where a large number of parameters are employed, but can still be useful for smaller systems.
  The iteration by parameter class commences by solving for the direct dependence on a particular parameter class, and at subsequent iterations the cross-dependences between classes are introduced. The update affects only the right-hand side of the equations, and, because the same sets of equations have to be solved at each iteration, an efficient computational implementation can be made. The largest set of equations that has to be solved at a time corresponds to the number of variables in an individual parameter class rather than the full set of parameters, which confers substantial computational benefits for very large problems.     

Different discretization method used in coupled water and heat transport mode for soil under freezing conditions

A coupled water and heat transport mode is established based on the Richards equation to study water flow and heat transport in soil during freezing process. Both the finite difference and finite element method are used in the discretization, respectively. Two different computer programs are written and used to simulate an indoor unidirectional frozen test. The freezing depth, freezing rate and temperature variation are compared among lab tests, finite difference calculation simulation and finite element calculation simulation. Result shows that: the finite difference method has a better performance in freezing depth simulation while the finite element method has a better performance in numerical stability in one-dimensional freezing simulation.     

Comparison of methods to model the gravitational gradients from topographic data bases

A number of methods have been developed over the last few decades to model the gravitational gradients using digital elevation data. All methods are based on second-order derivatives of the Newtonian mass integral for the gravitational potential. Foremost are algorithms that divide the topographic masses into prisms or more general polyhedra and sum the corresponding gradient contributions. Other methods are designed for computational speed and make use of the fast Fourier transform (FFT), require a regular rectangular grid of data, and yield gradients on the entire grid, but only at constant altitude. We add to these the ordinary numerical integration (in horizontal coordinates) of the gradient integrals. In total we compare two prism, two FFT and two ordinary numerical integration methods using 1" elevation data in two topographic regimes (rough and moderate terrain). Prism methods depend on the type of finite elements that are generated with the elevation data; in particular, alternative triangulations can yield significant differences in the gradients (up to tens of Eötvös). The FFT methods depend on a series development of the topographic heights, requiring terms up to 14th order in rough terrain; and, one popular method has significant bias errors (e.g. 13 Eötvös in the vertical–vertical gradient) embedded in its practical realization. The straightforward numerical integrations, whether on a rectangular or triangulated grid, yield sub-Eötvös differences in the gradients when compared to the other methods (except near the edges of the integration area) and they are as efficient computationally as the finite element methods.