Stefan M. Schmalholz

Buckling versus folding: Importance of viscoelasticity

We present a dominant wavelength solution for a viscoelastic layer embedded in a low viscosity matrix under layer-parallel compression based on the thin-plate approximation. We show that the deformation mode approximates the elastic or viscous limits depending on a parameter, R, which is the ratio of dominant wavelength predicted by pure viscous theory to the one predicted by pure elastic theory. In contrast, conventional analyses based on the Deborah number incorrectly predict the deformation mode. The dominant viscoelastic wavelength closely follows the minimum out of viscous and elastic dominant wavelengths. The viscoelastic thin-plate theory is verified by two-dimensional modeling of large strain viscoelastic folding, for which we develop a new numerical algorithm based on a combined spectral/finite-difference method. The robustness of the numerical code is demonstrated by calculation, for the first time, of the pressure field evolution during folding of a viscoelastic layer with up to 100% strain.

A passive seismic survey over a gas field: Analysis of low-frequency anomalies

Passive seismic low-frequency from approximately 1–6 Hz data have been acquired at several locations around the world. Spectra calculated from these data, acquired over fields with known hydrocarbon accumulations, show common spectral anomalies. Verification of whether these anomalies are common to only a few, many, or all hydrocarbon reservoirs can be provided only if more and detailed results are reported. An extensive survey was carried out above a tight gas reservoir and an adjacent exploration area in Mexico. Data from several hundred stations with three-component broadband seismometers distributed over approximately 200 km 2 were used for the analysis. Several hydrocarbon reservoir-related microtremor attributes were calculated, and mapped attributes were compared with known gas intervals, with good agreement. Wells drilled after the survey confirm a predicted high hydrocarbon potential in the exploration area. A preliminary model was developed to explain the source mechanism of those microtremors. Poroelastic effects caused by wave-induced fluid flow and oscillations of different fluid phases are significant processes in the low-frequency range that can modify the omnipresent seismic background spectrum. These processes only occur in partially saturated rocks. We assume that hydrocarbon reservoirs are partially saturated, whereas the surrounding rocks are fully saturated. Our real data observations are consistent with this conceptual model.

Low-frequency reflections from a thin layer with high attenuation caused by interlayer flow

The 1D interlayer-flow (or White’s) model is based on Biot’s theory of poroelasticity and explains low-frequency seismic wave attenuation in partially saturated rocks by wave-induced fluid flow between two alternating poroelastic layers, each saturated with a different fluid. We have developed approximate equations for both the minimum possible value of the quality factor, Q , and the corresponding fluid saturation for which Q is minimal. The simple approximate equations provide a better insight into the dependence of Q on basic petrophysical parameters and allow for a fast assessment of the minimal value of Q . The approximation is valid for a wide range of realistic petrophysical parameter values for sandstones partially saturated with gas and water, and shows that values of Q can be as small as two. We ap-plied the interlayer-flow model to study the reflection coefficient of a thin (i.e., between 6 and 11 times smaller than the incident wavelength) layer that is partially saturated with gas and water. ...

Finite amplitude folding: transition from exponential to layer length controlled growth

Abstract A new finite amplitude theory of folding has been developed by the combined application of analytical, asymptotic and numerical methods. The existing linear folding theory has been improved by considering nonlinear weakening of membrane stresses, which is caused by the stretching of the competent layer during folding. The resulting theory is simple and accurate for finite amplitude folding and is not restricted to infinitesimal amplitudes, as is the classical linear theory of folding. Two folding modes relevant to most natural settings were considered: (i) both membrane and fiber stresses are viscous during folding (the ‘viscous’ mode); (ii) membrane stresses are viscous whereas fiber stresses are elastic (the ‘viscoelastic’ mode). For these two modes, the new theory provided a nonlinear, ordinary differential equation for fold amplification during shortening and an estimate for crossover amplitude and strain where the linear theory breaks down. A new analytical relationship for amplitude versus strain was derived for strains much larger than the crossover strain. The new relationship agrees well with complete 2D numerical solutions for up to threefold shortening, whereas the exponential solution predicted by the linear theory is inaccurate by orders of magnitude for strains larger than the crossover value. Analysis of the crossover strain and amplitude as a function of the controlling parameters demonstrates that the linear theory is only applicable for a small range of amplitudes and strains. This renders unreliable the large strain prediction of wavelength selection based on the linear theory, especially for folding at high competence contrasts. To resolve this problem, the new finite amplitude theory is used to calculate the evolution of the growth rate spectra during progressive folding. The growth rate spectra exhibited splitting of a single maximum (predicted by the linear theory) into two maxima at large strains. This bifurcation occurred for both deformation modes. In contrast, the spectra of the cumulative amplification ratio (current over initial amplitude) maintained a single maximum value throughout. The wavelength selectivity is found to decrease at large strains, which helps explain the aperiodic forms of folds commonly observed in nature and the absence of long dominant wavelengths for high competence contrast folding. Calculation of the cumulative amplification spectra for different initial amplitude distributions, ranging from white to red noise, showed that the initial noise has a strong influence on the amplitude spectra for small strains. For larger strains, however, the cumulative amplification spectra were similar despite the strong difference in the initial noise.

Stress-strength relationship in the lithosphere during continental collision

Lithospheric strength profiles generated for a shortening continental lithosphere generally predict excessively high differential stresses in the sub-Moho continental mantle; this seems inconsistent with the relative scarcity of earthquakes at this depth. This inconsistency was put forward as evidence for weak mantle rheology. However, this argument implicitly assumes that strength envelopes are valid in actively deforming regions. We test this assumption on two end-member model lithospheres having identical upper crustal rheologies, but with (1) a weak lower crust and strong mantle, and (2) a strong lower crust and weak mantle. For this purpose, we compare one-dimensional (1-D) with 2-D visco-elastoplastic numerical models of continental shortening. The 2-D models show that strongly heterogeneous deformation typically follows initially homogeneous deformation. Lithospheric-scale buckle folds and shear zones result in strain rate variations of as much as three orders of magnitude. Differential stresses in the upper crust are close to yield, as predicted by 1-D models. Stresses in deeper lithospheric regions, however, are significantly smaller than in 1-D models, especially in actively deforming regions. Systematic numerical simulations as a function of temperature and deformation rate reveal that 1-D models are reliable in hot and/or slowly deforming lithospheres only. The relative scarcity of earthquakes at mantle levels should thus be interpreted as an intrinsic consequence of strong lithospheric deformation rather than as evidence for a weak upper mantle.

Finite-difference modeling of wave propagation on microscale: A snapshot of the work in progress

Digital rock methodology combines modern microscopic imagingwithadvancednumericalsimulationsofthephysicalproperties of rocks. Modeling of elastic-wave propagation directly from rock microstructure is integral to this technology. We survey recent development of the rotated staggered grid RSG finite-difference FD method for pore-scale simulation of elastic wavepropagationindigitalrocksamples,includingthedynamic elastic properties of rocks saturated with a viscous fluid. Examination of the accuracy of this algorithm on models with known analytical solutions provide an additional accuracy condition for numerical modeling on the microscale. We use both the elastic and viscoelastic versions of the RSG algorithm to study gas hydratesandtosimulatepropagationofBiot’sslowwave.Weapply RSG method ology to examine the effect of gas hydrate distributions in the pore space of a rock. We compare resulting P-wave velocities with experimentally measured data, as a basis for buildinganeffective-mediummodelforrockscontaininggashydrates. We then perform numerical simulations of Biot’s slow wave in a realistic 3D digital rock model, fully saturated with a nonviscous fluid corresponding to the high-frequency limit of poroelasticity, and placed inside a bulk fluid.The model clearly demonstrates Biot’s slow curve when the interface is open between the slab and bulk fluid.We demonstrate slow wave propagation in a porous medium saturated with a viscous fluid by analyzing an idealized 2D porous medium represented alternating solid and viscous fluid layers. Comparison of simulation results withtheexactsolutionforthislayeredsystemshowsgoodagreementoverabroadfrequencyrange.

Strain and competence contrast estimation from fold shape

A new method to estimate strain and competence contrast from natural fold shapes is developed and verified by analogue and numerical experiments. Strain is estimated relative to the nucleation amplitude, AN, which is the fold amplitude when the amplification velocities caused by kinematic layer thickening and dynamic folding are identical. AN is defined as the initial amplitude corresponding to zero strain because folding at amplitudes smaller than AN is dominantly by kinematic layer thickening. For amplitudes larger than AN, estimates of strain and competence contrast are contoured in thickness-to-wavelength (H/l) and amplitude-to-wavelength (A/l) space. These quantities can be measured for any observed fold shape. Contour maps are constructed using existing linear theories of folding, a new nonlinear theory of folding and numerical simulations, all for single-layer folding. The method represents a significant improvement to the arc length method. The strain estimation method is applied to folds in viscous (Newtonian), power-law (non-Newtonian) and viscoelastic layers. Also, strain partitioning in fold trains is investigated. Strain partitioning refers to the difference in strain accommodated by individual folds in the fold train and by the whole fold train. Fold trains within layers exhibiting viscous and viscoelastic rheology show different characteristic strain partitioning patterns. Strain partitioning patterns of natural fold trains can be used to assess the rheological behaviour during fold initiation. D 2001 Elsevier Science B.V. All rights reserved.

3D finite amplitude folding: Implications for stress evolution during crustal and lithospheric deformation

viscosity contrast. The results show that particular combinations of horizontal wavelengths (particular modes) grow faster than others; these ‘‘dominant’’ modes are the ones expected to develop in nature. Different than in the 2D case, wheretypicallyonlyonedominantmodeexists,in3Darange of modes amplify at nearly equal rate. The linear stability analysis is strictly valid only for infinitesimal amplitudes and assumes that the fold amplitude grows exponentially with horizontal strain (or deformation time). We perform numericalsimulations(finiteelementmethod)of3Dviscoussinglelayer folding to calculate the evolution of the fold amplitude with progressive shortening. The 2D analytical solution describing finite amplitude folding [Schmalholz and Podladchikov, 2000] is modified to accurately describe 3D finite amplitude folding. The solution is extended by an additional parameter, which is the ratio of the shortening strain rates in both horizontal directions. The numerical simulations are further applied to calculate the evolution of the averaged differential stress during 3D folding.

Automated thermotectonostratigraphic basin reconstruction: Viking Graben case study

We present a generic algorithm for automating sedimentary basin reconstruction. Automation is achieved through the coupling of a two-dimensional thermotectonostratigraphic forward model to an inverse scheme that updates the model parameters until the input stratigraphy is fitted to a desired accuracy. The forward model solves for lithospheric thinning, flexural isostasy, sediment deposition, and transient heat flow. The inverse model updates the crustal- and mantle-thinning factors and paleowater depth. Both models combined allow for automated forward modeling of the structural and thermal evolution of extensional sedimentary basins. The potential and robustness of this method is demonstrated through a reconstruction case study of the northern Viking Graben in the North Sea. This reconstruction fits present stratigraphy, borehole temperatures, vitrinite reflectance data, and paleowater depth. The predictive power of the model is illustrated through the successful identification of possible targets along the transect, where the principal source rocks are in the oil and gas windows. These locations coincide well with known oil and gas occurrences. The key benefits of the presented algorithm are as follows: (1) only standard input data are required, (2) crustal- and mantle-thinning factors and paleowater depth are automatically computed, and (3) sedimentary basin reconstruction is greatly facilitated and can thus be more easily integrated into basin analysis and exploration risk assessment.

Brittle fracture during folding of rocks: A finite element study

The goal of the present work is the development of a novel computational analysis tool to elaborate folding-induced fracture of geological structures. Discrete failure of brittle rocks is characterised by three sets of governing equations: the bulk problem, the interface problem and the crack problem. The former two sets which define the deformation field are highly nonlinear and strongly coupled. They are solved iteratively within a Hansbo-type finite element setting. The latter set defines the crack kinematics. It is linear and solved in a single post-processing step. To elaborate the features of the computational algorithm, we define a unique benchmark problem of a single, geometrically nonlinear plate, which is subjected to layer-parallel in-plane compression combined with different levels of superposed in-plane shear. The resulting folding, or buckling, induces brittle failure in the tensile regime. By systematically increasing the shear strain at constant compression, we develop crack deviation angle versus shear-to-compression ratio tables. We determine the corresponding damage zones, analyse the folding modes and elaborate the force versus amplification diagrams. The proposed two-field folding-induced fracture algorithm can ultimately be applied to interpret natural folded rocks and understand their evolution, structural development and histology.