Advances in numerical algorithms and methods in computational geosciences with modeling characteristics of multiple physical and chemical processes

This paper aims to provide a brief introduction to recent advances in numerical algorithms and methods in the emerging computational geoscience filed with general simulation characteristics of modeling multiple chemical and physical processes that take place in ore-generating systems within the Earth's crust. Due to significant differences between Earth systems and engineering systems, the existing numerical algorithms and methods, which are designed for simulating realistic problems in the engineering fields, may not be straightforwardly used to simulate ore-generating problems without significant improvements. Thus, extensive and systematic studies have been conducted, in recent years, to develop new numerical algorithms and methods for simulating different aspects of ore-generating problems. Not only can the outcomes of these studies provide new simulation tools for better understanding the controlled dynamic mechanisms that take place in ore-generating systems, but also they have enriched the research contents of computational mechanics in the broad sense.

Computational modeling of free-surface slurry flow problems using particle simulation method

The particle simulation method is used to solve free-surface slurry flow problems that may be encountered in several scientific and engineering fields.The main idea behind the use of the particle simulation method is to treat granular or other materials as an assembly of many particles.Compared with the continuum-mechanics-based numerical methods such as the finite element and finite volume methods,the movement of each particle is accurately described in the particle simulation method so that the free surface of a slurry flow problem can be automatically obtained.The major advantage of using the particle simulation method is that only a simple numerical algorithm is needed to solve the governing equation of a particle simulation system.For the purpose of illustrating how to use the particle simulation method to solve free-surface flow problems,three examples involving slurry flow on three different types of river beds have been considered.The related particle simulation results obtained from these three examples have demonstrated that:1) The particle simulation method is a promising and useful method for solving free-surface flow problems encountered in both the scientific and engineering fields;2) The shape and irregular roughness of a river bed can have a significant effect on the free surface morphologies of slurry flow when it passes through the river bed.

Replacement of annular domain with trapezoidal domain in computational modeling of nonaqueous-phase-liquid dissolution-front propagation problems

In order to simulate the instability phenomenon of a nonaqueous phase liquid(NAPL) dissolution front in a computational model, the intrinsic characteristic length is commonly used to determine the length scale at which the instability of the NAPL dissolution front can be initiated. This will require a huge number of finite elements if a whole NAPL dissolution system is simulated in the computational model. Even though modern supercomputers might be used to tackle this kind of NAPL dissolution problem, it can become prohibitive for commonly-used personal computers to do so. The main purpose of this work is to investigate whether or not the whole NAPL dissolution system of an annular domain can be replaced by a trapezoidal domain, so as to greatly reduce the requirements for computer efforts. The related simulation results have demonstrated that when the NAPL dissolution system under consideration is in a subcritical state, if the dissolution pattern around the entrance of an annulus domain is of interest, then a trapezoidal domain cannot be used to replace an annular domain in the computational simulation of the NAPL dissolution system.However, if the dissolution pattern away from the vicinity of the entrance of an annulus domain is of interest, then a trapezoidal domain can be used to replace an annular domain in the computational simulation of the NAPL dissolution system. When the NAPL dissolution system under consideration is in a supercritical state, a trapezoidal domain cannot be used to replace an annular domain in the computational simulation of the NAPL dissolution system.

Computational simulation of convective flow in the Earth crust under consideration of dynamic crust-mantle interactions

The finite element method was used to solve fluid dynamic interaction problems between the crust and mantle of the Earth. To consider different mechanical behaviours, the lithosphere consisting of the crust and upper mantle was simulated as fluid-saturated porous rocks, while the upper aesthenospheric part of the mantle was simulated as viscous fluids. Since the whole lithosphere was computationally simulated, the dynamic interaction between the crust and the upper mantle was appropriately considered. In particular, the mixing of mantle fluids and crustal fluids was simulated in the corresponding computational model. The related computational simulation results from an example problem demonstrate that the mantle fluids can flow into the crust and mix with the crustal fluids due to the resulting convective flows in the crust-mantle system. Likewise, the crustal fluids can also flow into the upper mantle and mix with the mantle fluids. This kind of fluids mixing and exchange is very important to the better understanding of the governing processes that control the ore body formation and mineralization in the upper crust of the Earth.

Simulation of rock deformation and mechanical characteristics using clump parallel-bond models

To properly simulate hard rock with a high ratio of the uniaxial compressive strength to tensile strength(UCS/TS) and realistic strength-failure envelope,the rock deformation and mechanical characteristics were discussed in detail when the particle simulation method with the clump parallel-bond model(CPBM) was used to conduct a series of numerical experiments at the specimen scale.Meanwhile,the effects of the loading procedure and crack density on the mechanical behavior of a specimen,which was modeled by the particle simulation method with the CPBM,were investigated.The related numerical results have demonstrated that:1) The uniaxial compressive strength(UCS),tensile strength(TS) and elastic modulus are overestimated when the conventional loading procedure is used in the particle simulation method with the CPBM; 2) The elastic modulus,strength and UCS/TS decrease,while Poisson ratio increases with the increase of the crack density in the particle simulation method with the CPBM; 3) The particle simulation method with the CPBM can be used to reproduce a high value of UCS/TS(>10),as well as a high friction angle and reasonable cohesion strength; 4) As the confining pressure increases,both the peak strength of the simulated specimen and the number of microscopic cracks increase,but the ratio of tensile cracks number to shear cracks number decreases in the particle simulation method with the CPBM; 5) Compared with the conventional parallel-bond model,the CPBM can be used to reproduce more accurate results for simulating the rock deformation and mechanical characteristics.

Finite element modeling of pore-fluid flow in the Dachang ore district,Guangxi,China:Implications for hydrothermal mineralization

Convective heat transfer associated with the circulation of porefluid in porous rocks and fractures within the upper crust of the Earth is substantial when the temperature gradient is sufficiently high. In order to understand the process of Snpolymetallic mineralization in the Dachang ore district of Guangxi, a finite element method has been used in this study to simulate both pore-fluid flow and heat transfer in this district. On the basis of related geological, tectonic and geophysical constraints, a computational model was established. It enables a computational simulation and sensitivity analysis to be carried out for investigating ore-forming pore-fluid flow and other key factors that may affect hydrothermal ore genesis in the district. The related simulation results have indicated that： （1） permeable fault zones in the Dacbang ore district can serve as preferential pathways for pore-fluid flow on a regional-scale; and （2） the pore-fluid flow can affect the salinity distribution. This latter factor is part of the reason why Sn-polymetallic mineralization has taken place in this district.

Effective loading algorithm associated with explicit dynamic relaxation method for simulating static problems

Based on the fact that a static problem has an equivalent wave speed of infinity and a dynamic problem has a wave speed of finite value, an effective loading algorithm associated with the explicit dynamic relaxation method was presented to produce meaningful numerical solutions for static problems. The central part of the explicit dynamic relaxation method is to turn a time-independent static problem into an artificial time-dependent dynamic problem. The related numerical testing results demonstrate that: (1) the proposed effective loading algorithm is capable of enabling an applied load in a static problem to be propagated throughout the whole system within a given loading increment, so that the time-independent solution of the static problem can be obtained; (2) the proposed effective loading algorithm can be straightforwardly applied to the particle simulation method for solving a wide range of static problems.

Effects of porosity heterogeneity on chemical dissolution-front instability in fluid-saturated rocks

Homogeneity and heterogeneity are two totally different concepts in nature.At the particle length scale,rocks exhibit strong heterogeneity in their constituents and porosities.When the heterogeneity of porosity obeys the random uniform distribution,both the mean value and the variance of porosities in the heterogeneous porosity field can be used to reflect the overall heterogeneous characteristics of the porosity field.The main purpose of this work is to investigate the effects of porosity heterogeneity on chemical dissolution front instability in fluid-saturated rocks by the computational simulation method.The related computational simulation results have demonstrated that:1) since the propagation speed of a chemical dissolution front is inversely proportional to the difference between the final porosity and the mean value of porosities in the initial porosity field,an increase in the extent of the porosity heterogeneity can cause an increase in the mean value of porosities in the initial porosity field and an increase in the propagation speed of the chemical dissolution front.2) An increase in the variance of porosities in the initial porosity field can cause an increase in the instability probability of the chemical dissolution front in the fluid-saturated rock.3) The greater the mean value of porosities in the initial porosity field,the quicker the irregular morphology of the chemical dissolution front changes in the supercritical chemical dissolution systems.This means that the irregular morphology of a chemical dissolution front grows quicker in a porosity field of heterogeneity than it does in that of homogeneity when the chemical dissolution system is at a supercritical stage.

Effects of different numerical algorithms on simulation of chemical dissolution-front instability in fluid-saturated porous rocks

Many scientific and engineering problems need to use numerical methods and algorithms to obtain computational simulation results because analytical solutions are seldom available for them.The chemical dissolution-front instability problem in fluid-saturated porous rocks is no exception.Since this kind of instability problem has both the conventional(i.e.trivial)and the unconventional(i.e.nontrivial)solutions,it is necessary to examine the effects of different numerical algorithms,which are used to solve chemical dissolution-front instability problems in fluid-saturated porous rocks.Toward this goal,two different numerical algorithms associated with the commonly-used finite element method are considered in this paper.In the first numerical algorithm,the porosity,pore-fluid pressure and acid/solute concentration are selected as basic variables,while in the second numerical algorithm,the porosity,velocity of pore-fluid flow and acid/solute concentration are selected as basic variables.The particular attention is paid to the effects of these two numerical algorithms on the computational simulation results of unstable chemical dissolution-front propagation in fluid-saturated porous rocks.The related computational simulation results have demonstrated that:1)the first numerical algorithm associated with the porosity-pressure-concentration approach can realistically simulate the evolution processes of unstable chemical dissolution-front propagation in chemical dissolution systems.2)The second numerical algorithm associated with the porosity-velocity-concentration approach fails to simulate the evolution processes of unstable chemical dissolution-front propagation.3)The extra differential operation is the main source to result in the failure of the second numerical algorithm.

Finite element modeling of convective pore-fluid flow in fluid-saturated porous rocks within upper crust:An overview

Convective pore-fluid flow (CPFF) plays a critical role in generating mineral deposits and oil reservoirs within the deep Earth. Therefore, theoretical understanding and numerical modeling of the thermodynamic process that triggers and controls the CPFF are extremely important for the exploration of new mineral deposits and underground oil resources. From the viewpoint of science, the CPFF within the upper crust can be treated as a kind of thermodynamic instability problem of pore-fluid in fluid-saturated porous media. The key issue of dealing with this kind of problem is to assess whether a nonlinear thermodynamic system under consideration is supercritical. To overcome limitations of using theoretical analysis and experimental methods in dealing with the CPFF problems within the upper crust, finite element modeling has been broadly employed for solving this kind of problem over the past two decades. The main purpose of this paper is to overview recent developments and applications of finite element modeling associated with solving the CPFF problems in large length-scale geological systems of complicated geometries and complex material distributions. In particular, two kinds of commonly-used finite element modeling approaches, namely the steady-state and transient-state approaches, and their advantages/disadvantages are thoroughly presented and discussed.

Potential mechanisms of pore-fluid movement from continental lithospheric mantle into upper continental crust

Through integrating the state of the art scientific knowledge in different research fields, some potential mechanisms of large-scale movements of underground pore-fluids such as H2O and CO2 in the continental lithosphere were presented and discussed. The results show that the generation and propagation of porosity waves are important mechanisms to transport mass and heat fluxes from the continental lithospheric mantle into the lower continental crust; the generation and propagation of porosity waves, pore-fluid flow focusing through lower and middle crustal faults, advection of pore-fluids through the lower and middle crust, and whole-crust convection in some particular cases are important mechanisms to transport mass and heat fluxes from the lower into the upper continental crust; heat and mass transport through convective pore-fluid flow is the most effective mechanism of ore body formation and mineralization in hydrothermal systems; due to heat and mass exchange at the interface between the earth surface, hydrosphere and atmosphere, it is very important to consider the hydro-geological effect of the deep earth pore-fluids such as H2O and CO2 on the global warming and climate change in future investigations.

Existence of Weak Solutions for a Class of Quasilinear Parabolic Problems in Weighted Sobolev Space

In this paper, we investigate the existence and uniqueness of weak solutions for a new class of initial/boundary-value parabolic problems with nonlinear perturbation term in weighted Sobolev space. By building up the compact imbedding in weighted Sobolev space and extending Galerkin’s method to a new class of nonlinear problems, we drive out that there exists at least one weak solution of the nonlinear equations in the interval [0,T] for the fixed time T>0.

Numerical modeling of coupled fluid flow,heat transport and mechanical deformation:An example from the Chanziping ore district,South China

This paper presents numerical investigation on the ore-forming fluid migration driven by tectonic deformation and thermally-induced buoyancy force in the Chanziping ore district in South China.A series of numerical scenarios are considered to examine the effect of meteoric water precipitation, the dip angle of the faults,unconformity surface,and thermal input on the ore genesis.Our computations reveal that the downward basinal fluid flow driven by extensional stress mixes with the upward basal fluid driven by the thermal input from depth at the junction of two faults at a temperature of about 200℃,triggering the precipitation of the Chanziping uranium deposit.

CONVERGENCE OF THE CRANK-NICOLSON/NEWTON SCHEME FOR NONLINEAR PARABOLIC PROBLEM

In this paper, the Crank-Nicolson/Newton scheme for solving numerically second- order nonlinear parabolic problem is proposed. The standard Galerkin finite element method based on P2 conforming elements is used to the spatial discretization of the problem and the Crank-Nieolson/Newton scheme is applied to the time discretization of the resulted finite element equations. Moreover, assuming the appropriate regularity of the exact solution and the finite element solution, we obtain optimal error estimates of the fully discrete Crank- Nicolson/Newton scheme of nonlinear parabolic problem. Finally, numerical experiments are presented to show the efficient performance of the proposed scheme.

Pullback attractor of 2D nonautonomous g-Navier-Stokes equations with linear dampness

The pullback attractors for the 2D nonautonomous g-Navier-Stokes equations with linear dampness axe investigated on some unbounded domains. The existence of the pullback attractors is proved by verifying the existence of pullback D-absorbing sets with cocycle and obtaining the pullback ：D-asymptotic compactness. Furthermore, the estimation of the fractal dimensions for the 2D g-Navier-Stokes equations is given.

A New Closure to Slice Model for Slope Stability Analysis

This paper presents a new closure to slice models for evaluating slopes. The discussion is based on the minimal inter-slice action (MIA) hypothesis, which results in a new slice model without including artificially adjustable parameters. It has been realized that the new slice model predicts the minimum value of the safety factor, while all other slice models available always overestimate the value of the safety factor. Moreover, the gravity moment of each slice is found to be opposite to the overturning moment, which is different from the existing knowledge. In particular, the new slice model overcomes the situation where different assumptions of the inter-slice force function will give different safety factors to the same slope. The related numerical examples indicate that the new slice model can serve as a reliable tool for investigating geotechnical slope stability.

Velocity Projection with Upwind Scheme Based on the Discontinuous Galerkin Methods for the Two Phase Flow Problem

The upwind scheme is very important in the numerical approximation of some problems such as the convection dominated problem, the two-phase flow problem, and so on. For the fractional flow formulation of the two-phase flow problem, the Penalty Discontinuous Galerkin (PDG) methods combined with the upwind scheme are usually used to solve the phase pressure equation. In this case, unless the upwind scheme is taken into consideration in the velocity reconstruction, the local mass balance cannot hold exactly. In this paper, we present a scheme of velocity reconstruction in some H(div) spaces with considering the upwind scheme totally. Furthermore, the different ways to calculate the nonlinear coefficients may have distinct and significant effects, which have been investigated by some authors. We propose a new algorithm to obtain a more effective and stable approximation of the coefficients under the consideration of the upwind scheme.

Computational simulation of wave propagation problems in infinite domains

This paper deals with the computational simulation of both scalar wave and vector wave propagation problems in infinite domains. Due to its advantages in simulating complicated geometry and complex material properties, the finite element method is used to simulate the near field of a wave propagation problem involving an infinite domain. To avoid wave reflection and refraction at the common boundary between the near field and the far field of an infinite domain, we have to use some special treatments to this boundary. For a wave radiation problem, a wave absorbing boundary can be applied to the common boundary between the near field and the far field of an infinite domain, while for a wave scattering problem, the dynamic infinite element can be used to propagate the incident wave from the near field to the far field of the infinite domain. For the sake of illustrating how these two different approaches are used to simulate the effect of the far field, a mathematical expression for a wave absorbing boundary of high-order accuracy is derived from a two-dimensional scalar wave radiation problem in an infinite domain, while the detailed mathematical formulation of the dynamic infinite element is derived from a two-dimensional vector wave scattering problem in an infinite domain. Finally, the coupled method of finite elements and dynamic infinite elements is used to investigate the effects of topographical conditions on the free field motion along the surface of a canyon.

Computational Aspects of Multiscale Simulationwith the Lumped Particle Framework

First introduced in[2],the lumped particle framework is a flexible and numerically efficient framework for the modelling of particle transport in fluid flow.In this paper,the framework is expanded to simulate multicomponent particle-laden fluid flow.This is accomplished by introducing simulation protocols tomodel particles over a wide range of length and time scales.Consequently,we present a time ordering scheme and an approximate approach for accelerating the computation of evolution of different particle constituents with large differences in physical scales.We apply the extended framework on the temporal evolution of three particle constituents in sandladen flow,and horizontal release of spherical particles.Furthermore,we evaluate the numerical error of the lumped particle model.In this context,we discuss the Velocity-Verlet numerical scheme,and show how to apply this to solving Newton’s equations within the framework.We show that the increased accuracy of the Velocity-Verlet scheme is not lost when applied to the lumped particle framework.

Near-field performance of chemical dissolution-front instability around a circular acid-injection-well in fluid-saturated porous media

This paper investigates the near-field performance of chemical dissolution-front instability(CDFI) around a circular acidinjection-well in fluid-saturated porous media(FSPM) through using purely mathematical deductions.After the mathematical governing equations of the CDFI problem involving radially divergent flow are briefly described,both analytical base solutions and perturbation solutions for the considered problem are mathematically derived.These analytical solutions lead to the theoretical expression of the perturbation induced dimensionless growth-rate and the following two new findings.The first new finding is that the critical Peclet number of a chemical dissolution system(CDS) associated with radially divergent flow in FSPM is not only a function of the permeability ratio between the undissolved and dissolved regions as well as the dimensionless wavenumber,but also a function of the circular chemical dissolution-front location relative to the circular acid-injection-well in FSPM.The second new finding is that as the direct result of considering a nonzero radius of the circular acid-injection-well,there exits a critical closeness number,which may be used to assess where the circular chemical dissolution-front starts becoming unstable in the CDS associated with radially divergent flow.Based on these two new findings,a theoretical criterion of two parts has been established.The first part of the established theoretical criterion answers the scientific question when a circular chemical dissolution-front can become unstable,while the second part of the established theoretical criterion answers the scientific question where a circular chemical dissolution-front can become unstable.Through applying the established theoretical criterion,a long-term existing mystery why the wormhole pattern of fractal nature and the compact pattern of fingering nature are formed at different locations away from the circular acid-injection-well circumference in fluid-saturated carbonate rocks has been successfully revealed.