Geophysical Journal International

SummaryAccurate characterization of the magnitude-frequency distribution of seismicity, and its associated uncertainties, is essential for seismic hazard assessment. This distribution is commonly described by the Gutenberg–Richter (GR) relation, parameterized by the b-value, which has been identified as a potential proxy for investigating many spatiotemporally varying Earth phenomena. Estimating the spatiotemporal variability of b-values often requires windowing, forcing a trade-off between resolution and statistical reliability. New probabilistic methods circumvent this by inferring both the number and locations of change points directly from earthquake catalogs. Nevertheless, accurately determining the b-value remains difficult because the GR relation only holds over a limited range of magnitudes. This research develops a general statistical model to address several methodological challenges in estimating the magnitude-frequency distribution of observed seismicity, including variations in space or time. The approach simultaneously solves for the b-value and magnitude-range limits. This avoids potential bias due to inaccurate manual truncation of earthquake catalogs. The model considers the entire observed catalog and parameterizes the decay of the distribution at both low and high magnitudes. Consequently, robust uncertainties in estimated b-values reflect uncertainty in the range of magnitudes over which the GR relation is observed to be valid. Importantly, spatiotemporal variations in the parameters that define the magnitude range are considered to be independent from the b-value, as we assume the physical factors that influence the GR relation are independent of the factors that limit the observed earthquake catalog. We demonstrate this methodology through application to simulated and observed earthquake catalogs. In particular, the value of our approach is highlighted through application to observed records of induced seismicity associated with fluid-injection operations in western Canada. Our results demonstrate accurate b-value estimates and associated uncertainties. Furthermore, the additional parameters that define the magnitude range serve as proxies for other factors including seismic network performance, recording duration, potential geometric limitations on earthquake size, and potential injection characteristics (in induced seismicity cases). Our approach also allows for the investigation of how these other factors may vary in space/time. Results from this work contribute to rigorous propagation of accurate b-value estimates, including uncertainties, into subsequent analyses such as seismic hazard models and regulatory protocols that are applied to industrial activity.

SummarySeismic emissions from wind turbines (WTs) depend on the rotation of the WT blades and the wind direction-dependent movement of the WT. Mechanical coupling between the WT foundation and the subsurface generates complex seismic wavefields, making it challenging to manually separate the contributions of different signal sources, thus complicating data labelling. We address this challenge by applying unsupervised machine learning techniques that do not require labelled data. Our analysis focuses on seismic WT emissions recorded near Wind Farm Tegelberg in the eastern Swabian Alb, Southwest Germany. Specifically, we extract time-averaged wavelet features by temporal averaging the wavelet transformation of the continuous three-component seismic data and subsequently apply the clustering algorithm Hierarchical Density-Based Spatial Clustering of Applications with Noise (HDBSCAN). The resulting clusters not only capture the variations in the WT rotation rate but also reveal a clear dependency on wind direction, associated with the radiation pattern of different surface waves. Our results demonstrate the potential of HDBSCAN to uncover meaningful, source-related patterns in continuous seismic records.

SummaryThe Semail Ophiolite in Oman represents one of the well-preserved ophiolite complexes globally and provides a unique window into the processes of obduction. The emplacement of Semail ophiolite onto the Arabian lithosphere is a result of intra-oceanic subduction, was strongly influenced by inheritance features preserved from pre-obduction tectonic processes. Therefore, a detailed characterization of the crustal architecture and rheological properties of the lithosphere are essential for improving our understanding of obduction processes. In this study, we investigated the crustal structure and Moho depth beneath the central Oman Mountains through analysis of P-wave receiver functions (PRFs) and Bouguer gravity anomalies. We utilize broadband seismic data recorded at 12 seismic stations spanning the ophiolite belt and surrounding regions (Ghaba basin and Saih Hatat Dome). PRFs analysis reveals noticeable lateral variations in Moho depths ranging from ∼39 km beneath the sedimentary basin to ∼46 km beneath the ophiolite belt, and decreasing to ∼30 km underneath Saih Hatat Dome (SHD). 2D forward modeling of Bouguer gravity anomalies (−36 to 91 mGal) constraints with seismological results shows flexural bending of the Moho topography and thin crust (∼ 30 km) beneath the SHD. The 2D forward flexural modelling analysis suggests that lithospheric flexure is due to the emplacement of the ∼5 km thick Semail Ophiolite. The presence of a thin crust beneath the SHD is caused by Permian rifting and thinning of the continental lithosphere. The observed high value of Vp/Vs (1.75 – 1.87) also provides support for Permian mafic intrusions to the lower crust. The Arabian lithosphere exhibits lower mechanical strength in the southern region (Te = 25 km) relative to the northern area, a characteristic likely inherited from pre-obduction magmatic processes. These results provide new geophysical constraints on the crustal architecture of the southern Oman Mountains and emphasize the role of surface loading in shaping lithospheric structure during ophiolite emplacement.

SummaryThis study aims to evaluate the effectiveness of the remove-restore method applied to GRACE (Gravity Recovery and Climate Experiment) gravity solutions, in which climate-related signals are first removed to allow a more meaningful interpretation of residual gravity signals associated with dynamic processes in Earth’s deep interior. By removing seasonal cycles and long-term trends, the analysis focuses on non-seasonal variations where causal attribution is clearer. Results indicate that climate correction reduces GRACE signal variability by approximately 30% over both oceanic and continental regions, with the strongest impact observed in major river basins. The correction is most effective for temporal scales below 10 years and spatial scales up to spherical harmonic degree 25. While overall variability decreases, certain frequency bands exhibit increased variability, suggesting a potential degradation of the signal due to model or data limitations. Globally, correlations between corrected GRACE signals and key climate indices largely diminish, confirming substantial removal of climate-related variability. However, the climate contribution to time-variable gravity beyond seasonal scales likely exceeds 30%, indicating incomplete correction and occasional alteration of residual signals that complicate the interpretation of deeper Earth processes. Despite these challenges, climate model-based correction shows promise for advancing source separation and deepening understanding of Earth’s interior dynamics via time-variable gravity data, contingent on future improvements in climate modelling.

SummaryFull-waveform inversion has been broadly adopted for acoustic and elastic media, but it lacks widely accepted methods for robust uncertainty quantification. This lack is in part due to an absence of assessment of proposed uncertainty quantification strategies. Here, we investigate four relatively inexpensive uncertainty estimation approaches based on truncated singular value decomposition of the inverse problem Hessian and its inverse. We numerically test these approaches across a range of parameter scales and application problems. We find that uncertainty estimates based on truncated singular value decomposition of the Hessian outperform those based on singular values of the inverse Hessian, due to both favorable singular value spectra of the former, and the greater ease of sampling the Hessian.

SummaryPalaeomagnetic studies of impact glasses offer valuable insights into their magnetization processes and thermal histories associated with impact cratering events. Australasian tektites are broadly distributed in the largest and youngest strewn field of the Cenozoic, and they provide a unique opportunity to investigate the intensity of Earth’s magnetic field around 788 ka and potential impact-induced magnetic fields. The northern part of the Australasian strewn field covers South China, and it corresponds to the uprange zone of the impactor’s trajectory. Magnetic properties of Australasian tektites in South China may contain unique information about this impact event, but their palaeomagnetic characteristics remain poorly constrained. Here, we report the first palaeointensity data of Australasian tektites sampled from the Early-Middle Pleistocene strata in South China. The results show that Muong Nong-type tektites recorded palaeointensities of 30 ± 8 μT, consistent with the geomagnetic field intensity around 780–790 ka. These findings suggest that around 788 ka, Earth’s magnetic field had partially recovered from the earlier intensity decline associated with the precursor event of the Matuyama–Brunhes reversal. By contrast, the splash-form tektites in South China are characterized by extremely weak natural remanent magnetization and unstable magnetization components, posing challenges for deriving reliable palaeointensity data. Although strong impact-induced remanent magnetization was not detected in the samples, this study demonstrates that Australasian tektites, particularly Muong Nong-type, are well suited for palaeomagnetic studies that may reveal potential impact-induced magnetization.

SummaryFluid-rock interactions in geological reservoirs can influence pore pressure and induce ground deformation at rates from millimeters to centimeters per year. Elastic deformation models often simplify structural heterogeneity that controls pore pressure and strain distributions, leading to inaccurate interpretations of reservoir properties from geodetic data. Here we investigate how depth-varying rock hydromechanical properties affect the magnitude, rate, and spatiotemporal characteristics of poroelastic deformation and pore pressure. Motivated by the Salton Sea geothermal field, we develop finite-element models of multilayered reservoirs to assess their transient and steady-state behavior in single-well fluid-extraction scenarios. These cases include (1) caprock-reservoir systems with varying permeability and caprock thickness, (2) compaction-induced porosity variations following Athy’s law, and (3) depth-dependent Young’s modulus. While uniformly lower porosity or permeability produces higher rates and earlier onset of deformation and pore-pressure changes, a less permeable or thicker caprock reduces vertical surface displacements, with pressure change reversals near the surface. Young’s modulus varying in alternating or linear profiles generally produces larger vertical displacements and non-monotonic displacement rate histories due to cross-layer fluid migration. Regarding spatiotemporal patterns, porosity or permeability decreasing with depth, or a thicker caprock, accelerates radial expansion of the deformation signal. In contrast, only layered mechanical properties can substantially alter the initial crossover distance and peak-value ratio between the vertical and radial surface displacements, indicating distinct impacts on deformation signatures. Our findings highlight the importance of accounting for structural heterogeneity in predicting and inferring the evolution of poroelastic processes in reservoir systems.

SummaryElectrical resistivity tomography (ERT) is used to infer the subsurface resistivity structure. ERT requires solving a nonlinear inverse problem that is often approximated as linear to reduce computational time. However, the approximation requires assumptions that cause limitations for the data analysis. Most of the computational time is due to the forward problem that requires solving the Poisson equation. Recently, similar forward problems have been shown to be replaceable with a surrogate model of lower computational cost. We present a geoelectric surrogate based on Fourier Neural Operators (FNO) and demonstrate a successful application in nonlinear inversion. The standard deviation of the surrogate prediction error for unseen samples is <5%. Furthermore, the surrogate reduces computational time by over three orders of magnitude per realization, enabling ERT for previously intractable settings. We apply the surrogate in Markov chain Monte Carlo (MCMC) inversion of simulated data. The results resolve sharp resistivity changes with plausible uncertainties.

SummaryThe Geocenter Motion (GCM) time series captures periodic variations arising from diverse Earth system changes. This study pioneers the use of Successive Variational Mode Decomposition (SVMD) in GCM research, enabling the precise extraction and analysis of these meaningful geophysical signals. SVMD outperformed Singular Spectrum Analysis (SSA) by effectively isolating signals and minimizing interference from components with similar variance contributions. However, a high maximum penalty factor in SVMD may lead to noise-dominated Intrinsic Mode Functions (IMFs). To overcome this limitation, we propose an extraction criterion that utilizes the standard deviation of the correlation coefficient and mean kurtosis as thresholds. Validations with simulations and the real GCM time series demonstrate its superiority over traditional single- and dual-threshold criteria, effectively retaining valuable information while excluding most noise-dominated IMFs. This improved approach is further employed to explore the geophysical driving factors of key periodic variations in the GCM time series, focusing on the annual, semi-annual, 10.5-year, 451-day, ∼160-day, and ∼120-day periods. Multi-source GCM analyses combined with the fingerprint method reveal distinct contributions from the Antarctic and Greenland ice sheets, terrestrial water storage, continental glaciers, and atmosphere-ocean interactions to different periodic signals. This study provides a robust methodology for decomposing GCM and attributing its variations to underlying Earth system changes, advancing our understanding and interpretation of global mass redistribution.

SummaryFull waveform inversion (FWI) is a popular method for subsurface parameter estimation. Despite its effectiveness in building high-resolution velocity models, the quality of the inversion result is significantly dependent on a fairly accurate, smooth initial model, which is often challenging to build. To weaken the influence related to the inaccurate initial model, we propose a deep learning (DL) matching-based FWI framework, namely DLM-FWI, where multiple convolution neural networks (CNNs) are used to construct an adaptive matching filter to better pinpoint the discrepancies between the synthetic and observed data. With the help of the CNN-based matching filter, the synthetic data will be regularized first, leading to intermediate data, and the model update will be conducted by minimizing the misfit between the intermediate and the observed data for improved data-fitting. More importantly, we integrate the whole inversion process into an automatic differentiation (AD) framework, simplifying the implementation of classic FWI. We apply the proposed DLM-FWI method to both synthetic and field datasets to validate its effectiveness. The results demonstrate that compared with classic FWI, DLM-FWI performs better in subsurface model reconstruction when the initial model is far from the global minimum.

SummaryUnderstanding subsurface fluid distribution in volcanic reservoirs is critical for geothermal energy development, critical mineral exploration, and forecasting eruptions. Here, we use travel-time tomography to image the seismic velocity structure beneath Aluto volcano, the first pilot geothermal project in Ethiopia, located in the Main Ethiopian Rift. Using seismic data recorded from January 2012 to January 2014, we invert for the 3D P-wave (Vp), S-wave (Vs), and Vp/Vs ratio. To reduce the non-uniqueness in interpretation, we also compare our results with previously published work on attenuation tomography and magnetotelluric images. Elevated Vp/Vs ratios (at 0 km below sea level (bsl)) around productive geothermal wells suggest high fluid content and/or elevated temperature. Vp/Vs values above 1.8 are observed along the caldera rims and hydrothermal vents, indicating fault and fracture systems as primary fluid conduits. High Vp/Vs below 6 km bsl likely reflects high-temperature areas or the presence of partial melt. In contrast, low Vp/Vs (<1.5), low Vp, and average to high Vs beneath the caldera at around 5 km bsl is interpreted as a crystallised body with over-pressurised gas volume formed during phase separation and transported upward through fractures and fault systems, accumulating at shallower levels. These findings highlight fluid pathways through the caldera rims and faults, with volatile-rich partial melt at greater depth beneath the caldera centre. Travel-time tomography thus offers a valuable constraints on subsurface fluid distribution and is valuable tool in geothermal exploration.

SummaryThis study presents a workflow to monitor spatiotemporal variations of the secondary microseisms using multi-array analysis. We employ ambient-noise cross-correlation beamforming (CC beamforming) across three dense seismic networks with different instrument responses: ANTICS in Albania (nodal-geophone and broadband), Hi-net in Japan (short-period), and SCSN (broadband) in Southern California. Independent of their instrumentation, these networks enable us to track the spatial and temporal evolution of secondary microseism sources in the northern Hemisphere from autumn 2022 to spring 2023. The workflow involves continuous data preprocessing for different instrumented sensors, ambient-noise cross-correlation, beamforming, and beam-power back-projection into a global map. We also propose sliding-window raw-data beamforming (RA beamforming) for the continuous broadband data in this workflow to record the absolute amplitudes of secondary microseisms recorded by ANTICS. Joint CC beamforming analysis across the three different networks improves the resolution of ambient-noise source localization and displays high consistency with the equivalent vertical force at the ocean floor. The results indicate that secondary microseism sources in the northern Hemisphere are predominantly driven by winter storms in the northern Atlantic and northern Pacific. The relative and absolute amplitudes of the beam-power for the northern Atlantic are also extracted from CC beamforming based on geophone sensors and RA beamforming based on broadband instruments from ANTICS, respectively. Both approaches provide robust estimates of microseism strength in the northern Atlantic, with CC beamforming displaying a higher correlation with the modeled ocean floor equivalent forces. This study confirms the feasibility of using cost-effective nodal seismic arrays for detailed monitoring of secondary microseisms and highlights the potential for integrating multi-array seismic data with oceanographic models for an improved understanding of seismic noise generation and propagation.

SummaryMoho topography model of subduction zones is an important component of deep tectonics and an important basis for verifying geodynamic processes. As one of the main factors affecting the accuracy of Moho topography model inverted by gravity method, the selection of inversion parameters suffers from the effect of nonlinear terms, which need to be reduced by constraints. Therefore, we applied the differential evolutionary algorithm to compute the inversion parameters and obtained a refined Moho topographic model of North Manila subduction on this basis. Synthetic tests show that the differential evolution algorithm can effectively mitigate the impact of nonlinear terms. With or without noise, the differential evolution algorithm is effective in finding better inversion parameters compared to the linear regression method. Particularly in Moho density contrast, the average value obtained from multiple runs of differential evolution algorithm still achieved a 54.4% improvement in accuracy. In practical application, the comparison results show that the RMS of the difference between this paper’s model and all seismic control points is 2.37 km with an improvement of at least 35.1%, which proves that this paper’s method is reliable. Furthermore, we examined the impact of various parameters on the method to validate its robustness.

SummaryThe far-field impact of Tibetan Plateau (TP) expansion on cratonic blocks remains enigmatic. We address this for the Ordos Block (OB) by constructing a high-resolution 3-D shear-wave velocity model using ambient noise tomography from an unprecedented dense seismic array (461 stations). Our model reveals: (1) NE-trending high-velocity anomalies at 10–25 km depth correlating with crustal magnetic signatures, providing seismic evidence for late Archean amalgamation of micro-blocks (Jining, Ordos, Xuchang, Xuhuai); (2) TP-induced reactivation manifesting as southwestern OB crustal thickening (50 km) with a high-velocity lower-crustal layer (≥4.0 km/s; 100 km wide, 10 km thick), attributed to TP lower-crustal underthrusting beyond the plateau margin (35.5°–37.5°N), facilitating >200 km strain transfer into the OB interior; (3) Incipient rifting dynamics in the Daihai Rift, where upper-crustal high-Vs (preserved rigidity) overlies mid-lower crust/uppermost mantle low-Vs anomalies (mantle-sourced thermal modification), indicating early-stage rifting driven by combined Pacific plate retreat and TP far-field stresses; (4) Craton-wide segmentation across a fundamental 37.5°N lithospheric boundary demarcating mantle upwelling/crust-mantle interaction (north) from passive TP push-dominated deformation (south). These findings redefine the OB as a strain-partitioned system, where lithospheric heritage controls differential response to plate-boundary forces.

SummaryThe Tanlu fault zone is an NNE-SSW oriented, large and deep strike-slip fault system running through eastern China. To investigate seismotectonics in and around the Tanlu fault zone, we adopt the LOC-FLOW approach to build a high-precision earthquake catalog. Our seismic data were recorded at 120 broadband TanluArray temporary stations with a sampling rate of 40 Hz and 76 broadband permanent stations with a sampling rate of 10 Hz from July 2019 to March 2023. We first conduct a series of experiments around the Luxi uplift and find that a higher sampling rate and a denser array of stations can significantly enhance the earthquake detection ability. Then we utilize both the temporary and permanent stations to conduct earthquake detection and location for the entire study area. As a result, 9648 earthquakes are detected and located in the REAL catalog, 6543 earthquakes in the HypoInverse catalog, and 5619 earthquakes in the HypoDD catalog, representing an increase of 20%, 22%, and 22%, respectively, as compared with the cases when only the TanluArray stations are used. Our location results show that earthquakes are mainly distributed in the Tanlu fault zone and active faults in relevant tectonic units. We collect 322 focal mechanism solutions (M > 2.0) of previous results from 2000 to 2020 to invert the stress field of the whole study region. The results show that the maximum principal stress axis of the whole study area is in the NEE-SWW direction, except that the Huoshan region is in the E-W direction. Along the Tanlu fault zone, the highest seismicity occurs in the Suqian-Weifang segment, and the Suqian seismic gap may be due to the aseismic slip along the fault planes of the Tanlu fault zone.

SummaryIn recent years, physics-informed neural networks (PINNs) have emerged as a powerful tool for seismic traveltime modeling and tomography. However, conventional PINNs do not consider applicable physical priors or quantify the uncertainty of the inverse problem, which is critical for reliable geological interpretation with topographical complexity. Thus, we propose a comprehensive PINN-based framework designed to tackle the critical challenge of inverting for the velocity in models with irregular topography, while also quantifying the inherent uncertainty. Leveraging automatic differentiation, our mesh-free approach directly accommodates complex surfaces without the need for specialized grids. To enhance inversion accuracy and physical consistency, we uniquely incorporate additional physical priors, namely well-log velocity profile and the principle of reciprocity. Furthermore, to address the non-uniqueness of the inverse problem, we integrate Monte Carlo (MC) Dropout to efficiently quantify model uncertainty without architectural modifications. Through 2D and 3D experiments on synthetic and real-world geological models, we demonstrate that our method accurately inverts for velocity structures with highly irregular topography. Results show that the inclusion of physical priors significantly improves model performance, while uncertainty quantification via MC dropout successfully highlights regions of higher uncertainty in the inverted velocity field, aligning with geological complexities in the velocities. This work establishes a robust and practical methodology for accurate and reliable seismic tomography in challenging geological settings.

SummaryShear-wave splitting measurement returns two parameters related to the fabric of the upper mantle: the orientation of the fast polarisation (fast direction), and a measure of the intensity and thickness of the fabric known as the split time. Spatial statistics of compiled splitting measurements indicate that the fast direction is spatially coherent, while the split time is not. We show, through modelling large numbers of noisy measurements, that single-earthquake splitting measurements exhibit a prominent upward bias in split time, the degree of which depends on specifics of the measurement process. Averaging single-event splitting parameters over many measurements does not mitigate this bias; however, stacking of error surfaces from individual measurements does, given sufficient back-azimuthal coverage, while also greatly reducing scatter. Published splitting results use a mix of these two averaging techniques, and this inconsistent bias between studies is likely responsible for the lack of spatial coherence in compiled split-time measurements. We demonstrate this in real data by examining a data set from Alberta, Canada and surrounding areas, for which a recent study published parameter-averaged results. By examining a comparable data set using error-surface stacking, we are able to greatly increase the coherence of the split times while obtaining highly similar fast directions. Our coherent split times are mapped to reveal a zone of strong splitting beneath the active Cordillera, and three zones of moderate to low split time within the cratonic lithosphere.

SummaryUnderstanding the spatial variability of ice content in frozen ground is key to design adequate measures to manage different ecosystems in frame of climate change. To-date investigations in frozen ground require the analysis of borehole data or the collection of multiple geophysical data. Here we propose the use of spectral induced polarization (SIP) as a technique that provides in quasi real-time information about changes in ice content in the subsurface. We demonstrate that exploring the frequency dependence in electrical conductivity and polarization (capacitive) properties in the frequency range between 0.1 and 75 Hz provides direct information about their relative variations in ice content. Our study is based on measurements conducted at nine representative permafrost sites in the European Alps with varying landforms and ice contents, including a pure ice and an unfrozen reference. We use the phase frequency effect (ϕFE) parameter as a parameter describing both the amplitude of the polarization and its frequency dependence to compare the response associated to the different sites. Our results show the lowest ϕFE in sites with low ice content, while increases in this parameter are associated with higher ice content. We evaluate the correlation between SIP parameters and validation ground ice contents for all sites and observe a clear correlation between ϕFE and volumetric ice content. The ϕFE results exhibit distinct landform-specific patterns, with the highest values found in rock glaciers, intermediate values in frozen talus slopes, lower values in bedrock permafrost, and the lowest in unfrozen talus slopes, reducing interpretation ambiguities in electrical resistivity results for assessing ice content.

SummaryTime-domain induced polarization (TDIP) data carry spectral information that can be used for petrophysical interpretation. At the same time, TDIP data can be collected in the field more efficiently than frequency-domain induced polarization (FDIP) data, thanks to the use of square-wave signals. However, TDIP field data are prone to noise, particularly strong near industrial installations and urban areas, above conductive media and in cases where little current is injected. The integral chargeability is a useful parameter to smoothen out the signal but it precludes any spectral interpretation. Debye decomposition (DD) is recognized as one of the best methods for spectral interpretation but the extracted parameters are particularly affected by data noise. More generally, processing TDIP data before further analysis, such as inversion or spectral analysis, is usually necessary for any quantitative interpretation. We propose here an automatic processing algorithm, based on the Kohlrausch-Williams-Watts (KWW) function, which is very close to the Havriliak-Negami model in frequency-domain, that fulfills this need. The processing is completed by an empirical handling of early-time electromagnetic coupling effects to improve the overall performance. The resulting procedure, tested and validated on three datasets that cover a large range of contexts, electrode configurations and acquisition settings, is available as open-source MATLAB scripts. The proposed approach is especially useful for further extracting spectral information from TDIP data through DD. Thanks to the theoretical framework offered by the KWW function, the behavior of the integral chargeability could be investigated in a systematic manner, using both synthetic and field TDIP data. Recommendations could be formulated on how to make use of the spectral information, while keeping the automatic processing transparent and accessible to unexperienced users. This work advances the use of TDIP in the field of environmental geophysics.

SummaryJoint petrophysical inversion is a powerful technique for using multiple geophysical modalities to estimate petrophysical or geotechnical parameters of the subsurface. A precise knowledge of the petrophysical laws for the full model domain is imperative to enable petrophysical coupling. In this work, we investigate the effect of partially invalid petrophysical laws on the inversion of a synthetic data set, using electrical resistivity tomography (ERT) and seismic traveltime data to image a CO2 plume in a Carbon Capture and Storage (CCS) setup. We consider a model consisting of a reservoir and a caprock in which only the reservoir can be described by a petrophysical law. We first apply a conventional (joint) petrophysical inversion (JPI) and show that the use of wrong petrophysical laws leads to systemic artefacts within the parts of the model in which the petrophysical relations are invalid. We then present a new hybrid partially petrophysically coupled joint inversion (P-JPI) approach that combines petrophysical coupling for regions with valid petrophysical laws, and structural coupling, whenever no reliable petrophysical laws are available. The P-JPI approach outperforms tomography based on the individual ERT or seismic data set, as well as joint structural inversion (JSI) based on the cross-gradient functional. The partially petrophysically coupled joint inversion thus enables petrophysical coupling and provides a unique, quantitatively interpretable saturation model for the CO2-plume. We further show that it is possible to detect zones with incorrect petrophysical relations by analysing the difference of the model updates based on the stand-alone data sets. Finally, we combine the detection of zones of incorrect petrophysical laws with the P-JPI to derive an inversion scheme that is independent of prior knowledge of the validity of petrophysical laws. Our novel methods facilitate direct estimation of the petrophysical subsurface parameters from multiple geophysical measurements if petrophysical relations are only available for parts of the model domain and provide means to quantify the spatial extent of regions where the petrophysical relations are valid.