This study investigates the impact of different slip-dependent fault failure parametrizations on the pre-seismic and coseismic phases of the earthquake cycle using 2-D finite element simulations. Various failure laws are considered: linear slip-weakening, double slip-weakening, parabolic cohesive zone and exponential cohesive zone laws. The pre-seismic phase is modelled using a quasi-static approach, while the coseismic phase is modelled dynamically. Results demonstrate the importance of accounting for the different shapes of failure laws during the pre-seismic phase, as they lead to variations in pre-seismic phase duration, associated pre-seismic slip and lateral extent of slip front. Failure laws with gentler initial slopes require more time to initiate dynamic rupture compared to those with steeper initial weakening slopes. The presence of an initial strengthening segment in the failure law affects the amount and lateral extent of pre-seismic slip. During dynamic rupture propagation, the specific details of the failure law become less significant due to the Lorentz contraction of the process zone. However, variations in the duration of the pre-seismic phases result in different dynamic stress drops and, consequently, variations in rupture acceleration times and earthquake magnitudes for different failure laws.

The failure law prescribed along the fault surface and the elastic stiffness of the surround-ing medium play important roles in determining the characteristics of earthquakes. Here we use a 1D spring-slider model that includes inertia, along with a simple poly-linear fail-ure law composed of multiple linear segments to provide insight into earthquake initiation and growth. The poly-linear failure law, which parameterizes shear resistance as a function of slip, allows analytical solutions describing the system for each failure law segment. Analytical solutions facilitate investigation of the effects of the slopes of the different fail-ure law segments in relation to the slope of the elastic loading curve determined by the spring stiffness. Depending on the relation between the slope of the failure law segment and the elastic loading slope, there are three stability regimes in the system: harmonic oscillations, exponential growth, and cubic growth. By combining the different solution regimes within one earthquake cycle, we observe a wide range of behaviors of this simple system: interseismic oscillatory creep, precursory signals before the main event, a shorter or a much longer acceleration phase before the onset of instability, and varying durations of the preseismic and coseismic phases. These results provide a potential explanation for some seismic observations, including increased levels of "seismic noise" prior to an earth-quake, precursory events, tremor and low-frequency earthquakes.

SUMMARY By anatomizing the classic Liouville equation (LE), an alternative theoretical framework is established for the Earth polar motion where the turbulent time derivative of the fluid forcing, Lforce, is eliminated. The observed polar motion is found a lumped signal of two physical variables, one of which is the forced polar motion, mforce, that balances out the fluid forcing through a simple algebraic identity. The second component is the inertial polar motion, minert, which conserves the total angular momentum by the restoring power of the equatorial bulge. Aside from its numerical difficulties, the time derivative dLforce/dt has proved to be a redundant artefact that mixes the physical signals in the solution of the LE. By an analytical appraisal, we find that for non-Chandler forced polar motions, the dLforce/dt term in the excitation function of the classic LE is 100 per cent responsible for the forced term, mforce, and contributes 0.3 per cent of the inertial term, minert. The Chandler signal originates exclusively from the inertial polar motion minert. The artefact dLforce/dt perturbs the true Chandler signal from the classic LE by 0.3 per cent. Anatomy of the LE allows us to show that a linearized mechanism for the Chandler Wobble excitation is valid, since the linearized governing equation for the inertial polar motion minert stands as part of the exact solution of the full nonlinear equation. Anatomized polar motion equations are also considered in the presence of ocean tides, under the lunisolar torque, and on an elastically deformable Earth. Inaccuracies and misinterpretations associated with the classic LE under those circumstances are clarified in formulating the anatomized polar motion equation group.

There is growing concern about seismicity triggered by human activities, whereby small increases in stress bring tectonically loaded faults to failure. Examples of such activities include mining, impoundment of water, stimulation of geothermal fields, extraction of hydrocarbons and water, and the injection of water, CO 2 and methane into subsurface reservoirs 1 . In the absence of sufficient information to understand and control the processes that trigger earthquakes, authorities have set up empirical regulatory monitoring-based frameworks with varying degrees of success 2 , 3 . Field experiments in the early 1970s at the Rangely, Colorado (USA) oil field 4 suggested that seismicity might be turned on or off by cycling subsurface fluid pressure above or below a threshold. Here we report the development, testing and implementation of a multidisciplinary methodology for managing triggered seismicity using comprehensive and detailed information about the subsurface to calibrate geomechanical and earthquake source physics models. We then validate these models by comparing their predictions to subsequent observations made after calibration. We use our approach in the Val d’Agri oil field in seismically active southern Italy, demonstrating the successful management of triggered seismicity using a process-based method applied to a producing hydrocarbon field. Applying our approach elsewhere could help to manage and mitigate triggered seismicity.

Acoustic emission (AE) is a widely used technology to study source mechanisms and material properties during high-pressure rock failure experiments. It is important to understand the physical quantities that acoustic emission sensors measure, as well as the response of these sensors as a function of frequency. This study calibrates the newly built AE system in the MIT Rock Physics Laboratory using a ball-bouncing system. Full waveforms of multibounce events due to ball drops are used to infer the transfer function of lead zirconate titanate (PZT) sensors in high pressure environments. Uncertainty in the sensor transfer functions is quantified using a waveform-based Bayesian approach. The quantification of in situ sensor transfer functions makes it possible to apply full waveform analysis for acoustic emissions at high pressures.

We present a crosslink constraint method for numerically modeling dynamic slip on intersecting faults, without prescribing slip (dis-)continuation directions. The fault intersections are constrained by crosslinked split nodes, such that the slip can only be continuous on one of the two intersecting faults at a time and location. The method resolves the episodic intersection offset by examining the dynamic fault traction resulting from two sets of constraint equations, one for each slip direction. To verify this method, we modify two benchmark problems, hosted at Southern California Earthquake Center (SCEC), by allowing a branching fault to step across a main fault. The modified SCEC problem results agree with our expectations that the intersection offset scenarios are dictated by the nucleation patch location and initial fault traction. This new method comes with an open-source finite-element code Defmod.

We present a fundamental solution‐based finite‐element (FE) method to homogenize heterogeneous elastic medium, that is, fault zone, under static, and dynamic loading. This method incorporates Eshelby’s strain perturbation into FE weak forms. The resulting numerical model implicitly considers the existence of inhomogeneity bodies within each element, without introducing additional degrees of freedom. The new method is implemented within an open‐source FE package that is applicable to alternating seismic and aseismic cycles. To demonstrate this method, we modify a dynamic fault‐slip problem, hosted at Southern California Earthquake Center (SCEC), by introducing a fault zone that contains different microstructures than the host matrix. The preliminary results suggest that the fault‐zone microstructure orientation has effects on fault slip, seismic arrivals and waveform frequency contents.

We extend the Eshelby’s (equivalent) inclusion method to consider ellipsoidal volumes occupied by fluids of limited compressibility. The new method decomposes the problem into two parts – a stress free transformation problem, and a volume compatibility problem. Implemented in an open source computer code, Esh3D, the new method allows one to model interacting and fluid–solid mixed inclusions in whole and truncated spaces. We verify the new method against finite element (FE) method, showing it to be accurate and efficient.

Microseismic monitoring is generally the most reliable method for estimating stimulated fractured volume. Receivers used in microseismic monitoring measure only seismic events. That limitation explains why only a small portion of the energy budget during hydraulic fracturing can be estimated by information obtained from microseismic monitoring. We performed a series of numerical experiments to investigate the effects of rock mechanical properties and fracture friction characteristics on seismic efficiency and rupture velocity. We conducted numerical experiments using acoustic emission for saw-cut samples under triaxial loads and applied slip-weakening constitutive modeling for natural fractures to study how the Young's modulus and slip-weakening distance affect seismic efficiency and rupture velocity. Perhaps surprisingly, our results show that rocks with higher values of the Young's modulus have lower seismic efficiency generated from sliding on pre-existing natural fractures, while lower rigidity leads to higher seismic efficiency. These results do not contradict general beliefs about the effect of rigidity on fracability. More rigid rocks are more favorable for hydraulic fracturing and generate larger fracture networks; however, compared with less rigid rocks, fewer events would be detected seismically. The results also give insight into how to connect geomechanical numerical modeling of hydraulic fractures in naturally fractured reservoirs with microseismic data from the field and actual subsurface-generated fractured networks.

This paper presents a mortar-based finite element formulation for modeling the dynamics of shear rupture on rough interfaces governed by slip-weakening and rate and state (RS) friction laws, focusing on the dynamics of earthquakes. The method utilizes the dual Lagrange multipliers and the primal–dual active set strategy concepts, together with a consistent discretization and linearization of the contact forces and constraints, and the friction laws to obtain a semi-smooth Newton method. The discretization of the RS friction law involves a procedure to condense out the state variables, thus eliminating the addition of another set of unknowns into the system. Several numerical examples of shear rupture on frictional rough interfaces demonstrate the efficiency of the method and examine the effects of the different time discretization schemes on the convergence, energy conservation, and the time evolution of shear traction and slip rate.

We study the response to slow tectonic loading of rough faults governed by velocity weakening rate and state friction, using a 2‐D plane strain model. Our numerical approach accounts for all stages in the seismic cycle, and in each simulation we model a sequence of two earthquakes or more. We focus on the global behavior of the faults and find that as the roughness amplitude, b r , increases and the minimum wavelength of roughness decreases, there is a transition from seismic slip to aseismic slip, in which the load on the fault is released by more slip events but with lower slip rate, lower seismic moment per unit length, M 0,1 d , and lower average static stress drop on the fault, Δ τ t . Even larger decreases with roughness are observed when these source parameters are estimated only for the dynamic stage of the rupture. For b r ≤ 0.002, the source parameters M 0,1 d and Δ τ t decrease mutually and the relationship between Δ τ t and the average fault strain is similar to that of a smooth fault. For faults with larger values of b r that are completely ruptured during the slip events, the average fault strain generally decreases more rapidly with roughness than Δ τ t .

We study numerically the effects of fault roughness on the nucleation process during earthquake sequences. The faults are governed by a rate and state friction law. The roughness introduces local barriers that complicate the nucleation process and result in asymmetric expansion of the rupture, nonmonotonic increase in the slip rates on the fault, and the generation of multiple slip pulses. These complexities are reflected as irregular fluctuations in the moment rate. There is a large difference between first slip events in the sequences and later events. In the first events, for roughness amplitude br ≤ 0.002, there is a large increase in the nucleation length with increasing br. For larger values of br, slip is mostly aseismic. For the later events there is a trade‐off between the effects of the finite fault length and the fault roughness. For br ≤ 0.002, the finite length is a more dominant factor and the nucleation length barely changes with br. For larger values of br, the roughness plays a larger role and the nucleation length increases significantly with br. Using an energy balance approach, where the roughness is accounted for in the fault stiffness, we derive an approximate solution for the nucleation length on rough faults. The solution agrees well with the main trends observed in the simulations for the later events and provides an estimate of the frictional and roughness properties under which faults experience a transition between seismic and aseismic slip.

由流体的注入和采出所诱发的地震活动已成为一个围绕地下水资源和能源的科学讨论热点.本文展示了流动和地质力学耦合的模拟技术在2012年5月意大利北部卡沃内油田附近发生的破坏性地震(M W 6.0和M W 5.8)序列的震后分析中的应用.该序列引发了一个问题:这些地震是否可能由石油和天然气生产活动诱发.本文的分析强有力地表明,卡沃内油田流体的采出和注入的联合效应并不是观测到的地震活动的诱因.更普遍地,本文研究表明耦合流动和地质力学的计算模型可将地质、地震构造、测井、流体压力、流动速率和大地测量数据整合,并为评估和管理与诱发型地震活动相关的危害提供一种有希望的方法.

Seismicity induced by fluid injection and withdrawal has emerged as a central element of the scientific discussion around subsurface technologies that tap into water and energy resources. Here we present the application of coupled flow-geomechanics simulation technology to the post mortem analysis of a sequence of damaging earthquakes (M-w=6.0 and 5.8) in May 2012 near the Cavone oil field, in northern Italy. This sequence raised the question of whether these earthquakes might have been triggered by activities due to oil and gas production. Our analysis strongly suggests that the combined effects of fluid production and injection from the Cavone field were not a driver for the observed seismicity. More generally, our study illustrates that computational modeling of coupled flow and geomechanics permits the integration of geologic, seismotectonic, well log, fluid pressure and flow rate, and geodetic data and provides a promising approach for assessing and managing hazards associated with induced seismicity.

Ceres, the largest body in the asteroid belt (940 km diameter), provides a unique opportunity to study the interior structure of a volatile-rich dwarf planet. Variations in a planetary body's subsurface rheology and density affect the rate of topographic relaxation. Preferential attenuation of long wavelength topography (>= 150 km) on Ceres suggests that the viscosity of its crust decreases with increasing depth. We present finite element (FE) geodynamical simulations of Ceres to identify the internal structures and compositions that best reproduce its topography as observed by the NASA Dawn mission. We infer that Ceres has a mechanically strong crust with maximum effective viscosity similar to 10(25) Pas. Combined with density constraints, this rheology suggests a crustal composition of carbonates or phyllosilicates, water ice, and at least 30 volume percent (vol.%) low-density, high-strength phases most consistent with salt and/or clathrate hydrates. The inference of these crustal materials supports the past existence of a global ocean, consistent with the observed surface composition. Meanwhile, we infer that the uppermost >= 60 km of the silicate-rich mantle is mechanically weak with viscosity <10(21) Pas, suggesting the presence of liquid pore fluids in this region and a low temperature history that avoided igneous differentiation due to late accretion or efficient heat loss through hydrothermal processes. (C) 2017 Elsevier B.V. All rights reserved.

Faults are rough at all scales and can be described as self-affine fractals. This deviation from planarity results in geometric asperities and a locally heterogeneous stress field, which affect the nucleation and propagation of shear rupture. We study this effect numerically at the scale of small earthquakes, in which realistic geometry and friction law parameters can be incorporated in the model. We aim to understand the effect of roughness on faults with L ~ 10 – 200 m. At this scale we can choose the minimum roughness wavelength, λmin, to be the size of lab samples (5 – 10 cm) and thus use the observed lab-scale slip or rate based friction laws without modifying the constitutive parameters, assuming that the experimental data already include the effects of smaller wavelengths of roughness. Moreover, using a variable time step size, we gradually increase the remote stress and let the rupture nucleate spontaneously, rather than introducing artificial procedures to nucleate a seismic event. Numerically, maintaining λmin and consequently the smallest element size, Δx, fixed while increasing the fault length poses two challenges. First, keeping Δx fixed is computationally expensive. Second, the slip increases with increasing fault length and the assumption of small slip relative to the size of the elements is not valid. To overcome the first challenge, we refine the mesh near the fault, using hanging nodes. To overcome the second challenge, we use the Mortar finite element method [Bernardi et al., 1994; Wohlmuth, 2000], in which non-matching meshes are allowed across the fault and the contacts are continuously updated. We introduce slip weakening and rate and state friction laws into the method and study both the nucleation and propagation of shear rupture, using variable time steps with a quasistatic scheme for the inter-seismic stage and a dynamic implicit Newmark scheme for the co-seismic stage. For a static benchmark, we demonstrate that the method predicts accurately the stresses and displacements along a fault with a non-matching grid due to a uniform stress drop. We also design a benchmark problem to show that the method accurately models the behavior of the friction coefficient in response to a change in the slip rate on a fault governed by a rate and state friction law. Simulations of a 10 meter long horizontal fault with different amplitude roughness and a slip-weakening friction law show the significant effect of roughness on: (1) Slip on the fault and consequently the seismic moment; (2) Stress drop; (3) Rupture properties, such as rupture velocity, breakdown zone, and the observed relation between shear stress and slip; and (4) Different stages in the nucleation processes. For example, with the adopted spontaneous nucleation procedure, the experiment-based nucleation model of Ohnaka [2000] is observed also in the simulations (Fig. 1), and important quantities regarding the nucleation and propagation of the rupture can be measured.

SummaryCharacterization of reservoir properties like porosity and permeability in reservoir models typically relies on history matching of production data, well pressure data, and possibly other fluid‐dynamical data. Calibrated (history‐matched) reservoir models are then used for forecasting production and designing effective strategies for improved oil and gas recovery. Here, we perform assimilation of both flow and deformation data for joint inversion of reservoir properties. Given the coupled nature of subsurface flow and deformation processes, joint inversion requires efficient simulation tools of coupled reservoir flow and mechanical deformation. We apply our coupled simulation tool to a real underground gas storage field in Italy. We simulate the initial gas production period and several decades of seasonal natural gas storage and production. We perform a probabilistic estimation of rock properties by joint inversion of ground deformation data from geodetic measurements and fluid flow data from wells. Using an efficient implementation of the ensemble smoother as the estimator and our coupled multiphase flow and geomechanics simulator as the forward model, we show that incorporating deformation data leads to a significant reduction of uncertainty in the prior distributions of rock properties such as porosity, permeability, and pore compressibility. Copyright © 2015 John Wiley & Sons, Ltd.