Search

This is data and code for “Mesh size effect on finite source inversion with 3D finite element modeling” contents 1. Green's function matrix and synthetic observation data for checkerboard test from simplified model 2. Green's function matrix and synthetic observation data for checkerboard test from regional model 3. Green's function matrix and observed coseismic displacements at the GPS stations for 2011 Mw 9.0 Tohoku-Oki earthquake The details of each file are described in readme.txt

Rifting is a large-scale planetary evolution process that forms a new oceanic crust with mid-ocean ridges in an extensional environment. In this process, mantle convection occurred and material circulates, forming a volcano in the surrounding area. It is well known that mantle flow of rifting causes volcanism, but most of the volcanic processes are concentrated in the mid-ocean ridge and rift center axis. Recently, many of theory (lithosphere delamination, edge-driven convection, slab tearing, etc.) have been discussed to explain intra-plate volcanic mechanisms at non-plume and non-extension conditions. However, has been rarely studying the correlation between intra-plate volcanism and distant rifting. To identify the origin of occurring volcanoes in the continental margin and intra-plate is necessary studying to evolution mechanism and lower mantle. The formation of rifting and continental margin is closely related, and it is assumed that mantle convection significantly affects intra-plate volcanism. Well known through previous studies that mantle convection Along the rifting axis affects various evolution such as rift propagation, ridge jump, rift fail or end tip with transform fault. In this study, we estimate the possibility mantle convection can induce intra-plate volcanism at rifting end tip and continental beyond the margin. We adopt the open-source finite element geodynamics software, ASPECT, which makes a 3-D rifting model for observing the evolution process and mantle convection below the continental margin.

Glacier mass loss related to ice flow is a quantitative factor that controls sea level rise and global warming. Ice flow is attributed to the complicated interaction between grounding line migration, calving basal friction, and topography. Previous studies using massive geophysical observation, including IPR (Ice Penetrating Radar), InSAR (Interferometric Synthetic Aperture Rada), and heat flow exploration, highlight the need to develop a sophisticated numerical model. Because it is difficult to numerically solve the full-stokes equation obtaining ice velocities, the simplified governing equation (i.e., HO, Higher Order model; SSA, Shallow Shelf Approximation; SIA, Shallow Ice Approximation) is widely used in the ice sheet dynamics community. The SSA approach, which assumes a depth-independent velocity model, has computational cost reduction based on simplified descriptions of the full-stokes equation. Here we developed the two-dimensional SSA numerical model to better understand ice dynamics using COMSOL Multiphysics® (hereafter COMSOL), which is a user-friendly finite element software providing convenient GUI, mesh generation, post-processing tools, various types of elements, and the order of shape function. To verify the application of COMSOL to ice sheet dynamics, we compared it with an open-source finite element package, Ice Sheet System Model (ISSM). The concise toy model successfully simulated the distribution of viscosity and velocity and the evolution of surface topography in both COMSOL and ISSM. Furthermore, we applied realistic bathymetry data (i.e., Bedmap2) of Pine Island Glacier, where the large ice mass loss occurs, to clarify more similarity of each model. The surface velocity patterns calculated by COMSOL and ISSM are significantly similar for various physical properties (e.g., ice viscosity, friction, and rate of accumulation and melt). We propose that COMSOL, which efficiently handles mesh generation and visualization, and various weak forms, can sufficiently be applied in ice dynamics.

Subsurface fluid injection of environmental (e.g., carbon capture storage) and industrial (e.g., enhanced geothermal system) projects change pore pressure and underground stress which may induce fault slip. The temperature of injected fluid controls pore pressure and underground stress by thermo-poroelastic effect of fluid injection, that explains the interaction between pore fluid flow and elastic deformation in a porous medium. Since the perturbed subsurface stress distribution increases seismic uncertainty, a numerical modeling for varying injection conditions (e.g., injection scenario and fluid temperature) is helpful for understanding thermo-poroelastic behavior before the fluid injection. In this study, we build 2-dimensional finite element fluid injection models that simulate thermo-poromechanical processes using a COMSOL Multiphysics®. The thermal equation (Fourier’s Law) is coupled with the poroelastic theory to investigate the role of thermohydraulic-convection and -stress in a porous medium. We confirm that these injection conditions may change pore pressure, subsurface stress, and surface displacement, which supports the necessity of monitoring during/after fluid injection.

Marine ice sheet dynamics plays a key role in the mechanical interaction between glaciers and the ocean. A steady-state marine ice sheet, due to the buttressing force of ice shelves, suppresses the glaciers from discharging into the ocean. Various atmospheric-oceanic-glacial mechanisms such as ice melting, buttressing, and snow accumulation lead to the instability of marine ice sheets. The unstable marine ice sheet cause thinning ice thickness and the resultant ice shelf collapse. Ice shelf collapse accelerates the rapid glaciers discharge and retreats, directly affecting sea level rise. The positive feedback loop of these processes significantly influences the ice sheets dynamics, especially the advance and retreat of the grounding line. The grounding line can be defined as a transition zone between the grounded ice sheet in bedrock and the floating ice shelf due to water buoyancy. Understanding grounding line dynamics may help us predict marine ice sheet instability and future sea level rise. It has been challenging to determine an accurate grounding line in the numerical simulation of the marine ice sheet due to mesh resolution. Here we performed two-dimensional ice flow modeling using open-source finite element software (i.e., Elmer/ice), to quantitatively evaluate the effect of diverse mesh sizes on the position of the grounding line. We tested a series of numerical models to precisely define varying the grounding line position influenced by bedrock topography, snow accumulation, and ice melting. The models with high mesh resolution required short time steps to obtain the accurate grounding line position. To consider both calculation efficiency and the accuracy of the position of the grounding line, we found adequate time steps corresponding to each mesh size ambient the grounding line. Our systematic results can provide inspiration for choosing a suitable mesh size and time step to determine a more accurate grounding line position.

Estimating temperature in subduction zones is not easy because the tectonic structures in the area are complicated and it is not possible to reach and directly measure the temperature of subducting slab. In addition, existing water flows in the highly permeable sediment layer interfere with temperature measurements on the seafloor. The observed systematic decays of the seafloor spreading magnetic anomalies in subduction zone give insight into using it as a non-contact geothermometer. We interpret this first decay in the magnetic anomaly amplitude as reflecting the opening of normal faults and fissures from the outer flexural rise, and the renewed circulation of seawater into the oceanic crust, resulting in additional alteration of the magnetic mineral. And the systematic decays are happened in subducted oceanic crust following two steps: (1) the thermal demagnetization of the extrusive basalt layer where titanomagnetite and titanomaghemite reach their blocking temperature between 150-350°C; and (2) the magnetization of the deeper crustal layers slowly decreases, reflecting the progressive heating of the slab toward the Curie temperature of magnetite, 580 °C. These results indicate that the temperature of the subduction zone might be higher than the previous understandings. To identify it, we developed a numerical geodynamic model for subduction zones. Of course, the thermal structure in the subduction zone can be explained in large part by the heat conduction as in the existing geodynamic studies. However, our result support that hydrothermal circulations have a significant affect that changes in the physical properties and thermal structures of subduction zones., The 28th IUGG General Assembly (IUGG2023) (Berlin 2023)

Carbon dioxide emissions from dolomite decarbonation play an essential role in the weakening of carbonate faults by lowering the effective normal stress, which is thermally activated at temperatures above 600–700 °C. However, the mechanochemical effect of low-crystalline ultrafine fault gouge on the decarbonation and slip behavior of dolomite-bearing faults remains unclear. In this study, we obtained a series of artificial dolomite fault gouges with systematically varying particle sizes and dolomite crystallinities using a high-energy ball mill. The laboratory-scale pulverization of dolomite yielded MgO at temperatures below 50 °C, indicating that mechanical decarbonation without significant heating occurred due to the collapse of the crystalline structure, as revealed by X-ray diffraction and solid-state nuclear magnetic resonance results. Furthermore, the onset temperature of thermal decarbonation decreased to ∼400 °C. Numerical modeling reproduced this two-stage decarbonation, where the pore pressure increased due to low-temperature thermal decarbonation, leading to slip weakening on the fault plane even at 400–500 °C; i.e., 200–300 °C lower than previously reported temperatures. Thus, the presence of small amounts of low-crystalline dolomite in a fault plane may lead to a severely reduced shear strength due to thermal decomposition at ∼400 °C with a small slip weakening distance.

The growth mechanism of the Tibetan Plateau, postulated by a number of hypotheses, remains under intense debate. Our analysis of recent satellite-based gravity model reveals that Tibetan lithosphere has been decoupled and folded. It is further evidenced by the existence of crustal melts and channel flow that have been observed by seismic and magnetotelluric explorations. Based on 3D geodynamic simulations, we elucidate the exact buckling structures in the upper crust and lithospheric mantle: at mixed wavelengths between ∼240 and ∼400 km, the lower crustal viscosity is smaller than ∼10 19 Pa·s, implicating weak lower crustal flow beneath the Plateau. This mixed wavelength is consistent with the result of our inverse gravity modeling. Our results facilitate a new plausible hypothesis that the decoupled lithospheric folding mechanism can explain the growth mechanism of the anomalously thick and wide Tibetan Plateau by conflating our idea and contemporary hypothesized scientific findings.

SUMMARY Our understanding of the seismicity of continental interiors, far from plate margins, relies on the ability to account for behaviours across a broad range of time and spatial scales. Deformation rates around seismic faults range from the slip-on-fault during earthquakes to the long-term viscous deformation of surrounding lithosphere, thereby presenting a challenge to modelling techniques. The aim of this study was to test a new method to simulate seismic faults using a continuum approach, reconciling the deformation of viscoelastoplastic lithospheres over geological timescales. A von Mises yield criterion is adopted as a proxy for the frictional shear strength of a fault. In the elastoplastic fault models a rapid change in strength occurs after plastic yielding, to achieve stress–strain equilibrium, when the coseismic slip and slip velocity from the strain-rate response and size of the fault are calculated. The cumulative step-function shape of the slip and temporally partitioned slip velocity of the fault demonstrated self-consistent discrete fault motion. The implementation of elastoplastic faults successfully reproduced the conceptual models of seismic recurrence, that is strictly periodic and time- and slip-predictable. Elastoplastic faults that include a slip velocity strengthening and weakening with reduction of the time-step size during the slip stage generated yield patterns of coseismic stress changes in surrounding areas, which were similar to those calculated from actual earthquakes. A test of fault interaction captured the migration of stress between two faults under different spatial arrangements, reproducing realistic behaviours across time and spatial scales of faults in continental interiors.

Experimental and numerical methods.

To understand back-arc basin dynamics of the western Pacific, constraining the crustal type and rheology of the Ulleung Basin in the East Sea (Sea of Japan) is essential. We performed finite elemen...