Sensitivity of observational data is important in the study of Glacial Isostatic Adjustment (GIA). However, depending on whether sensitivity is used for the Inverse Problem or the Forward Problem, the final formulation and display of the sensitivity kernel will be different. Unfortunately, in the past, both perspectives give the same name to their quantity computed/displayed, and that has caused some confusion. To distinguish between the two, their perspective should be added to the names. This paper focuses only on the perspective of the Forward Problem where the input parameters are known. The Perturbation method has been successfully used in the computation of the sensitivity kernels of observations on 1D and 3D viscosity variations from the Forward perspective. One aim of this paper is to review and clarify the physics of the Perturbation method and bring out some important aspects of this method that have been misunderstood or neglected. Another aim is to present sensitivity kernels from the Perturbation method using 3D (both radially and laterally heterogeneous) Earth models with realistic ice history. These new results are now suitable for future comparison with those from new methods using the Forward perspective. Finally, the sensitivity computations for realistic ice histories on a 3D Earth is reviewed and used to search for optimal locations of new GIA observations.

A glacial forebulge is a bending-related upheaval of the lithosphere outside a glaciated area that co-occurs to the depression of the lithosphere below an ice sheet. The forebulge of the last glaciation attracted attention over more than one century, but quantitative descriptions on the geometry of the forebulge are rare. While many studies mention forebulge dynamics as a possible cause for a certain observation, very few studies provide a detailed and systematic exploration of the forebulge's precise dynamics. That way the forebulge became occasionally a rather mysterious structure with many unknowns. We aim to shed light into the forebulge discussion. After reviewing the history of forebulge research, we outline the theory behind the spatio-temporal forebulge development including controlling factors, and present forebulge observations in geological and geodetic records of North America and the northern parts of Central Europe. We use a state-of-the art finite-element model that can fit multiple observations of the last glaciation simultaneously, to illustrate forebulge development in North America and northern Central Europe and address the issue of whether the zero-uplift hinge line is a good indication of the location of the forebulge front. Finally, we discuss effects of the forebulge on the sea-level change pattern and the evolution of lithospheric stresses, which can induce intraplate earthquakes. We also show that the existence of a glacial forebulge outside the ice margin is not consistent with the assumption of isostatic equilibrium at the Last Glacial Maximum, and there is no strain rate-stress paradox.

Intra-plate faults are a special challenge in seismology, because of the long intervals between individual seismic events and the fact that such faults are often hidden below young sediments. This makes such faults difficult to detect and thus they can be the source of unexpected and fatal earthquakes. The Børglum fault is located in a slowly deforming area in northern Denmark and represents one of the northern boundary faults of the Sorgenfrei-Tornquist Zone. With a length of at least 250 km, it is capable to produce significant seismic events. Previous studies indicated that the Børglum fault is seismically active and this fuelled the demand for further analysis of the fault structure and its seismic hazard potential. Due to excellent coastal outcrops and available high-resolution DEMs, the Børglum fault is a perfect natural laboratory to analyse a hidden active fault. We present a multi-method approach based on outcrop analyses, shear-wave seismic reflection surveys, DEM analysis and numerical simulations of deglaciation-induced Coulomb failure stress change. The 2D seismic surveys show that the analysed segment of the Børglum fault is a complex fault system with a strike-slip component. This interpretation is based on positive flower structures on the seismic surveys, the presence of elongated mini-basins and the geometry of the drainage pattern in the study area. On the basis of soft-sediment deformation structures and disaggregation bands developed in Late Pleniglacial to Lateglacial marine and lacustrine deposits, we derive repeated phases of fault activity with earthquake magnitudes of up to M=7. The geometry of the drainage pattern in the study area indicates a close relationship between fault activity and topography. Based on the timing of fault activity and results from numerical simulations of deglaciation-related lithospheric stress build-up, it is likely that the Børglum fault is a glacially triggered fault and that the analysed part of the Sorgenfrei-Tornquist Zone is susceptible to glacially triggered fault reactivation.

SUMMARYThis paper presents a method that modifies commercial engineering-oriented finite element packages for the modelling of Glacial Isostatic Adjustment (GIA) on a self-gravitating, compressible and spherical Earth with 3-D structures. The approach, called the iterative finite element body and surface force (FEMIBSF) approach, solves the equilibrium equation for deformation using the ABAQUS finite element package and calculates potential perturbation consistently with finite element theory, avoiding the use of spherical harmonics. The key to this approach lies in computing the mean external body forces for each finite element within the Earth and pressure on Earth's surface and core–mantle boundary (CMB). These quantities, which drive the deformation and stress perturbation of GIA but are not included in the equation of motion of commercial finite element packages, are implemented therein. The method also demonstrates how to calculate degree-1 deformation directly in the spatial domain and Earth-load system for GIA models. To validate the FEMIBSF method, loading Love numbers (LLNs) for homogeneous and layered earth models are calculated and compared with three independent GIA methodologies: the normal-mode method, the iterative body force method and the spectral-finite element method. Results show that the FEMIBSF method can accurately reproduce the unstable modes for the homogeneous compressible model and agree reasonably well with the Love number results from other methods. It is found that the accuracy of the FEMIBSF method increases with higher resolution, but a non-conformal mesh should be avoided due to creating the so-called hanging nodes. The role of a potential force at the CMB is also studied and found to only affect the long-wavelength surface potential perturbation and deformation in the viscous time regime. In conclusion, the FEMIBSF method is ready for use in realistic GIA studies, with modelled vertical and horizontal displacement rates in a disc load case showing agreement with other two GIA methods within the uncertainty level of GNSS measurements.

The mid-Holocene sea-level highstand refers to the development of higher-than-present relative sea levels (RSLs) in far-field regions between 7,000 and 4,000 years ago because of equatorial ocean syphoning and continental levering. The timing, magnitude and spatial variability of the highstand are uncertain and the highstand parameterization in Glacial Isostatic Adjustment (GIA) modelling is understudied. Here, we use the RSL records of Southeast Asia to investigate the sensitivity of the mid-Holocene highstand properties to ice and Earth model parameters, including lithospheric thickness, mantle viscosity (both 1D and 3D), and deglaciation history of Antarctica and global ice sheets. We found that the Earth model variation only affects the magnitude of the mid-Holocene highstand unless low upper mantle viscosity is used. The timing of the highstand moves towards present and there is an absence of the highstand if upper mantle viscosity is <4.0 ×1019 Pa s or ≤1.0 ×1019 Pa s, respectively. Ice model variation changes both the timing and magnitude of the mid-Holocene highstand. Delaying the ice-equivalent sea level will shift the timing of the highstand later and result in a lower highstand magnitude. We produced a mid-Holocene highstand “treasure map” that considers topography change and accommodation space to guide future RSL data collection efforts in Southeast Asia. The highstand “treasure map” indicates the northern east coast and central west coast of Malay-Thai Peninsula, east coast of Sumatra, north coast of Java, and southwest coast of Borneo are very likely (90% probability) to preserve mid-Holocene RSL data.

The western Russian Arctic was partially covered by the Eurasian ice sheet complex during the Last Glacial Maximum (~26 ka BP) and is a focus area for Glacial Isostatic Adjustment (GIA) studies. However, there have been few GIA studies conducted in the Russian Arctic due to the lack of high quality deglacial relative sea-level (RSL) data. Recently, Baranskaya et al. (2018) released a quality-controlled deglacial RSL database for the Russian Arctic that consists of ~400 sea-level index points and ~250 marine and terrestrial limiting data that constrain RSL since 20 ka BP. Here, we use the RSL database to constrain the 3D Earth structure beneath the Russian Arctic, with consideration of the uncertainty in ice model ICE-7G_NA, which is assessed via iteratively refining the ice model with fixed 1D Earth model to achieve a best fit with the RSL data. Also, the uncertainties in 3D Earth parameters and RSL predictions are investigated. We find an optimal 3D Earth model (Vis3D) improves the fit with the deglacial RSL data compared with the VM7 1D model when fixed with the ICE-7G_NA ice model. Similarly, we show improved fit in the White Sea area, where 1D model shows notable misfits, with the refined ice model ICE-7G_WSR when fixed with VM7 Earth model. The comparable fits of ICE-7G_NA (Vis3D) and ICE-7G_WSR (VM7) implies that the uncertainty in the ice model might be improperly mapped into 3D viscosity structure when a fixed ice model is employed. Furthermore, fixed with refined ice model ICE-7G_WSR, we find an optimal 3D Earth model (Vis3D_R), which fits better than ICE-7G_WSR (VM7), and the magnitude of lateral heterogeneity decreases significantly from Vis3D to Vis3D_R. We conclude that uncertainty in the ice model needs to be considered in 3D GIA studies.

GRACE-based estimates for groundwater storage (GWS) changes in North America substantially depend upon correction of glacial isostatic adjustment (GIA) effects, which are usually removed with GIA models. In this study, GIA effects are eliminated by employing an independent separation approach with the aid of Global Navigation Satellite System (GNSS) vertical velocity data. Our goal is to provide an independent estimate for monthly GWS changes within North America in 1-degree-grids and their trends over the whole GRACE mission lifetime from April 2002 to June 2017. This estimate is derived from the release-6 version of GRACE monthly level-2 data, GNSS data, land surface models for soil moisture and snow water equivalent, and satellite altimetric lake level data. We find a GWS anomaly in form of an increasing trend in Saskatchewan, which affects the Saskatchewan Province and the states of Montana, North Dakota and Minnesota, and 4 GWS anomalies with declining trends in Nevada, California, Arizona and Texas, respectively. The monthly changes of these GWS anomalies, except for the one in Nevada, are validated by well level data. We provide results for average monthly GWS changes and the trends for the 5 anomalies but also in separate form for the 13 affected states or provinces. The increasing trends of the Saskatchewan GWS anomaly and the affected 3 states are related to increasing precipitation and can be elucidated by the decreasing drought intensity level. On the contrary, the declining trends in GWS can be explained by weakening precipitation and are mostly supported by the increasing drought intensity level in the other 4 anomalies and the affected states, which are Nevada, California, Arizona, New Mexico, Texas, Oklahoma, Kansas, and Colorado. Our estimates of monthly GWS changes and their trends can serve as alternative and beneficial input for the sustainable management of groundwater resources in North America.

Analyses of glacial isostatic adjustment (GIA) and deglacial relative sea‐level (RSL) change in the Russian Arctic deliver important insights into the Earth's viscosity structure and the deglaciation history of the Eurasian ice sheet complex. Here, we validate the 1D GIA models ICE‐6G_C (VM5a) and ICE‐7G_NA (VM7) and select new 3D GIA models in the Russian Arctic against a quality‐controlled deglacial RSL database of >500 sea‐level data points from 24 regions. Both 1D models correspond to the RSL data along the southern coast of the Barents Sea and Franz Josef Land from ∼11 ka BP to present but show notable misfits (>50 m at 10 ka BP) with the White Sea data. We find 3D model predictions of deglacial RSL resolve most of the misfits with the observed data for the White Sea while retaining comparable fits in other regions of the Russian Arctic. Our results further reveal: (a) RSL in the western Russian Arctic is sensitive to elastic lithosphere with lateral thickness variation and 3D viscosity structure in the upper mantle; and (b) RSL in the whole Russian Arctic is less sensitive to 3D viscosity structure in the lower mantle compared to the upper mantle. The 3D models reveal a compromise in the upper mantle between the background viscosity and scaling factor to best fit the RSL data, which needs to be considered in future 3D GIA studies.

By far the most prescient insights into the interior structure of the planet have been provided on the basis of elastic wave seismology. Analysis of the travel times of shear or compression wave phases excited by individual earthquakes, or through analysis of the elastic gravitational free oscillations that individual earthquakes of sufficiently large magnitude may excite, has been the central focus of Earth physics research for more than a century. Unfortunately, data provide no information that is directly relevant to understanding the solid state ‘flow’ of the polycrystalline outer ‘mantle’ shell of the planet that is involved in the thermally driven convective circulation that is responsible for powering the ‘drift’ of the continents and which controls the rate of planetary cooling on long timescales. For this reason, there has been an increasing focus on the understanding of physical phenomenology that is unambiguously associated with mantle flow processes that are distinct from those directly associated with the convective circulation itself. This paper reviews the past many decades of work that has been invested in understanding the most important of such processes, namely that which has come to be referred to as ‘glacial isostatic adjustment’ (GIA). This process concerns the response of the planet to the loading and unloading of the high latitude continents by the massive accumulations of glacial ice that have occurred with almost metronomic regularity over the most recent million years of Earth history. Forced by the impact of gravitational n -body effects on the geometry of Earth’s orbit around the Sun through the impact upon the terrestrial regime of received solar insolation, these surface mass loads on the continents have left indelible records of their occurrence in the ‘Earth system’ consisting of the oceans, continents, and the great polar ice sheets on Greenland and Antarctica themselves. Although this ice-age phenomenology has been clearly recognized since early in the last century, it was for over 50 years considered to be no more than an interesting curiosity, the understanding of which remained on the periphery of the theoretical physics of the Earth. This was the case in part because no globally applicable theory was available that could be applied to rigorously interpret the observations. Equally important to understanding the scientific lethargy that held back the understanding of this phenomenon involving mantle flow processes was the lack of appreciation of the wide range of observations that were in fact related to GIA physics. This paper is devoted to a review of the global theories of the GIA process that have since been developed as a means of interpreting the extensive variety of observations that are now recognized as being involved in the response of the planet to the loading and unloading of its surface by glacial ice. The paper will also provide examples of the further analyses of Earth physics and climate related processes that applications of the modern theoretical structures have enabled.

Glacial Isostatic Adjustment (GIA) induced by the melting of the Pleistocene Ice Sheets causes differential land uplift, relative sea level and geoid changes. Thus, GIA in North America may affect water flow-accumulation and the rate of sedimentation and erosion in the South Saskatchewan River Basin (SSRB), but so far this has not been well investigated. Our aim here is to use surface topography in the SSRB and simple models of surface water flow to compute flow-accumulation, wetness index, stream power index and sediment transport index - the latter two affect the rates of erosion and sedimentation. Since the river basin became virtually ice-free around 8 ka BP, we shall study the effects of GIA induced differential land uplift during the last 8 ka on these indexes. Using the present-day surface topography ETOPO1 model, we see that the stream power index and sediment transport index in the SSRB may not be high enough to alter the surface topography significantly today and probably during the last 8 ka except for places around the Rocky Mountains. The effect of using 1 and 3 arc minute grid resolution of the ETOPO1 model does not significantly alter the value of these indexes. However, we note that using 1 arc minute grid is much more computationally intensive, so only a smaller area of the SSRB can be included in the computation. Next, we assume that sedimentation and erosion did not occur in the SSRB during the last 8 ka BP, and the change in surface topography is only due to GIA induced differential uplift. We use land uplift predicted by a large number of GIA models to study the changes in stream power & sediment transport indexes in the last 8 ka BP. Our base GIA model is ICE6G_C(VM5a). Then we investigate the effects of using uplift predicted by other GIA models that can still fit the observed relative sea level (RSL), uplift rate and gravity-rate-of-change data in North America reasonably well. These alternate GIA models have lateral heterogeneity in the mantle and lithosphere included – in particular we test those that give the largest differential uplift in the SSRB. We found that the effect of these other GIA earth models is not large on the stream power & sediment transport indexes. Finally, we investigate the sensitivity of these indexes on the ice models that are consistent with GIA observations. The results of this study will be useful to our understanding of water flow accumulation, sedimentation and erosion in the past, present and future and for water resource management in North America.

Holocene relative sea-level (RSL) records from regions distal from ice sheets (far-field) are commonly characterized by a mid-Holocene highstand, when RSL reached higher than present levels. The magnitude and timing of the mid-Holocene highstand varies spatially due to hydro-isostatic processes including ocean syphoning and continental levering. While there are open questions regarding the timing, magnitude and source of ice-equivalent sea level in the middle to late Holocene. Here, we compare Glacial Isostatic Adjustment (GIA) model predictions to a standardized database of sea-level index points (SLIPs) from Southeast Asia where we have near-complete Holocene records. The database has more than 130 SLIPs that span the time period from ~9.5 ka BP to present. We investigate the sensitivity of mid-Holocene RSL predictions to GIA parameters, including the lateral lithospheric thickness variation, mantle viscosity (both 1D and 3D), and deglaciation history from different ice sheets (e.g., Laurentide, Fennoscandia, Antarctica). We compute gravitationally self-consistent RSL histories for the GIA model with time dependent coastlines and rotational feedback using the Coupled Laplace-Finite Element Method. The preliminary results show that the timing of the highstand is mainly controlled by the deglaciation history (ice-equivalent sea level), while the magnitude is dominated by Earth parameters (e.g., lithospheric thickness, mantle viscosity). We further investigate whether there is meltwater input during middle to late Holocene and whether the RSL records from Southeast Asia can reveal the meltwater source, like Antarctica.

The reactivation of glacially induced faults is linked to the increase and decrease of ice mass. But, whether faults are reactivated by glacially induced stresses depends to a large degree on the crustal stress field, fault properties and fluid pressures. The background (tectonic and lithostatic) stress field has a major effect on the potential for reactivation, as the varying stresses induced by the ice sheet affects the state of stress around the fault, bringing the fault to more stable or more unstable conditions. Here, we describe the effect of glacially induced stresses on fault reactivation under three potential background stress regimes of normal, strike-slip and thrust/reverse faulting. The Mohr diagram is used to illustrate how glacially induced stresses affect the location and the size of the Mohr circle. We review these different cases by applying an analysis of the stress state at different time points in the glacial cycle. In addition, we present an overview of fault properties that affect the reactivation of glacially induced faults, such as pore-fluid pressure and coefficient of friction.

To model glacial triggering of earthquakes, it is necessary to obtain the spatio-temporal variation of glacial isostatic adjustment-induced stress during a glacial cycled. This can be computed efficiently using commercial Finite Element codes with appropriate modifications to include the important effects of 'pre-stress advection', 'internal buoyancy' and 'self-gravity'. The modifications described in Wu (2004) are reviewed for incompressible and so-called materially compressible flat-earths. When the glacial isostatic adjustment-induced stress is superimposed on the background tectonic stress and overburden pressure, the time variation of earthquake potential at various locations in the Earth can be evaluated for any fault orientation. To model more complex slip and fault behavior over time, the three-stage Finite Element model approach of Steffen et al. (2014) is reviewed. Finally, selected numerical examples and their results from both modelling approaches are shown.

We quantify GIA prediction uncertainties of 250 1D and 3D glacial isostatic adjustment (GIA) models through comparisons with deglacial relative sea-level (RSL) data from North America and rate of vertical land motion (<(U)over dot>) and gravity rate of change (<(G)over dot>) from GNSS and GRACE data, respectively. Spatially, the size of the RSL uncertainties varies across North America with the largest from Hudson Bay and near previous ice margins along the northern Atlantic and Pacific coasts, which suggests 3D viscosity structure in the lower mantle and laterally varying lithospheric thickness. Temporally, RSL uncertainties decrease from the Last Glacial Maximum to present except for west of Hudson Bay and the northeastern Pacific coast. The uncertainties of both these regions increase from 30 to 45 m between 15 and 11 ka BP, whichmay be due to the rapid decrease of surface loading at that time. Present-day <(U)over dot> and <(G)over dot>uncertainties are largest in southwestern Hudson Bay with magnitudes of 2.4 mm/year and 0.4 mu Gal/year, mainly due to the 3D viscosity structure in the lower mantle.

The Canadian landmass of North America and the Russian Arctic were covered by large ice sheets during the Last Glacial Maximum, and have been key areas for Glacial Isostatic Adjustment (GIA) studies. Previous GIA studies have applied 1D models of Earth’s interior viscoelastic structure; however, seismic tomography, field geology and recent studies reveal the potential importance of 3D models of this structure. Here, using the latest quality-controlled deglacial sea-level databases from North America and the Russian Arctic, we investigate the effects of 3D structure on GIA predictions. We explore scaling factors in the upper mantle (βUM) and lower mantle (βLM) and the 1D background viscosity model (ηo) with predictions of of the ICE-6G_C (VM5a) glaciation/deglaciation model of Peltier et al (2015, JGR) in these two regions, and compare with the best fit 3D viscosity structures.We compute gravitationally self-consistent relative sea-level histories with time dependent coastlines and rotational feedback using both the Normal Mode Method and Coupled Laplace-Finite Element Method. A subset of 3D GIA models is found that can fit the deglacial sea-level databases for both regions. These databases cover both the near and intermediate field regions. However, North America and Russian Arctic prefer different 3D structures (i.e., combinations of (ηo, βUM, βLM)) to provide the best fits. The Russian Arctic database prefers a softer background viscosity model (ηo), but larger scaling factors (βUM, βLM) than those preferred by the North America database.Outstanding issues include the uncertainty of the history of local glaciation history. For example, preliminary modifications of the ice model in Russian Arctic reveal that the misfits of 1D models can be significantly reduced, but still fit less well than the best fit 3D GIA model.An additional issue concerns the extent to which the 3D models are able to improve both fits in North America and Russian Arctic when compared with 1D internal structure (ICE-6G_C VM5a & ICE-7G VM7), will be assessed in a preliminary fashion.

We study the terrestrial water storage (TWS) and groundwater storage (GWS) changes in Canada and United States. We employ the separation approach from Wang et al. (2013) together with the improved GRACE data of Release 6 for a longer time span until December, 2016. The TWS signals from lake levels are derived from satellite altimetry data over the lakes while TWS signals due to soil moisture (SM) and snow water equivalent (SWE) changes from hydrology models. There are four significant trend anomalies in North America for both TWS and GWS changes. Two positive anomalies are found in Canada with their centers in the provinces of Saskatchewan and Quebec, respectively, due to increased precipitation and/or increased runoff in their surroundings. Two negative anomalies are shown in the United States with their centers in California and the northwest of Texas, respectively, which are due to decreased precipitation and, especially for California, high water usage for agriculture.

Global reconstructions of the Earth's Quaternary ice sheet history such as ICE-6G assume that mantle rheology is linear and Earth properties are laterally homogeneous. However, these assumptions are unrealistic because high temperature/pressure creep experiments show that both linear and non-linear creep laws operate simultaneously in the mantle. Also surface geology and seismic tomography show that mantle properties vary laterally. This study uses the Coupled Laplace-Finite Element Method to find an ice model that is consistent with an Earth model with composite rheology where both linear and non-linear creep laws operate simultaneously in the mantle, and the effective viscosity is laterally heterogeneous. We use the ICE-6G ice model as initial load on a composite rheology Earth model and iteratively improve both the ice and Earth model parameters until they fit observations of Glacial Isostatic Adjustment or GIA (relative sea level and land uplift rate) well simultaneously. We thereby focus on North America and northern Europe. Our best combination model, called ICE-C(M_C), also fits gravity rate-of-change data in both areas. This finding removes the obstacle in using composite rheology in GIA studies since previous work found it impossible for composite rheology to fit both relative sea level and land uplift rate data simultaneously. It shows that composite rheology is clearly an alternative to linear rheology in future GIA studies. The ICE-C model is at the last glacial maximum (LGM) position from 26 ka BP (thousand years before present) to 19 ka BP, which implies that composite rheology also supports the timing of LGM proposed by Clark et al. (2009). The model is found to be significantly thicker in the center of rebound in Fennoscandia than ICE-6G from the LGM to about 8 ka BP, but significantly thinner around the Fennoscandian coastline from LGM to about 12 ka BP. These changes in the lower boundary conditions in northern Europe may affect the reconstruction of past climate there. Regional differences in ice thickness between ICE-C and ICE-6G are also found in Laurentia. During LGM there are negligible differences in total ice volume between ICE-C and ICE-6G, thus it remains questionable if a composite rheology may help in solving the missing ice problem. ICE-C is also consistent with meltwater pulses 1A and 1B. Mantle viscosities and creep parameters of Earth model M_C are consistent with previous findings from GIA modeling and microphysics experiments. (C) 2019 Elsevier B.V. All rights reserved.

A new generation of numerical models is being developed to model Glacial Isostatic Adjustment in a self-gravitating spherical earth with lateral heterogeneity and/or nonlinear rheology in the mantle. Of special interest is the Iterative Stress Transform (IST) method (also known as Coupled Laplace-Finite Element method) because it uses commercial finite-element packages which are readily available, well tested and reliable. Although the IST method is efficient and can produce very accurate results, it is mainly developed for incompressible earth models. So there are efforts to generalize the IST method for more realistic compressible earths. Here, we will extend the finding of Bangtsson & Lund to confirm that the IST method is not likely to be generalized for a compressible earth. Next, a new approach, called the Iterative Body Force (IBF) is presented, which aims to replace all body forces in each element by their volumetric average and solve the governing differential equations iteratively. The result of the IBF approach is then benchmarked with the conventional normal mode method (NMM) for laterally homogeneous axisymmetric earth models forced by Heaviside harmonic loads. For incompressible earth models, the IBF approach gives excellent agreement with NMM that uses analytical propagation method. For compressible earth models, good agreement is also obtained with NMM that uses the numerical integration method, provided that the timescale of compressional instability is long compared with the loading period. However, the agreement deteriorates if the effect of gravitational instability is significant because the numerical errors grow differently for each method. Finally, the IBF approach is used to study the spatiotemporal evolution of body forces and understand the development of instability. It is shown that compression in a uniform layer can result in the top becoming denser than the bottom of the same layer which promotes Rayleigh-Taylor instability. If this can be suppressed by the stabilization forces of pre-stress advection and internal buoyancy of the layer, then stability remains. However, if the compressional instability becomes large enough to change the direction of the pre-stress advection force, then the deformation can grow so large that convection instability can be triggered.