JGR–Solid Earth

The quartz α → β transition is a displacive phase transition associated with a significant change in elastic properties. However, the elastic properties of quartz at high-pressure and temperature remain poorly constrained experimentally, particularly within the field of β-quartz. Here, we conducted an experimental study on the quartz α → β transition during which P-wave velocities were measured in-situ at pressure (from 0.5 to 1.25 GPa) and temperature (200–900°C) conditions of the continental lower crust. Experiments were carried out on samples of microcrystalline material (grain size of 3–6 μm) and single-crystals. In all these, the transition was observed as a minimum in P-wave velocities, preceded by an important softening while P-wave velocities measured in the β-quartz field were systematically lower than that predicted by thermodynamic databases. Additional experiments during which acoustic emission (AE) were monitored showed no significant peak of AEs near or at the transition temperature. Microstructural analysis nevertheless revealed the importance of microcracking while Electron Back-Scatter Diffraction (EBSD) imaging on polycrystalline samples revealed a prevalence of Dauphiné twinning in samples that underwent the transition. Our results suggest that the velocity change due to the transition known at low pressure might be less important at higher pressure, implying a change in the relative compressibilities of α and β quartz. If true, the velocity changes related to the α → β quartz transition at lower crustal conditions might be lower than that expected in thickened continental crust.

We analyze nearest-neighbor proximities of earthquakes in California based on the joint distribution (T, R) of rescaled time T and rescaled distance R between pairs of earthquakes (Zaliapin & Ben-Zion, 2013a, https://doi.org/10.1002/jgrb.50179), using seismic catalogs from several regions and several catalogs for the San Jacinto Fault Zone (SJFZ). The study aims to identify informative modes in nearest-neighbor diagrams beyond the general background and clustered modes, and to assess seismic catalogs derived by different methods. The results show that earthquake clusters with large and small-to-medium mainshocks have approximately diagonal and horizontal (T, R) distributions of the clustered mode, respectively, reflecting different triggering distances of mainshocks. Earthquakes in the creeping section of San Andreas Fault have a distinct “repeaters mode” characterized by very large rescaled times T and very small rescaled distances R, due to nearly identical locations of repeating events. Induced seismicity in the Geysers and Coso geothermal fields follow mostly the background mode, but with larger rescaled times T and smaller rescaled distances R compared to tectonic background seismicity. We also document differences in (T, R) distributions of catalogs constructed by different techniques (analyst-picks, template-matching and deep-learning) for the SJFZ, and detect a mode with very large R and small T in the template-matching and deep-learning based catalogs. This mode may reflect dynamic triggering by passing waves and/or catalog artifacts.

Weak serpentine minerals affect the mechanical behavior of serpentinized peridotites at depth, and may play a significant role in deformation localization within subduction zones, at local or regional scale. Mixtures of olivine with 5, 10, 20 and 50 vol. % fraction of antigorite, proxies for serpentinized peridotites, were deformed in axial shortening geometry under high pressures (ca. 2–5 GPa) and moderate temperatures (ca. 350°C), with in situ stress and strain measurements using synchrotron X-rays. We evaluate the average partitioning of stresses at the grains scale within each phase (mineral) of the aggregate and compare with pure olivine aggregates in the same conditions. The in situ stress balance is different between low antigorite contents up to 10 vol. %, and higher contents above 20 vol. %. Microstructure and stress levels suggest the deformation mechanisms under these experimental conditions are akin to (semi)brittle and frictional processes. Unlike when close to dehydration temperatures, hardening of the aggregate is observed at low serpentine fractions, due to an increase in local stress concentrations. Below and above the 10–20 vol. % threshold, the stress state in the aggregate corresponds to friction laws already measured for pure olivine aggregates and pure antigorite aggregates respectively. As expected, the behavior of the two-phase aggregate does not evolve as calculated from simple iso-stress or iso-strain bounds, and calls for more advanced physical models of two-phase mixtures.

The past decades have witnessed the increased applications of induced polarization (IP) method in the critical zone studies with ubiquitous clay minerals. Although IP outperforms traditional electrical and electromagnetic methods through its unique ability to measure quadrature conductivity, the nonlinearity that quadrature conductivity behaves with salinities and frequencies greatly tortures IP practitioners, as (a) salinity-dependency makes the quadrature conductivity a varyingly unstable parameter to quantitatively estimate hydraulic properties and clay content; (b) frequency-dependent Cole-Cole and Debye/Warburg decomposition models, although mathematically sound, physically mingle the properties of pore water and clay minerals and are empirical in nature. From basic principles, we demonstrate that quadrature conductivity remains a hybrid property involving both clay and water, and develop relevant models to distinguish them. Our models are validated by theories, experiments, simulations, and comparisons, all of which proclaim considerable advantages over previous models and offer the prospect of quantitative applications.

Fault-damage zones comprise multiscale fracture networks that may slip dynamically and interact with the main fault during earthquake rupture. Using 3D dynamic rupture simulations and scale-dependent fracture energy, we examine dynamic interactions of more than 800 intersecting multiscale fractures surrounding a listric fault, emulating a major listric fault and its damage zone. We investigate 10 distinct orientations of maximum horizontal stress, probing the conditions necessary for sustained slip within the fracture network or activating the main fault. Additionally, we assess the feasibility of nucleating dynamic rupture earthquake cascades from a distant fracture and investigate the sensitivity of fracture network cascading rupture to the effective normal stress level. We model either pure cascades or main fault rupture with limited off-fault slip. We find that cascading ruptures within the fracture network are dynamically feasible under certain conditions, including: (a) the fracture energy scales with fracture and fault size, (b) favorable relative pre-stress of fractures within the ambient stress field, and (c) close proximity of fractures. We find that cascading rupture within the fracture network discourages rupture on the main fault. Our simulations suggest that fractures with favorable relative pre-stress, embedded within a fault damage zone, may lead to cascading earthquake rupture that shadows main fault slip. We find that such off-fault events may reach moment magnitudes up to M w ≈ 5.5, comparable to magnitudes that can be otherwise hosted by the main fault. Our findings offer insights into physical processes governing cascading earthquake dynamic rupture within multiscale fracture networks.

Fluid-filled cracks sustain a slow guided wave (Krauklis wave or crack wave) whose resonant frequencies are widely used for interpreting long period (LP) and very long period (VLP) seismic signals at active volcanoes. Significant efforts have been made to model this process using analytical developments along an infinite crack or numerical methods on simple crack geometries. In this work, we develop an efficient hybrid numerical method for computing resonant frequencies of complex-shaped fluid-filled cracks and networks of cracks and apply it to explain the ratio of spectral peaks in the VLP signals from the Fani Maoré submarine volcano that formed in Mayotte in 2018. By coupling triangular boundary elements and the finite volume method, we successfully handle complex geometries and achieve computational efficiency by discretizing solely the crack surfaces. The resonant frequencies are directly determined through eigenvalue analysis. After proper verification, we systematically analyze the resonant frequencies of rectangular and elliptical cracks, quantifying the effect of aspect ratio and crack stiffness. We then discuss theoretically the contribution of fluid viscosity and seismic radiation to energy dissipation. Finally, we obtain a crack geometry that successfully explains the characteristic ratio between the first two modes of the VLP seismic signals from the Fani Maoré submarine volcano in Mayotte. Our work not only reveals rich eigenmodes in complex-shaped cracks but also contributes to illuminating the subsurface plumbing system of active volcanoes. The developed model is readily applicable to crack wave resonances in other geological settings, such as glacier hydrology and hydrocarbon reservoirs.

Numerous studies have reported the occurrence of aseismic slip or slow slip events along faults induced by fluid injection. However, the underlying physical mechanism and its impact on induced seismicity remain unclear. In this study, we develop a numerical model that incorporates fluid injection on a fault governed by rate-and-state friction to simulate the coupled processes of pore-pressure diffusion, aseismic slip, and dynamic rupture. We establish a field-scale model to emulate the source characteristics of induced seismicity near the Dallas-Fort Worth Airport (DFWA), Texas, where events with lower-stress drops have been observed. Our numerical calculations reveal that the diffusion of fluid pressure modifies fault criticality and induces aseismic slip with lower stress drop values (<1 MPa), which further influence the timing and source properties of subsequent seismic ruptures. We observe that the level of pore-pressure perturbation exhibits a positive correlation with aseismic-stress drops but a reversed trend with seismic-stress drops. Simulations encompassing diverse injection operations and fault frictional parameters generate a wide spectrum of slip modes, with the scaling relationship of moment (M 0) with ruptured radius (r 0) following an unusual trend, M0∝r04.4 ${M}_{0}\propto {r}_{0}^{4.4}$, similar to M0∝r04.7 ${M}_{0}\propto {r}_{0}^{4.7}$ observed in the DFWA sequence. Based on the consistent scaling, we hypothesize that the lower-stress-drop events in the DFWA may imply less dynamic ruptures in the transition from aseismic to seismic slip, located in the middle of the broad slip spectrum, as illustrated in our simulations.

The slow spreading Reykjanes Ridge overlies the Iceland hotspot and has undergone well ordered changes in crustal segmentation. Previous studies have attributed these changes to varying mantle plume thermal effects, rendering the lithosphere ductile or brittle. Here we use seafloor spreading magnetic anomalies to show that crustal accretion has been focused throughout its spreading history and to determine the detailed evolution of Reykjanes Ridge segments. By ∼53 Ma, organized spreading had developed on an orthogonally spreading linear axis following continental breakup. After a plate motion change at ∼38 Ma, orthogonally spreading offset ridge segments formed by ridge propagation forming varying length fracture zones. From then to the present, the offset segments diachronously migrated back to the original linear geometry from north to south replacing orthogonal with oblique spreading as the axis became linear again. Fracture zones were not terminated, however, simply by reducing segment offsets even to zero. Their termination involved the axial propagation of buoyant upwelling instabilities across the discontinuities, correspondingly extending V-shaped crustal ridges southward. This evolution was guided by a persistent linear deep damp mantle melting interval maintained by the episodic propagation of buoyant upwelling instabilities. Our study indicates that at slow spreading ridges, where buoyant upwelling instabilities govern crustal segmentation, spatial gradients in mantle melting properties may direct the behavior of the instabilities. Where ridges overlie regional hotspot gradients in mantle melting, buoyant instabilities may propagate systematically, and plate boundary evolution may follow an organized pattern.

We develop a 3-D isotropic shear velocity model for the Alaska subduction zone using data from seafloor and land-based seismographs to investigate along-strike variations in structure. By applying ambient noise and teleseismic Helmholtz tomography, we derive Rayleigh wave group and phase velocity dispersion maps, then invert them for shear velocity structure using a Bayesian Monte Carlo algorithm. For land-based stations, we perform a joint inversion of receiver functions and dispersion curves. The forearc crust is relatively thick (35–42 km) and has reduced lower crustal velocities beneath the Kodiak and Semidi segments, which may promote higher seismic coupling. Bristol Bay Basin crust is relatively thin and has a high-velocity lower layer, suggesting a dense mafic lower crust emplaced by the rifting processes. The incoming plate shows low uppermost mantle velocities, indicating serpentinization. This hydration is more pronounced in the Shumagin segment, with greater velocity reduction extending to 18 ± 3 km depth, compared to the Semidi segment, showing smaller reductions extending to 14 ± 3 km depth. Our estimates of percent serpentinization from VS reduction and VP/VS are larger than those determined using VP reduction in prior studies, likely due to water in cracks affecting VS more than VP. Revised estimates of serpentinization show that more water subducts than previous studies, and that twice as much mantle water is subducted in the Shumagin segment compared to the Semidi segment. Together with estimates from other subduction zones, the results indicate a wide variation in subducted mantle water between different subduction segments.

Models of near-field tsunamis and an extreme hurricane provide further evidence for a great precolonial earthquake along the Puerto Rico Trench. The models are benchmarked to brain-coral boulders and cobbles on Anegada, 125 km south of the trench. The models are screened by their success in flooding the mapped sites of these erratics, which were emplaced some six centuries ago. Among 25 tsunami scenarios, 19 have megathrust sources and the rest posit normal faulting on the outer rise. The modeled storm, the most extreme of 15 hurricanes of category 5, produces tsunami-like bores from surf beat. In the tsunami scenarios, simulated flow depth is 1 m or more at all the clast sites, and 2 m or more at nearly all, given either a megathrust rupture 255 km long with 7.5 m of dip slip and M8.45, or an outer-rise rupture 130 km long with 11.4 m of dip slip and M8.17. By contrast, many coral clasts lie beyond the reach of simulated flooding from the extreme hurricane. The tsunami screening may underestimate earthquake size by neglecting trees and shrubs that likely impeded both the simulated flows and the observed clasts; and it may overestimate earthquake size by leaving coastal sand barriers intact. The screening results broadly agree with those from previously published tsunami simulations. In either successful scenario, the average recurrence interval spans thousands of years, and flooding on the nearest Caribbean shores begins within a half-hour.

In this study, we used the Pn tomography method to obtain detailed velocity and anisotropy structures of the uppermost mantle beneath Sakhalin–Kuril–Kamchatka region for improving the understanding of plate subduction, arc–arc collision, and volcanism. We found low Pn velocities beneath volcanoes and areas characterized by pronounced tectonic activity and high Pn velocities with strong anisotropy in the subducting plate. Low Pn velocity anomalies beneath southern Sakhalin connected the low-velocity anomalies in the mantle and crust, indicating the ascent of fluid or melt, and may provide a magmatic source for the rear-arc Rishiri volcano. In the absence of plate subduction, low-velocity anomalies and north–south Pn anisotropy manifested beneath northern Kamchatka, revealing the lateral propagation of mantle flow beneath this northern region. We suggest that the eastern boundary of the slab window at the Kamchatka–Aleutian junction is likely to be located near the Komandorsky Islands.

The southeastern margin of the Tibetan Plateau has undergone complex deformation since the Cenozoic, resulting in a high level of seismicity and seismic hazard. Knowledge about the seismic anisotropy provides important insight about the deformation mechanism and the regional seismotectonics beneath this tectonically active region. In this study, we conduct fullwave multi-scale tomography to investigate the seismic anisotropy in the southeastern margin of the Tibetan Plateau. Broadband records at 111 permanent stations in the region from 470 teleseismic events are used to obtain 5,216 high-quality SKS splitting intensity (SI) measurements, which are then inverted in conjunction with 3D sensitivity kernels to obtain an anisotropic model with multi-scale resolution. Resolution tests show that our data set recovers anisotropy anomalies reasonably well on the scale of 1° × 1° horizontally and ∼100 km vertically. Our result suggests that in the southeastern margin of the Tibetan Plateau the deformation in the lithosphere and asthenosphere are decoupled. The anisotropy in the lithosphere varies both laterally and vertically as a result of dynamic interactions of neighboring blocks as well as lithospheric reactivation. The anisotropy in the asthenosphere largely follows the direction of regional absolute plate motion. The SKS splittings observed at the surface are shown to be consistent with the vertical integral of our depth-dependent anisotropy model over lithospheric and asthenospheric depths.

Ambient noise tomography has become a popular method in the past two decades to image the crust and uppermost mantle structure. To date, broadband Rayleigh wave signals can be obtained from ambient noise, which can be utilized to study the earth's interior structure from the surface down to ∼200–300 km depths. However, it is hard to extract intermediate- and long-period (>50 s) Love wave signals from ambient noise using conventional data processing techniques for ambient noise. Array-based data processing techniques can enhance weak signals. In this study, we adopt a modified algorithm of the double-beamforming method to extract broadband Love wave signals from ambient noise. We validate the accuracy of the dispersion curves measured from our method by comparing them to those measured from the conventional method. Then, we use a finite frequency ambient noise tomography method to construct broadband Love wave phase velocity maps across the contiguous USA. These phase velocity maps are consistent with those obtained from conventional methods at short periods (<40 s). Finally, we analyze the resolution of our double-beamforming method based on checkerboard tests and find that the resolution of phase velocity maps based on our method is close to the aperture of the subarrays used in our double-beamforming method.

Ocean bottom seismometers (OBSs) have been used to detect submarine structural and tectonic information for decades. According to signal source controllability, OBS data have generally been classified into active and passive source data categories. The former mainly focuses on the compressional wave (P-wave) velocity inversion and always lacks valid information about the shear wave (S-wave) velocity structure. While the latter provides structural information with limited resolution due to the aperture of the stations. Overcoming the barriers between processing these two data types will allow the reuse of a vast amount of data from active source experiments to explore the submarine S-wave velocity structural properties. Here, we creatively applied ambient noise interferometry to invert the S-wave velocity structure using data from active source OBS deployment conducted in the southernmost Mariana subduction zone, which had already been utilized to detect submarine P-wave velocity structure. Considering the short time duration and relatively low quality of this type of data, a combined method of short-segment cross-correlation and selected time-frequency domain phase-weighted stacking was adopted to obtain stable cross-correlation functions, which were subsequently used to invert S-wave velocity structures. Compared to previous studies using different methods, our result sheds new light on the crust and upper mantle structure of the southernmost Mariana subduction zone. This method could be used to detect more information based on the reutilization of existing active source OBS data.

Two adjacent groundwater wells on the North China Platform are used to study how earthquakes impacted aquifers. We use the response of water level to solid Earth tides to document changes after earthquakes and how aquifer and fracture properties recovered to pre-earthquake properties. We consider two models for the phase and amplitude of water level response to the lunar diurnal (O1) and semidiurnal (M2) tides: a leaky aquifer model, and a model in which fracture orientation determines the response. In the leaky aquifer model, changes arise from changes in permeability and storage; in the fracture model, changes are due to changes in apparent orientation of transmissive fractures. Responses in one well are best explained by the leaky aquifer model, and can explain the large amplitude coseismic water level and permeability changes and the non-recoverable changes after the largest earthquake. Responses in the other well are consistent with the fracture model and show little coseismic change in water level but changes in apparent fracture orientation. Larger ground motions lead to larger coseismic water level changes and longer recovery times. We propose that the well in the more permeable and shallow aquifer has less variable pore-pressures around the well. Larger coseismic strains from water level changes may enable longer-lasting changes in aquifer properties. We conclude that relatively high permeability aquifers are less susceptible to impacts from seismic waves, and thus have small changes in water levels and hydrogeological properties.

Deformation experiments on hematite characterize its slip-rate dependent frictional properties and deformation mechanisms. These data inform interpretations of slip behavior from exhumed hematite-coated faults and present-day deformation at depth. We used a rotary-shear apparatus to conduct single-velocity and velocity-step experiments on polycrystalline specular hematite rock (∼17 μm average plate thickness) at slip rates of 0.85 μm/s to 320 mm/s, displacements of primarily 1–3 cm and up to 45 cm, and normal stresses of 5 and 8.5 MPa. The average coefficient of friction is 0.70; velocity-step experiments indicate velocity-strengthening to velocity-neutral behavior at rates <1 mm/s. Scanning electron microscopy showed experimentally generated faults develop in a semi-continuous, thin layer of red hematite gouge. Angular gouge particles have an average diameter of ∼0.7 μm, and grain size reduction during slip yields a factor of 10–100 increase in surface area. Hematite is amenable to (U-Th)/He thermochronometry, which can quantify fault-related thermal and mechanical processes. Comparison of hematite (U-Th)/He dates from the undeformed material and experimentally produced gouge indicates He loss occurs during comminution at slow deformation rates without an associated temperature rise required for diffusive loss. Our results imply that, in natural fault rocks, deformation localizes within coarse-grained hematite by stable sliding, and that hematite (U-Th)/He dates acquired from ultracataclasite or highly comminuted gouge reflect minor He loss unrelated to thermal processes. Consequently, the magnitude of temperature rise and associated thermal resetting in hematite-bearing fault rocks based on (U-Th)/He thermochronometry may be overestimated if only diffusive loss of He is considered.

We discuss the tectonic structure, seismic stratigraphy and evolution of the NW Sulu Sea using reprocessed 2D reflection profiles. The NW Sulu Sea is located between the Palawan continental shelf and the Cagayan Ridge and represents the northern part of the Sulu Sea, a marginal sea resulting from Paleogene extension and subsequent Neogene contraction due to convergence between the Palawan and the Philippine blocks. The basin consists of six seismo-stratigraphic units overlying crystalline basement. Syn-orogenic depocenters contain calibrated Middle Miocene to, possibly, Lower Miocene units, while rift-related depocenters consist of uncalibrated but tentatively dated Paleogene to Lower Miocene units. Thickness and depth maps of the main units and bounding horizons differentiate the Piedra-Blanca and the Rasa domains, separated by the NW-Sulu-Break major tectonic structure. Fault-bounded rift-related depocenters are strongly segmented. We interpret that NW-SE and NE-SW trending zones accommodate shape and trend variations of these depocenters. We suggest that these zones may link rift segments, recording different extensional deformation. Miocene thrusting and folding in the Piedra-Blanca Domain and mudflow with associated gravitational structures in the Rasa Domain influenced the deposition of syn-orogenic units. Rift structures inherited from rift segmentation may have conditioned the style and distribution of contractional deformation during the subsequent incipient reactivation during contraction. In the context of SE Asia, our results support that the timing of rifting of the NW Sulu Sea overlaps with the opening of the South China Sea and the North Palawan margin, which may indicate a common geodynamic driving force triggering extension.

Synthesized polycrystalline samples composed of enstatite and olivine with different volumetric ratios were deformed in compression under anhydrous conditions in a Paterson gas-medium apparatus at 1150–1300°C, an oxygen fugacity buffered at Ni/NiO, and confining pressures of 300 or 450 MPa (protoenstatite or orthoenstatite fields). Mechanical data suggest a transition from diffusion to dislocation creep with increasing differential stress for all compositions. Microstructural analyses by optical and scanning electron microscopy reveal well-mixed aggregates and homogeneous deformation. Crystallographic preferred orientations measured by electron backscatter diffraction are consistent with activation of the slip systems (010)[100] and (010)[001] for olivine and (100)[001] and (010)[001] for enstatite, as expected at these conditions. Nonlinear least-squares fitting to the full data set from each experiment allowed the determination of dislocation creep flow laws for the different mixtures. The stress exponent is 3.5 for all compositions, and the apparent activation energies increase slightly as a function of enstatite volume fraction. Within the limits of experimental uncertainties, all two-phase aggregates have strengths that lie between the uniform strain rate (Taylor) and the uniform stress (Sachs) bounds calculated using the dislocation creep flow laws for olivine and enstatite. Calculation of the Taylor and Sachs bounds at strain rate and temperature conditions expected in nature (but not extrapolating in pressure) indicates that using the dislocation creep flow law for monomineralic olivine aggregates provides a good estimate of the viscosity of olivine-orthopyroxene rocks deforming by dislocation creep in the deeper lithosphere and asthenosphere.

The start of the Paleoproterozoic supercontinent cycle is typically taken as the initiation of orogenesis at ca. 2.1 Ga leading to the assembly of Earth's first supercontinent, Columbia. However, the dearth of ca. 2.5–2.2 Ga geological records makes it difficult to deduce tectonic factors during the onset of the Paleoproterozoic supercontinent cycle. The petrogenesis of tonalite–trondhjemite–granodiorite (TTG) provides useful proxies for tracing prevailing geodynamic regimes of early continental evolution. However, marked decreases of TTG and other magmatism occurred across the Archean–Paleoproterozoic transition and have previously precluded forming testable hypotheses. Early Paleoproterozoic TTGs have been identified in the North China Craton (NCC) and other cratons, which may represent the last major pulse of TTGs globally. Here we present low δ18O and δ30Si ca. 2.3 Ga TTGs from the NCC, together with thermodynamic modeling and compilation of stable O and Si isotopes for TTGs globally through time. The ca. 2.3 Ga TTGs were derived from the partial melting of Archean basaltic crust and give lighter average zircon δ18O (3.15 ± 0.35‰) and whole-rock δ30Si values (−0.17 ± 0.08‰) than most Archean TTGs. Considering coeval mafic-felsic igneous rocks, and lithospheric thinning since ca. 2.5 Ga based on estimated crustal thickness through the Neoarchean–Paleoproterozoic, we posit the onset of an intraplate rifting consistent with the anomalous low-δ18O magmatism. Continental rifting of Archean cratons/supercratons plausibly hints at the formation of rifts driving subduction initiation as the veritable onset of the Paleoproterozoic supercontinent cycle.

Geodetic methods can monitor changes in terrestrial water storage (TWS) across large regions in near real-time. Here, we investigate the effect of assumed Earth structure on TWS estimates derived from Global Navigation Satellite System (GNSS) displacement time series. Through a series of synthetic tests, we systematically explore how the spatial wavelength of water load affects the error of TWS estimates. Large loads (e.g., >1,000 km) are well recovered regardless of the assumed Earth model. For small loads (e.g., <10 km), however, errors can exceed 75% when an incorrect model for the Earth is chosen. As a case study, we consider the sensitivity of seasonal TWS estimates within mountainous watersheds of the western U.S., finding estimates that differ by over 13% for a collection of common global and regional structural models. Errors in the recovered water load generally scale with the total weight of the load; thus, long-term changes in storage can produce significant uplift (subsidence), enhancing errors. We demonstrate that regions experiencing systematic and large-scale variations in water storage, such as the Greenland ice sheet, exhibit significant differences in predicted displacement (over 20 mm) depending on the choice of Earth model. Since the discrepancies exceed GNSS observational precision, an appropriate Earth model must be adopted when inverting GNSS observations for mass changes in these regions. Furthermore, regions with large-scale mass changes that can be quantified using independent data (e.g., altimetry, gravity) present opportunities to use geodetic observations to refine structural properties of seismologically derived models for the Earth's interior structure.