ALBERT

Description: Glacial isostatic adjustment is largely governed by the rheological properties of the Earth's mantle. Large mass redistributions in the ocean–cryosphere system and the subsequent response of the viscoelastic Earth have led to dramatic sea level changes in the past. This process is ongoing, and in order to understand and predict current and future sea level changes, the knowledge of mantle properties such as viscosity is essential. In this study, we present a method to obtain estimates of mantle viscosities by the assimilation of relative sea level rates of change into a viscoelastic model of the lithosphere and mantle. We set up a particle filter with probabilistic resampling. In an identical twin experiment, we show that mantle viscosities can be recovered in a glacial isostatic adjustment model of a simple three-layer Earth structure consisting of an elastic lithosphere and two mantle layers of different viscosity. We investigate the ensemble behaviour on different parameters in the following three set-ups: (1) global observations data set since last glacial maximum with different ensemble initialisations and observation uncertainties, (2) regional observations from Fennoscandia or Laurentide/Greenland only, and (3) limiting the observation period to 10 ka until the present. We show that the recovery is successful in all cases if the target parameter values are properly sampled by the initial ensemble probability distribution. This even includes cases in which the target viscosity values are located far in the tail of the initial ensemble probability distribution. Experiments show that the method is successful if enough near-field observations are available. This makes it work best for a period after substantial deglaciation until the present when the number of sea level indicators is relatively high.

Description: Viscoelastic deformations of an earth structure in response to a time-varying surface load are analyzed in glacial isostatic adjustment (GIA). When solving this problem, aspects like flexure of the lithosphere and retarded response of mantle material become evident. Quantified are these by flexural rigidity and relaxation times. The concepts partly lose their relevance when changing from a 1D earth structure (only radial variations) to a 2D or a 3D earth structure (lateral variations). In regions like Fennoscandia and Laurentide, which are affected by GIA, lateral variations of the lithosphere and mantle structure are moderate and, so, the application of a 1D earth structure is widely accepted. But, also for these two regions one has to keep in mind that the respective 1D earth structures differ and that such an approximation mainly holds in the central part of the respective region. In contrast, lateral variations or a local structure of different viscosity have to be considered in areas like Patagonia, Antarctica or Alaska which is located above tectonic activity or covers a region with significant lateral changes in earth structure. But, already for the two former examples one has to keep in mind that the respective 1D earth structures inferred from GIA modelling differ between the two regions. Focusing on the relaxation behavior and the mantle-material transport, we discuss the effect of lateral variations on the deformation process. We will assess to which extent a 1D earth structure can represent lateral variability in structural features, and, at which point a 3D earth structure has to be considered. Such questions are of concern, when discussing GIA for geodetic applications as well as in earth system modeling as this study contributes to the climate modeling initiative Palmod.

Description: Glacial-isostatic adjustment (GIA) models simulate the viscoelastic response of the solid earth due to loading. During the last glacial maximum, large areas in the northern and southern hemisphere were covered by km-thick ice sheets. Although most of the ice has been melted already 8,000 year ago, the time-delayed response of the viscoelastic earth is still a significant contribution to present-day uplift rates. The implementation of GIA models in global climate models is an essential part of the current research. Hereby, the choice of an appropriate earth structure in the GIA model plays an important role and has to be constrained by observational data. Here, we apply present-day uplift data to constrain a set of GIA models that differ in 3D earth structure. To this end, these different GIA models are validated against GPS uplift rates provided by Schumacher et al. (2019). The GPS stations are globally distributed and not necessarily clustered in regions with strong GIA signal. For validation, regions with the largest gradient present in the GIA signal are most crucial. Thus, we use a weighting scheme, where those GPS stations get a higher weight that are less correlated to all other stations. Additionally, uncertainties in the GPS rates appear due to the length of the GPS time series and due to station specifics such as the used GPS receiver, and are provided together with the rates as standard deviations. Thence, the weighting used for the validation is the sum of the correlation derived weights and the uncertainty derived weights. With this weighting in place, different GIA models can be validated against present day uplift rates by means of root mean square errors or mean absolute error.

Description: The interaction between ice sheets and the solid Earth plays an important role for ice-sheet stability and sea-level change and hence for global climate models. Glacial-isostatic adjustment (GIA) models enable simulation of the solid Earth response due to variations in ice-sheet and ocean loading and prediction of the relative sea-level change. Because the viscoelastic response of the solid Earth depends on both ice-sheet distribution and the Earth’s rheology, independent constraints for the Earth structure in GIA models are beneficial. Seismic tomography models facilitate insights into the Earth’s interior, revealing lateral variability of the mantle viscosity that allows studying its relevance in GIA modeling. Especially, in regions of low mantle viscosity, the predicted surface deformations generated with such 3D GIA models differ considerably from those generated by traditional GIA models with radially symmetric structures. But also, the conversion from seismic velocity variations to viscosity is affected by a set of uncertainties. Here, we apply geodynamically constrained 3D Earth structures. We analyze the impact of conversion parameters (reduction factor in Arrhenius law and radial viscosity profile) on relative sea-level predictions. Furthermore, we focus on exemplary low-viscosity regions like the Cascadian subduction zone and southern Patagonia, which coincide with significant ice-mass changes.

Description: The Antarctic Ice Sheet rests on a bed that is characterized by tectonical activity and hence by a heterogeneous rheology. Spots of extremely weak lithosphere structure could have strong impacts on the Glacial Isostatic Adjustment and hence on the stability of the ice sheet, possibly also for confined glacier regions and on timescales of decades down to even years (Barletta et al., 2018). We coupled the VIscoelastic Lithosphere and MAntle model (VILMA) to the Parallel Ice Sheet Model (PISM) and ran simulations over the last two glacial cycles. In this framework, VILMA considers both viscoelastic deformations of the solid Earth and gravitationally consistent mass redistribution in the ocean by solving for the sea-level equation (Martinec et al., 2018). In turn, PISM interprets this as a vertical shift in bed topography that directly affects the stress balance within the ice sheet and hence the grounding line dynamics at the interface of ice, ocean and bedrock. Here we present first results of the coupled Antarctic glacial-cycle simulations and investigate technical aspects, such as optimal coupling time steps, iteration schemes and convergence, for both one-dimensional and three-dimensional Earth structures. This project is part of the German Climate Modeling Initiative, PalMod2.

Description: We suggest to apply data assimilation in glacial isostatic adjustment (GIA) to constrain the mantle viscosity structure based on sea level observations. We apply the Parallel Data Assimilation Framework (PDAF) to assimilate sea level data into the time-domain spectral-finite element code VILMA in order to obtain better estimates of the mantle viscosity structure. In a first step, we reduce to a spherically symmetric earth structure and prescribe the glaciation history. A particle filter is used to propagate an ensemble of models in time. At epochs when observations are available, each particle's performance is estimated and the particles are resampled based on their performance to form a new ensemble that better resembles the true viscosity distribution. Using this algorithm, we show the ability to recover mantle viscosities from a set of synthetic relative sea level observations. Those synthetic observations are obtained from a reference run with a given viscosity structure that defines the target viscosity values in our experiments. The viscosity estimation is applied to a three-layer model with an elastic lithosphere and two mantle layers, and to a multi-layer model with a smoother viscosity profile. We use various subsets of realistic observation locations (e.g. only observations from Fennoscandia) and show that it is possible to obtain the target viscosity values in those cases. We also vary the time from which observations are available to evolve the test cases towards a realistic scenario for the availability of relative sea level observations. The most relevant cases start at 26.5ka BP and at 10ka BP as they mark the beginning of the maximum glaciation and the end of deglaciation with a larger amount of observations following, respectively, and end at present day.

Description: The analysis of the Coulomb stress changes has become an important tool for seismic hazard evaluation because such stress changes may trigger or delay next earthquakes. Processes that can cause significant Coulomb stress changes include coseismic slip, earthquake-induced poroelastic effects as well as transient postseismic processes such as viscoelastic relaxation. In this study, we investigate the spatial and temporal evolution of pore fluid pressure changes and fluid flow during the seismic cycle, their dependency on the permeability in the crust and the interaction with postseismic viscoelastic relaxation. To achieve this, we use 2D finite-element models for intra-continental normal and thrust faults, which include coseismic slip, poroelastic effects, postseismic viscoelastic relaxation and interseismic stress accumulation. In different experiments, we vary (1) the permeability of the upper and lower crust while keeping the viscosity structure constant and (2) the viscosity of the lower crust and lithospheric mantle, while we keep the permeabilities constant. (1) The modelling results show that the highest changes in pore fluid pressure during and after the earthquake occur within a distance of ~ 1 km around the lower fault tip at the transition between upper and lower crust. The evolution of pore pressure and fluid flow depends primarily on the permeability in the upper crust. With decreasing permeability, the possibility of the pore fluids to flow decreases and thus, in the postseismic phase, the duration of the poroelastic relaxation increases, from a few days to several years, until the pore pressure reaches the initial pressure of the preseismic phase. In contrast, the influence of variations of the permeability in the lower crust on the pore pressure changes is negligible. For high upper-crustal permeabilities, postseismic vertical velocities are high and decreases rapidly with time, from around 120 mm/a after the first year by two orders of magnitude after 10 years, whereas for low permeabilities they remain consistently low over the years after the earthquake. (2) Models with low viscosity of the lower crust show that the timescales of poroelastic effects and viscoelastic relaxation overlap and affect the postseismic velocity already in the early postseismic phase and that both processes decay within a few years after the earthquake. For higher viscosities, the velocity is initially dominated by pore pressure changes during the first few years, whereas viscoelastic relaxation lasts for decades. Both processes also show differences in their spatial scale. Poroelastic effects occur within a few kilometers around the fault, whereas viscoelastic relaxation acts on tens to hundreds of kilometers. As both processes can cause Coulomb stress changes on faults in the vicinity of the earthquake source fault, it is important to understand the spatial and temporal evolution, the effects on the individual faults and the interaction of both processes during the earthquake cycle. Future work will therefore include the calculation and examination of Coulomb stress changes on intra-continental normal and thrust faults using 3D models that include poroelastic effects and viscoelastic relaxation.

Description: Earthquakes on faults in the brittle upper crust evoke sudden changes in pore fluid pressure as well as postseismic viscoelastic flow in the lower crust and lithospheric mantle but the relative importance of these processes during the postseismic phase has not been systematically studied. Here, we use two-dimensional finite-element models to investigate how pore fluid pressure changes and postseismic viscoelastic relaxation interact during the earthquake cycle of an intracontinental dip-slip fault. To isolate the effects from pore fluid flow and viscoelastic relaxation from each other, we performed experiments with and without pore fluid flow and viscoelastic relaxation, respectively. In different experiments, we further varied the permeability of the crust and the viscosity of lower crust or lithospheric mantle. Our model results show poroelastic effects dominate the velocity field in the first months after the earthquake. In models considering poroelastic effects, the surfaces of both hanging wall and footwall of the normal fault subside at different velocities, while they move upwards in the thrust fault model. Depending on the permeability and viscosity values, viscoelastic relaxation dominates the velocity field from about the second postseismic year onward although poroelastic effects may still occur if the permeability of the upper crust is sufficiently low. With respect to the spatial scales of poroelastic effects and viscoelastic relaxation, our results show that pore fluid pressure changes affect the velocity field mostly within 10–20 km around the fault, whereas the signal from viscoelastic relaxation is recognizable up to several tens of kilometres away from the fault. Our findings reveal that both poroelastic effects and viscoelastic relaxation may overlap earlier and over longer time periods than previously thought, which should be considered when interpreting aftershock distributions, postseismic Coulomb stress changes and surface displacements.

Description: This is a synthetic dataset. It was created from the outputs of the glacial isostatic adjustment model VILMA (Klemann et al. 2008). It consists of realtive sea level (RSL) data on a global regular grid. The resolution is 256 x 512 points (Lat x Lon). The tomporal range is from 123 ka BP until present day. Time steps vary between 2.5 kyrs at the beginning and 0.5 kyrs towards the end. The data were created for a specific configuration of the GIA model: lithosphere thickness = 60 km, lithosphere viscosity = 1.0E31 Pa s, upper mantle thickness = 610 km, upper mantle viscosity = 1.0E20 Pa s, lower mantle thickness = 3,221 km, lower mantle viscosity = 1.0E21 Pa s. The RSL data are accompanied by a observation locations mask. This mask was used to identify those locations in the global RSL dataset where real observations are available. The dataset consists of realtive sea level (RSL) data on a global regular grid. The resolution is 256 x 512 points (Lat x Lon). The temporal range is from 123 ka BP until present day. Time steps vary between 2.5 kyrs at the beginning and 0.5 kyrs towards the end. The data were created for a specific configuration of the GIA model: lithosphere thickness = 60 km, lithosphere viscosity = 1.0E31 Pa s, upper mantle thickness = 610 km, upper mantle viscosity = 1.0E20 Pa s, lower mantle thickness = 3,221 km, lower mantle viscosity = 1.0E21 Pa s. The RSL data are accompanied by observation locations masks. These masks were used to mark those locations in the global RSL dataset where real-life observations are available in order to restrict usage of the synthetic data to those locations.