Numerical Simulations of Near-Surface Stress Perturbation Caused by Erosion-Induced Isostatic Rebound of Viscoelastic Lithosphere
In plate tectonics, the relative motion between plates at a plate boundary generates tectonic stress, and most earthquakes occur at these plate boundaries. However, only a small proportion of earthquakes occur in plate interiors, such as the 1811–12 New Madrid Earthquake (8.1 Mw) in the U.S.A., and the 2008 Sichuan Earthquake (7.9 Mw) in China. These earthquakes, which occur within a plate, are called intraplate earthquakes. Numerous hypotheses have been suggested as the mechanisms that explain intraplate earthquakes, such as local stress perturbation due to short wavelength factors (e.g., gravitational body forces, isostatic rebound) and stress interactions between faults. Among these factors, isostatic rebound caused by the erosion of sedimentary layers induces stress perturbation at the surface. The magnitude of the stress variation ranges from tens to hundreds of Pa/yr, which can yield displacement of several meters for tens of thousands of years. In this study, numerical simulations of this stress perturbation using a commercial finite element package, known as COMSOL Multiphysics®, confirmed that this model can accurately simulate the values of stress variation and displacement in a natural environment. In our two-dimensional (plane strain) viscoelastic model, we divided the spatial domain into the upper crust, lower crust, and mantle, and assigned mechanical properties to these domains, including the Young’s modulus and viscosity. While the upper crust and lower crust is a purely elastic material, the upper mantle is characterized as a Maxwell viscoelastic material with temperature dependent power-law viscosity. We used the Heat Transfer Module in COMSOL Multiphysics® to set the initial temperature structure as a general geotherm inferred from conventional surface heat flux over the continental crust. Using the Interpolation Function provided by COMSOL Multiphysics®, we generated a temperature-dependent power-law function to calculate the viscosity. Since viscosity depends on temperature, viscosity decreases as depth increases from the surface to upper mantle by a factor of 10,000. To simulate the sedimentary erosion process, the stress that corresponds to the weight of a 12 m high sediment layer was gradually removed during 16,000 yr from the top of the model using the 'step function'. To quantify how the rheological structure of the mantle contributes to the deformation of the elastic upper crust, we multiplied the mantle viscosity in the reference model by 0.001–1,000 using a Parametric sweep.
Our results show that the displacement and stress variation rates at the surface are consistent with previous studies, which indicates that stress variations with a small magnitude can be accurately simulated using COMSOL Multiphysics®. Our model showed that the different values of mantle viscosity affect the magnitude of the displacement and stress rate at the surface. For instance, models at 0.001 and 1,000 times the reference case exhibit 50 m of displacement at 0.193 kPa/yr and 31 m of displacement at 0.157 kPa/yr, respectively. These results indicate the potential that our model has to predict mantle viscosity with geophysical and geodetic observations of surficial displacement.
