Geophysical Journal International

SummaryThe probability that any given time interval of duration τ is occupied by one earthquake, or more, characterizes the temporal statistics of seismic sequences and, therefore, the temporal clustering of events. The occupation probability reveals the fractal behaviour of seismic sequences, $\Phi (\tau ) \sim \tau ^{1 - D_\tau }$, at short times, defining a temporal fractal dimension, Dτ, for seismic events. We introduce an empirical model of the occupation probability, parametrized by the fractal dimension and two other parameters. We use the mathematical relationship between the occupation probability and the inter-event time probability density to develop intuition about the model parameters and to compare the proposed model with the conventional gamma model for inter-event times. Our model captures the statistical properties of a wide range of seismic sequences where the gamma model fails, from slow-slip driven swarms to burst-like episodes at the bottom of the seismogenic zone and low-frequency earthquakes. Using real and synthetic catalogues, we find that the model parameters are related to the degree of intermittency of the seismic activity, the characteristics of the bursts and the time scale of the quiescent periods. Measuring and modelling the occupation probability constitutes a valuable tool to categorize seismic sequences over a wide spectrum of seismic occurrence patterns. When applied to the new generation of earthquake catalogues, this empirical method highlights intermittency and fractal bursts, challenging conventional seismicity models and emphasizing the need to refine them to better capture these key features.

SummaryNowadays, many joint inversions are carried out to understand the earthquake source process. In the joint inversion analysis, we have to determine the relative weights among different datasets in addition to the regularization term, such as smoothing. Akaike’s Bayesian Information Criterion (ABIC) is known to be useful to find the appropriate values of such hyperparameters. This study proposes a method to jointly invert tsunami, GNSS, and SAR data using ABIC to construct a finite fault model. We demonstrate our inversion scheme in the case of the 2024 Noto Peninsula earthquake, whose fault geometry is still under discussion. Since the dip angle of the fault can also be considered as a hyperparameter, we evaluate three types of dip angles and estimate an appropriate value based on ABIC. In other words, our inversion scheme utilizes ABIC to determine the dip angle, the weights among datasets, and the spatial smoothness of fault slip. Our fault model indicates that (1) listric fault, varying the dip angle with depth, is the most appropriate among the ones we proposed, (2) the largest slip is on the fault under the northwestern corner of the peninsula, and (3) coseismic fault slip extends to offshore faults east of the peninsula. In the case of the listric fault, ABIC values GNSS and SAR data, which improves the agreement of the on-land coseismic displacement while also reproducing tsunami data. We also find that analyzing tsunami records in the frequency domain helps to obtain a robust inversion result when employing ABIC.

SummaryThe elastic deformation of Earth’s surface caused by internal mass distribution varies significantly across loading models, especially in high-precision applications. Although several studies have applied loading corrections to Global Navigation Satellite System (GNSS) time series in mainland China, discrepancies between models, particularly those involving Gravity Recovery and Climate Experiment (GRACE) data and its downscaled derivatives, remain insufficiently explored. Moreover, previous research has not comprehensively assessed vertical crustal deformation after applying different environmental loading corrections. This study systematically evaluates the correction effects of various environmental loading models on GNSS vertical displacements across mainland China, generating vertical velocity maps along with their associated uncertainties. The results show that hydrological loading (HYDL) has the most significant impact on GNSS vertical displacements, whereas non-tidal oceanic loading (NTOL) has the least effect. Substantial differences exist between various HYDL models, while discrepancies between non-tidal atmospheric loading (NTAL) and NTOL models are relatively minor. A comparison of correction effects between the HYDL model and GRACE data reveals that the HYDL model offers more accurate corrections, whereas downscaled GRACE data demonstrates improved performance, underscoring its potential advantages. After applying loading corrections and filtering common mode error (CME), the uncertainty in GNSS vertical velocity is notably reduced, although velocity variation remains small. This effect is also evident in seasonal variations. Furthermore, a comparison of vertical land motion (VLM) constrained by different HYDL models and downscaled GRACE data with VLM constrained by an independent land-water fusion model reveals higher consistency between the downscaled GRACE data and the independent model in the North China and northwestern Tianshan regions, suggesting that VLM derived from downscaled GRACE data may be more reliable. We also quantify the combined impact of geocenter motion and glacial isostatic adjustment (GIA) on the VLM trend across mainland China, estimated it at approximately 0.13 mm/yr. While the spatial characteristics of the VLM trend show minimal changes after correction, its intensity is significantly affected. This study provides crucial insights into the correcting of environmental loading effects in GNSS vertical displacements and contributes the latest observational results on vertical crustal deformation in mainland China.

SummaryThis study presents an enhanced method for integrating Global Navigation Satellite System (GNSS) and Interferometric Synthetic Aperture Radar (InSAR), referred to as Iterative Decomposition-based GNSS-enhanced InSAR (ID-GInSAR), to address both spatially correlated components (SCCs) and topographically correlated components (TCCs) in interferogram errors. While traditional GInSAR (GNSS-enhanced InSAR) is effective in mitigating long-wavelength SCCs, it often overlooks TCCs, which are particularly significant in regions with steep topographic gradients. The proposed ID-GInSAR approach employs an iterative decomposition process to decouple and independently model SCCs and TCCs, utilizing a combination of exponential and statistical models. The method is validated using Sentinel-1 SAR and GNSS data from California’s southern Central Valley. Results demonstrate that ID-GInSAR significantly lowers noise in interferograms, enhances the robustness of displacement time series, and improves the accuracy of co-seismic deformation and velocity field estimates. Specifically, ID-GInSAR reduces the root mean square (RMS) between GNSS and InSAR by up to 55% in individual interferograms and by an average of 30.4% in displacement time series compared to traditional GInSAR methods. Furthermore, ID-GInSAR effectively highlights subtle transient deformation, such as coseismic offsets, and provides more robust velocity fields over shorter time spans (less than three years). Finally, we compare our method with other approaches, including Remove/Filter/Restore (RFR) and GACOS, and discuss their applicability scenarios. Collectively, ID-GInSAR provides an alternative integration method for regions with complex topography where ground-based GNSS observations are available.

SummaryRupture directivity significantly increases horizontal peak ground acceleration, elongates aftershock clouds, and enlarges meizoseismal areas beyond the fault end in front of the direction of rupture propagation. In this study, we examine the directivity of 25 moderate to large earthquakes (Mw ≥ 6) from 1968 to 2017 in the Iranian plateau by employing relocated earthquake clusters, mapped surface ruptures, focal mechanisms of earthquakes, slip distribution models, spatial distribution of Peak Ground Acceleration (PGA) amplitudes and macroseismic effects. The methodology overcomes the lack of dense seismic networks required to study directivity using methods based on the azimuthal variation of the spectrum of seismic waves. We show that 16 out of the 25 (i.e., 64%) of the earthquakes investigated have mostly unidirectional rupture. This implies that unidirectional ruptures in a slow deforming continental collision zone such as the Iranian Plateau is only slightly less common than those observed globally. With the understanding that unidirectional rupture increases the probability of ground shaking off the termination of the causative faults, our findings highlight the importance of considering the directivity effect in earthquake hazard assessment in Iran and also in other slow deforming continental regions.

SummaryEvidence from seismic studies, mineral physics, thermal evolution models and geomagnetic observations is inconclusive about the presence of a stably stratified layer at the top of the Earth’s fluid outer core. Such a convectively stable layer could have a strong influence on the internal fluid waves propagating underneath the core-mantle boundary (CMB) that are used to probe the outermost region of the core through the wave interaction with the geomagnetic field and the rotation of the mantle. Here, we numerically investigate the effect of a top stable layer on the outer core fluid waves by calculating the eigenmodes in a neutrally stratified sphere permeated by a magnetic field with and without a top stable layer. We use a numerical model, assuming a flow with an m-fold azimuthal symmetry, that allows for radial motions across the lower boundary of the stable layer and angular momentum exchanges across the CMB through viscous and electromagnetic coupling. On interannual time-scales, we find torsional Alfvén waves that are only marginally affected by weak to moderate stratification strength in the outer layer. At decadal time-scales similarly weak stable layers promote the appearance of waves that propagate primarily within the stable layer itself and resemble Magneto-Archimedes-Coriolis (MAC) waves, even though they interact with the adiabatic fluid core below. These waves can exert viscous and electromagnetic torques on the mantle that are several orders of magnitude larger than those in the neutrally stratified case.

SummaryA novel time-lapse modelling scheme for Airborne Electromagnetics (AEM) monitoring datasets is presented, using data from multiple surveys applied to study the hydro-related evolution of the Bookpurnong floodplain in South Australia. Additionally, it introduces a new wide-ranging approach for this type of study, incorporating new processing, validation, and interpretation tools.Time-Lapse studies are widespread in the literature but are not commonly applied to model EM data, particularly AEM data. This is linked to the challenges of performing overlapping data acquisition with inductive systems. The key features of the new time-lapse scheme, which address these issues, include the definition of independent forward and model meshes, essential for considering discrepancies in the location of soundings which arise in multitemporal AEM data acquisition, and the incorporation of system flight height in the inversion. This proved crucial for achieving satisfactory data fitting and limiting artifact propagation in the time-lapse models.Additionally, a novel processing workflow for AEM multitemporal datasets is presented. This has proven important for effectively processing the multitemporal datasets, which presents new challenges in identifying noise coupling arising from the use of different systems across vintages of data, possible variations in acquisition settings operated by different field crews, and changes in subsurface resistivity in the survey area. Results generated from the time-lapse modelling are evaluated with an Independent Hydrogeological Validation (IHV), designed to support the geophysical models validation and interpretation by providing a first-step hydrogeological evaluation.At Bookpurnong, along a sector of the Murray River floodplain, multitemporal AEM surveys were collected in 2015, 2022 and 2024, to study natural and engineered changes in the groundwater system over time. The time-lapse models show significantly smaller variations compared to those determined with individually modelled survey data sets, while delineating sharply bounded changes in resistivity across the floodplain. This demonstrates the effectiveness of the new time-lapse scheme in minimizing inversion variations typically encountered with independently modelled results affected by larger equivalence issues.Here, AEM models are first compared with resistivity borehole measurements, revealing a strong match between the two methodologies and spatial variations in resistivity consistent with a meandering river across the floodplain. These variations are further validated and interpreted using the IHV approach, which revealed a direct correlation between the hydrological stress of the Murray River and the response of shallow aquifers. Additionally, time-lapse geophysical models, combined with a hydrostratigraphic analysis, allow for a direct correlation between shallow and deep hydrogeological responses.We believe that the time-lapse methodology described here can be widely applied to multitemporal studies using AEM datasets, enabling the study of a broad range of natural processes with great accuracy and at the basin scale.

SummaryThe complex three-dimensional (3D) geometry of active faults plays a crucial role in controlling earthquake location, extent, and rupture behavior, making the accurate representation of fault models essential. Fault structures are typically interpreted manually from relocated hypocenters and interpolated to generate 3D fault surfaces. However, this process is often non-unique and uncertain due to the uneven spatial distribution of earthquake hypocenters, the subjectivity of manual interpretation, and the complexity of non-planar faults. To address these challenges, we developed a method that combines adaptive threshold hierarchical clustering with quantitative evaluation to automatically and effectively construct 3D models of seismogenic faults. This method utilizes the nearest neighbor index (NNI) to determine whether seismic activity exhibits clustering characteristics indicative of fault structures. Adaptive threshold hierarchical clustering is subsequently applied to identify small earthquake clusters associated with each fault. High-density 3D automatic slicing ensures robust fitting of fault lines, and in combination with surface rupture data, discrete smooth interpolation (DSI) is used to construct a 3D fault model. For each fault, we calculate distances from small earthquake clusters to the 3D fault structure and analyze their spatial distribution using kernel density estimation (KDE) to optimize the model for a near-symmetric distribution of small earthquake clusters on both sides of the fault. We applied this method to the 2013 Ms 7.0 and 2022 Ms 6.1 earthquakes in southern Longmenshan, Sichuan, China, refining the 3D seismogenic fault models for both events. Additionally, we constructed a 3D fault model for the 2019 Ridgecrest Mw 7.1 earthquake sequence using the same approach. The results indicate that this method is applicable to both individual faults and multiple intersecting fault systems. Compared to traditional manual modeling approaches, our method significantly enhances the identification of small earthquake clusters, reduces reliance on manual interpretation, increases modeling efficiency, and minimizes errors. This innovative modeling technique advances the 3D geometric construction of complex active faults and is adaptable to a wide range of seismic research applications.

SummaryAn 18-layer, U-shaped convolutional neural network (CNN) was trained to predict Vs models and identify near-surface void locations. To enhance seismic inversion accuracy for real world applications, the model is trained on synthetic datasets augmented with field noise. While models trained on noise-augmented data showed poorer performance on synthetic testing datasets, they achieved lower root mean square error values and the best results on field data. The velocity model resulting from full-waveform inversion based on noise-augmented model accurately detected a void-like low velocity zone near the known void location. This approach shows that training with field-noise-augmented data allows machine learning models to generalize better to real-world conditions, increasing their reliability for velocity inversion in noisy environments. The results highlight the strong potential of this strategy, particularly if a diverse range of real noise samples is incorporated during training.

SummaryCoseismic slip models for large (Mw > 7) strike-slip earthquakes present a variety of shallow slip deficit (SSD). Accurate estimate of SSD is difficult, and it has been suggested that SSD are to some degree associated with fault-zone characteristics, incompleteness of data coverage as well as simplified model assumptions. Furthermore, SSD can also be sensitive to the amount of model smoothness adopted. Since phase gradient are sensitive to the missing shallow slip from our simulated data, we performed a synthetic test and presented a case study of the 2019 Ridgecrest earthquake sequence to validate that phase gradient from radar interferometry could help reveal the actual SSD for kinematic slip models even without enough near-fault observation. Our results indicate that even in the presence of a greater degree of observational gaps, the phase gradient can still nearly substitute for near-fault observations in constraining the shallow slip. Lastly, we provide a preferred coseismic slip model constrained by all available observations including phase gradient, but with 4-km data gap near the fault trace. This model results in ∼35% SSD for the Ridgecrest earthquakes, matching previous estimates that incorporate near-field data. Considering the phase gradient approach is a straightforward mathematical operation, this approach may be applicable to other types of earthquakes. Notably, due to the smaller amplitude and lower signal-to-noise ratio for the phase gradient data, one needs to carefully balance the trade-offs among weights of different datasets and model smoothness.

SummaryElectrokinetic signals, such as surface self-potential (SP) variations, offer a unique window into coupled fluid–electrical processes in the Earth’s crust, yet their quantitative interpretation remains challenging in complex geological settings. In this study, we develop an electrokinetic modeling framework by extending the modified Luco-Apsel-Chen generalized reflection and transmission method to simulate SP responses due to a fluid-injection source in layered geological media. After simulating electric signals, we apply location-specific amplification factors—derived from prior numerical investigations—that account for the effects of steel well casings. This post-processing step enables rigorous comparison with field observations. Using the well-documented deep fluid injection experiment at the Soultz-sous-Forêts geothermal site, we calibrate simulated pore pressure against downhole measurements to derive a realistic source function for direct comparison between modeled and observed SP signals. The model reproduces key spatiotemporal features of the mV-scale SP anomalies and, importantly, captures the observed slower decay of surface SP signals after shut-in despite the rapid decrease in deep pore pressure. Previous field-scale studies have qualitatively attributed this phenomenon to sustained ionic transport; here, our simulation results provide a quantitative demonstration that this process—driven by continued pore-fluid movement—is responsible for the slower SP decay, a mechanism not captured in earlier electrokinetic simulations. These findings provide new mechanistic insight into SP generation in stratified media, demonstrate the essential role of casing effects in field-scale interpretation, and establish a transferable framework for monitoring subsurface fluid flow in geothermal, hydrocarbon, and groundwater systems.

SummaryTsunami warning systems implemented worldwide rely on the fast characterization of earthquake sources, in particular on the estimation of the moment magnitude Mw. Reliable estimation of Mw typically takes 10 to 20 minutes for large events based on conventional seismic signals. A promising alternative is the use of Prompt Elasto-Gravity Signals (PEGS), which are very low-amplitude gravitational perturbations induced by earthquakes that travel at the speed of light and can be recorded by broadband seismometers at time scales of a few minutes after origin time. We propose here a first implementation of real-time PEGS analysis to enhance the Peruvian earthquake monitoring system by enabling rapid magnitude estimation for large and potentially tsunamigenic earthquakes. We train a graph neural network to recognize the space-time structure of PEGS over a large international set of broadband seismic stations, even when their amplitudes are smaller than the noise level, and to rapidly estimate the magnitude and location of the source. Our results indicate that the PEGS-based system can estimate the magnitude of Mw ≥ 8.2 earthquakes, within 5 minutes after the event’s initiation, with sufficient accuracy for tsunami warning purposes. Simulated real-time tests confirm the viability of the PEGS-based approach for operational early warning, providing robust source estimations of large magnitude events to the Peruvian earthquake monitoring system that are valuable for tsunami warning.

SummaryMagnetotelluric (MT) data inversion seeks to recover resistivity models of the subsurface. Solving the inversion problem is a non-trivial task, as multiple plausible solutions can be recovered due to the non-linearity of the problem. To reduce this non-linearity, we propose a data-driven approach where a 1D cumulative resistance model is estimated from MT data via a direct data transformation. We define the cumulative representation of layered models as the weighted sum of layer thickness divided by resistivity from surface to any depth level, which is the cumulative conductance. Its inverse, cumulative resistance, is directly related to the real part of the impedance computed from MT data. We train a neural network to transform the MT impedance into a resistance model. The corresponding 1D resistivity model is obtained without a priori information. We validate our approach using synthetic and real data, opening the discussion for future developments of this new perspective.

SummaryIn this work, we address the characterization of porosity and water saturation in a synthetic model of a shallow alluvial subsurface using frequency electromagnetic and seismic data. The inversion method employs a Gauss-Newton scheme, where the Jacobian of the merged seismic and electromagnetic data is formulated with respect to the spatially heterogeneous petrophysical parameters. This is made possible by introducing realistic petrophysical relationships, which significantly reduce the number of unknowns in the inverse problem and incorporate a strong prior correlation between the information contained in both data types regarding the subsurface composition. The results obtained show that this Simultaneous Joint Petrophysical Inversion (SiJPI) produces reconstructions with clear improvements compared to Independent Petrophysical Inversion (IPI). Indeed, it greatly enhances the spatial resolution of subsurface mapping, as well as the quantitative estimation of porosity and saturation.

SummaryThe Xiaojiang fault system (XJFS) is located in the southeastern margin of the Tibetan Plateau, which has been considered as an ideal site to explore the traces and effects of past tectonic activity. In this study, we obtain a high-resolution P-wave velocity and azimuthal anisotropy model of the upper crust beneath the XJFS, utilizing the 2D Pg wave tomography method including both station and event depth corrections. The upper crust displays obvious heterogeneity of both azimuthal anisotropy and P-wave velocity underlying the XJFS. The large azimuthal anisotropy beneath the XJFS, especially the regions where several faults interact, suggests the cracks are widely distributed. Generally, the upper crust is featured by several high-velocity bodies separated by low-velocity anomalies. The high-velocity bodies are speculated to be related to the remnant magmatic rocks of the Emeishan large igneous province. The low-velocity anomalies are interpreted to represent fault damage zones which could be attributed to the strike-slip faulting/shearing along the faults and upwelling of deeply-sourced partial melts and fluids. The present tectonic activity in the XJFS is characterized by rigid block extrusion along strike-slip faults in the upper crust, which is consistent with the rigid block extrusion model. We further propose a tectonic model to display the evolution of the upper crust beneath the XJFS, in which the faults and fluids activity plays an essential role.

AbstractGeophysical measurements such as induced polarization (IP) are invaluable for understanding the physical properties of rocks, including pore structure, hydraulic properties, and mineral content. However, collecting reliable IP measurements from low-permeability rocks poses substantial challenges due to the difficulty of saturating their tight pore spaces. Additionally, IP measurements on rocks that are not cored to fit conventional sample holders, or are irregularly shaped, are particularly difficult to obtain. In this work, we address these challenges through (1) the use of reliable saturation procedures developed for low-permeability samples, and (2) a molding procedure designed to overcome the difficulties of measuring IP on irregularly shaped or broken rock cores. Core-scale gravimetric porosity measurements closely match values obtained from destructive mercury intrusion porosimetry (MICP) on rock fragments, confirming the effectiveness of the saturation procedure. Direct comparisons of IP measurements between molded and unmolded cores demonstrate that the molding process does not significantly alter the intrinsic electrical response of the samples. Fully saturated mudstones exhibit strong statistically significant relationships between the mean relaxation time (τmean) and permeability (k), and between effective porosity (1/formation factor, F) and interconnected porosity (ϕ) (Archie’s law). Conversely, partial saturation due to ineffective saturation methods introduces substantial scatter to these petrophysical correlations. Overall, these findings underscore the potential of these methods to enhance the reliability and accuracy of SIP measurements on challenging rock samples.

SummaryOn January 18, 2021, a blind mid-crustal Mw ∼6.4 earthquake occurred near San Juan, Argentina. The observation of associated ground deformation with single interferograms is obscured by strong tropospheric signals. We apply appropriate corrections to the data and reconstruct the deformation field associated to the event through InSAR time series approach. We show it is possible to retrieve this signal to invert the fault parameters. The observed ground deformation is consistent with a high angle NW-dipping fault plane at a centroid depth of ∼19 km. The geometry of this fault supports the reactivation of pre-existing structures within the Cuyania Terrane, suggesting a direct structural connection and strain transfer to the actively deforming, east-vergent Precordillera front. We analyze our findings to deduce a static friction coefficient ≤0.3 for mid-crustal faults of the region.

SummaryWe refined the Attenuated ProPagation of Local Earthquake Shaking (APPLES) ground-motion-based earthquake early warning (EEW) approach, and directly compare APPLES performance with that of the source-characterization-based U.S. ShakeAlert EEW system for a suite of historical earthquakes in the U.S. West Coast and Japan. APPLES is an extension of the Propagation of Local Undamped Motion (PLUM) algorithm in which observed shaking intensity at seismic stations is used to forward-predict intensity distributions to surrounding areas using an attenuation model derived from an intensity prediction equation. We test new configuration options within APPLES, such as using the second highest estimated ground motion rather than the maximum, to better match median ground-motion observations and reduce alerts for small magnitude earthquakes, both of which are key alerting priorities within ShakeAlert. We evaluate these configurations alongside ShakeAlert by comparing the ground-motion estimation accuracy and available warning times relative to station observations and ShakeMap distributions. Our preferred APPLES configuration produces accurate ground-motion estimates and corresponds better with median observations compared to ShakeAlert’s estimates. This preferred configuration substantially reduces alert issuance for M < 5.0 earthquakes compared to the previous APPLES configuration, and alert-release criteria can further restrict alerts to primarily M ≥ 5.5 earthquakes without requiring magnitude estimation. Prioritizing matching median-observed ground motions may reduce APPLES warning times compared to configurations that were tuned to avoid missed alerts (such as those that use the maximum estimated ground motions), which can lead to shorter warning times compared to ShakeAlert for the same alert threshold. However, station-based warning time assessments demonstrate that APPLES can outperform ShakeAlert for high target thresholds. APPLES is a simple, independent EEW approach that may improve the robustness of EEW for the West Coast of the U.S.

SummaryMetallic infrastructure, such as steel sheet situated within landfills, poses significant challenges to accurate tracking of leachate using induced polarization (IP) methods. The application of IP method is efficient to delineate leakage; however, the presence of metallic structures can cause an interference on the survey and generate high-chargeability anomalies as observed in field survey. To comprehensively validate the interference caused by steel sheets, both numerical and empirical field tests were conducted. As expected, both results demonstrate that interference diminishes as the distance between survey line and metallic structure increases. Additionally, at consistent intervals, the chargeability values inverted using integral chargeability (IC) exhibit a monotonic increase with depth. Moreover, the interference induced by metallic structures is also affected by the controlling factors (i.e. depth, width and thickness) of the structure alongside the intrinsic resistivity and chargeability. Strategic utilization of the size, chargeability, and spatial positioning of metallic structures relative to survey lines can significantly enhance background polarization. This approach offers a promising framework for improving the spatial resolution of subsurface targets exhibiting low polarization effects. The optimization of survey line placement, which must consider the dimensions and electrical properties of metallic structures such as steel sheets, is essential for accurately characterizing landfill leachate using the IP method.

SummaryNorthern Europe experiences vertical land motion and sea level changes that deviate from the average as a consequence of past changes in ice sheet cover in Fennoscandia and the British Isles. The process, called Glacial Isostatic Adjustment (GIA), is controlled by the subsurface structure. Numerical models of GIA can be compared to observations of uplift or past sea level changes to constrain the subsurface structure, and such models can also be used to correct present-day sea level observations to reveal sea level changes due to climate change. GIA models for northern Europe usually adopt a homogeneous upper mantle viscosity even though seismic studies indicate contrasting elastic lithosphere thickness and upper mantle structure between Northwestern Europe and Eastern Europe. This raises the question whether the effect of lateral variations in structure (3D viscosity) can be detected in observations of GIA and whether including such variations can improve GIA model predictions. In this study we compare model output from a finite element GIA model with 3D viscosity to observations of paleo sea level and current vertical land motion. We use two different methods to derive 3D viscosities, based on seismic velocity anomalies and upper mantle temperature estimates. We use three different reconstructions of the Eurasian ice sheet, one based on an inversion using a 1D model, and two others based on glacial geology and modelling. When we use these two reconstructions, we find that the data are fit better using 3D viscosity models. Models with two separate 1D viscosities for Fennoscandia and for the British Isles cannot replicate a 3D model because a 3D model redistributes GIA-induced stresses differently from a combination of models with separate 2D viscosities. The fit to data across Fennoscandia is improved when, as indicated by seismic models, the upper mantle viscosity is higher than for the rest of Northern Europe. The best fit is obtained with a model with dry olivine rheology, in agreement with other evidence from Fennoscandia.