S. L. Butler

Forward modeling of applied geophysics methods using Comsol and comparison with analytical and laboratory analog models

Forward modeling is useful in geophysics both as a tool to interpret data in a research setting and as a tool to develop physical understanding in an educational setting. Gravity, magnetics, resistivity, and induced polarization are methods used in applied geophysics to probe Earths subsurface. In this contribution, we present forward models of these geophysical techniques using the finite-element modeling package Comsol. This package allows relatively easy implementation of these models and, as part of the AC/DC module, allows exterior boundaries to be placed at infinity, a boundary condition that is frequently encountered in geophysics. We compare the output of the finite-element calculations with analytical solutions and, for the resistivity method, with laboratory-scale analog experiments and demonstrate that these are in excellent agreement.

On scaling relations in time-dependent mantle convection and the heat transfer constraint on layering

We present a series of simulations of the mantle convection process based upon an axisymmetric numerical model and highlight a wide range of results in which scaling emerges. For the more challenging simulations it was found necessary to employ a finite difference mesh with uneven grid spacing in the radial coordinate, and we present the appropriate transformed field equations that are required to implement a model of this kind. The statistics of mass flux events transiting the 660-km phase transition are calculated for a large number of high-resolution calculations, and some of these are shown to display scale invariance properties in the high Rayleigh number regime. We also present a new parameterized model of convection and demonstrate its success in predicting the manner in which many of the bulk properties of the convection process scale with convection control parameters. Results are also presented which demonstrate that quantities such as heat flow, characteristic velocity, and thermal boundary layer thickness scale with the mean viscosity even in time-dependent simulations in which the effects of phase transitions, depth-varying viscosity, and internal heating are active. The heat flow scaling exponent is seen to decrease in magnitude with increasing internal heating rate and Clapeyron slope of the 660-km phase transition, but it is shown to be insensitive to depth variation of viscosity. Heat flow is seen to be reduced only modestly as the degree of layering increases unless layering is extreme. These calculations clearly demonstrate that in order for the surface heat flow predicted by the model to equal that characteristic of Earth, the mean viscosity of the mantle that controls the convection process must be considerably higher than the viscosity inferred on the basis of postglacial rebound and/or the flow must be significantly layered by the endothermic phase transition at 660 km depth. If mantle viscosity may be assumed to be Newtonian, in which case the creep resistance that controls rebound and convection must be the same, this constitutes a strong argument for the importance of layering. The force of this argument depends upon the existence of an accurate estimate of the temperature at the core-mantle boundary which has only recently become available.

Numerical modeling of fluid and electrical currents through geometries based on synchrotron X-ray tomographic images of reservoir rocks using Avizo and COMSOL

The use of numerical simulations to model physical processes occurring within subvolumes of rock samples that have been characterized using advanced 3D imaging techniques is becoming increasingly common. Not only do these simulations allow for the determination of macroscopic properties like hydraulic permeability and electrical formation factor, but they also allow the user to visualize processes taking place at the pore scale and they allow for multiple different processes to be simulated on the same geometry. Most efforts to date have used specialized research software for the purpose of simulations. In this contribution, we outline the steps taken to use commercial software Avizo to transform a 3D synchrotron X-ray-derived tomographic image of a rock core sample to an STL (STereoLithography) file which can be imported into the commercial multiphysics modeling package COMSOL. We demonstrate that the use of COMSOL to perform fluid and electrical current flow simulations through the pore spaces. The permeability and electrical formation factor of the sample are calculated and compared with laboratory-derived values and benchmark calculations. Although the simulation domains that we were able to model on a desk top computer were significantly smaller than representative elementary volumes, and we were able to establish Kozeny-Carman and Archies Law trends on which laboratory measurements and previous benchmark solutions fall. The rock core samples include a Fountainebleau sandstone used for benchmarking and a marly dolostone sampled from a well in the Weyburn oil field of southeastern Saskatchewan, Canada. Such carbonates are known to have complicated pore structures compared with sandstones, yet we are able to calculate reasonable macroscopic properties. We discuss the computing resources required. Graphical abstractStreamlines of simulated fluid flow through pore spaces of a carbonate rock imaged with synchrotron X-ray tomography. Colors indicate the magnitude of the flow velocity.Display Omitted HighlightsAvizo and Comsol are used to segment pore space and model fluid and electrical flows in rock pores.Pores are extracted from a synchrotron X-ray tomographic image of a carbonate rock core.Permeabilities and electrical formation factors agree with lab measurements and benchmarks.We describe the processing steps taken and the required computer resources.

Internal thermal boundary layer stability in phase transition modulated convection

The stability of a horizontal thermal boundary layer embedded within a very viscous fluid is investigated using the formalism of linear stability analysis. Thin thermal boundary layers in deep fluid regions and in the absence of phase transition and dynamical effects are thereby shown to be unstable at extremely long wavelengths. The stability of the internal thermal boundary layer which may exist at 660 km depth in the Earths mantle as a consequence of the dynamical influence of the endothermic phase transition from γ spinel to a mixture of perovskite and magnesiowustite, recently discussed in some detail by Solheim and Peltier [1994a], is investigated in order to better understand the “avalanche effect” observed in this and similar nonlinear, time dependent simulations of the mantle convection process. It is demonstrated that if the stability problem is treated as purely thermal, then the boundary layer is predicted to be extremely unstable and the presence of the 660-km endothermic phase transition at middepth within the boundary layer is further destabilizing. When the kinematic effect of flow convergence onto the boundary layer and phase transition region is active, however, it is shown that the layer may be strongly stabilized. In the regime of physically realistic velocity convergence, the critical Rayleigh number is predicted to lie in the range suggested by the numerical simulations of Solheim and Peltier [1994a]. A threshold value of the magnitude of the Clapeyron slope of the endothermic phase transition for a given velocity convergence is also shown to exist, beyond which the fastest-growing mode of instability changes from avalanche type to layered type.

On the origin and significance of subadiabatic temperature gradients in the mantle

[1] It is well established that the temperature gradients in the interiors of internally heated mantle convection models are subadiabatic. The subadiabatic gradients have been explained as arising because of a balance between vertical advection and internal heating; however, a detailed analysis of the energy balance in the subadiabatic regions has not been undertaken. In this paper, we examine in detail the energy balance in a suite of simple, two-dimensional convection calculations with mixed internal and basal heating, depthdependent viscosity, and continents. We find that there are three causes of subadiabatic gradients. One is the above mentioned balance, which becomes significant when the ratio of internal heating to total surface heat flow is large. The second mechanism involves the growth of the ‘‘overshoot’’ of the geotherm near the lower boundary where the dominant balance is between vertical and horizontal advection. The latter mechanism is significant even in relatively weakly internally heated calculations. For time-dependent calculations, we find that local secular cooling can be a dominant term in the energy equation and can lead to subadiabaticity. However, it does not show its signature on the shape of the time-averaged geotherm. We also compare the basal heat flow with parameterized calculations based on the temperature drop at the core-mantle boundary, calculated both with and without taking the subadiabatic gradient into account, and we find a significantly improved fit with its inclusion.

The shape distribution of splash-form tektites predicted by numerical simulations of rotating fluid drops

Splash-form tektites are glassy rocks ranging in size from roughly 1 to 100 mm that are believed to have formed from the splash of silicate liquid after a large terrestrial impact from which they are strewn over thousands of kilometres. They are found in an array of shapes including spheres, oblate ellipsoids, dumbbells, rods and possibly fragments of tori. It has recently become appreciated that surface tension and centrifugal forces associated with the rotation of fluid droplets are the main factors determining the shapes of these tektites. In this contribution, we compare the shape distribution of 1163 measured splash-form tektites with the results of the time evolution of a 3D numerical model of a rotating fluid drop with surface tension. We demonstrate that many aspects of the measured shape distribution can be explained by the results of the dynamical model.

Effective transport rates and transport-induced melting and solidification in mushy layers

Mushy layers are known to occur in magma chambers, sea-ice, and metal castings. They are often modeled as a porous layer in which a fluid and solid matrix exist in thermal and compositional equilibria. In nonreactive porous media, both advective and diffusive transport rates for heat and solute differ. In mushy layers, however, the temperature and composition of the fluid phase are constrained by the liquidus relationship giving rise to effective transport rates that are intermediate between those for heat and solute in passive porous layers. The transport of heat and solute, even if the invading fluid is itself in equilibrium, is also accompanied by a degree of solidification or melting due to the difference in the transport rates for these two quantities. In this paper, analytical expressions for the effective velocity and diffusivity in a mushy layer and for the degree of melting or solidification accompanying the passage of a front with a different temperature and composition are derived and compared ...

Linear analysis of melt band formation in a mid‐ocean ridge corner flow

Imposing external shear on systems of partial melt can cause compaction of the solid matrix and concentration of the interstitial liquid melt thereby generating bands of contrasting high and low porosity. These shear-induced porosity bands have been proposed to channel melt beneath a mid-ocean ridge (MOR). In this contribution, we impose a shear flow that evolves as the bands are transported along streamlines of MOR corner flow. We evaluate the suitability of porosity band formation as a mechanism for melt channeling beneath a MOR using a linear instability analysis with three different matrix viscosity conditions: isotropic strain rate independent, isotropic strain rate dependent, and anisotropic. Our analysis shows that the largest amplitude bands channel melt away from the ridge axis toward the base of the plate at the lithosphere-asthenosphere boundary, while the effectiveness of channeling through the bands is highly dependent on the mantle bulk viscosity.