Indirect joint petr oph ysical inv ersion of shallow-seismic and multi-offset ground-penetrating radar ﬁeld data

data simultaneously and provide consistent multiparameter models, which are hard to achieve by FWIs and individual petrophysical inversions. We also ﬁnd that the method is robust when there are uncertainties in petrophysical a priori information. Ov erall, the ﬁeld e xample prov es the great potential of using indirect JPI to solv e real-world problems.


I N T RO D U C T I O N
The near-surface area, a few tens of metres below free surface, are closely related to social development and life safety.A detailed characterization of this area is essential for urban construction, engineering exploration, environmental assessment, archaeological investigation, hydrological monitoring, polar research and so on (Everett 2013 ;Romero-Ruiz et al. 2018 ;Killingbeck et al. 2020 ).In near-surface surv e ys, geophysical techniques such as shallowseismic and ground-penetrating radar (GPR) methods are widely used (Doetsch et al. 2020 ;Leong & Zhu 2021 ;Liu et al. 2022a , b ).Shallow-seismic data are sensitive to the mechanical parameters in the subsurface but cannot identify the moisture distribution (Gassmann 1951 ).GPR data are highl y sensiti ve to the water content, but the depth of penetration is limited by the electrical conductivity (Annan 2005 ;Bradford & Deeds 2006 ).Individual inversions of these data may lead to inconsistent interpretations and not fully exploit their advantages.Combining the two data via joint inversion can provide complementary information for each inversion, reducing uncertainty and avoiding ambiguity (Linde et al. 2008 ;Domenzain et al. 2022 ;Moorkamp 2022 ;Huang et al. 2023 ).
The joint inversion methods can be divided into two main classes: joint structural inversion (JSI) and joint petrophysical inversion (JPI) (Abubakar et al. 2012 ;Linde & Doetsch 2016 ;Lan et al. 2018 ).JSI assumes that different geophysical parameters have similar spatial distributions (Gallardo & Meju 2011 ;Feng et al. 2017 ).The cross-gradient function is one of the most commonly used methods to quantify structural similarity (Gallardo & Meju 2003 ;Jordi et al. 2020 ).On the other hand, JPI supposes that different geophysical parameters are connected via petrophysical relations (Ghose & Slob 2006 ;Wagner et al. 2019 ).In general, JSI has a broader application than JPI because it does not strictly require accurate a priori petrophysical relations.Fur ther more, one can make petrophysical inferences based on JSI using methods like the scatter plot (Linde et al. 2006 ;Doetsch et al. 2010 ;Linde & Doetsch 2010 ).Ho wever , the structural constraint (soft link) is weaker than the petrophysical constraint (solid link) (Wagner & Uhlemann 2021 ).Therefore, JPI is also frequently used, especially in estimating petrophysical 974 C The Author(s) 2024.Published by Oxford University Press on behalf of The Royal Astronomical Society.This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( https://creati vecommons.org/licenses/by/4.0/ ), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
Indirect JPI of seismic and GPR field data 975 properties such as porosity, saturation and water content (Heincke et al. 2017 ;Wagner et al. 2019 ;Mollaret et al. 2020 ).
JPI is a valuable technique that links all geophysical parameters through their petrophysical relations.Ho wever , the applicability of JPI is limited by the accuracy of a priori information, including the assumed petrophysical relations and the rock matrix properties.To address this problem, we analysed the sensitivity of shallowseismic and GPR data to petrophysical parameters and proposed an indirect JPI method (Qin et al. 2022 ).The method indirectly uses petrophysical parametrization so that shallow-seismic and GPR data contribute to porosity and saturation updates, respecti vel y.Compared to the conventional JPI (Abubakar et al. 2012 ), indirect JPI increases the robustness of the inversion under imprecise a priori information and thus improves the applicability of the algorithm.We have validated the feasibility and robustness of the algorithm by synthetic examples with different acquisition geometries and wave types (Qin et al. 2022 ), but as of now it has not been used in real cases.
Indirect JPI w as de veloped from full-w av eform inv ersion (FWI) of shallow-seismic and GPR data.FWI is a multiparameter reconstruction technique that exploits the full information content of signals and provides high-resolution multiparameter subsurface models of geological interest (Tarantola 1984 ).FWI has been extensi vel y de veloped to solve problems at different scales, in different observation systems, and with different types of waves (Gao et al. 2021 ;Wang et al. 2021 ;Sun et al. 2023 ;Xu et al. 2023 ).In nearsurface scale, shallow-seismic wavefield is dominated by Love or Ra yleigh wa ves, which are the result of the interference of SH or P -SV w aves, respecti vel y.Unlike Rayleigh w av e FWI, Lov e wav e FWI has fewer model parameters, lower computational cost and less trade-offs between multiple parameters (Dokter et al. 2017 ;Wittkamp et al. 2019 ).Here, we apply indirect JPI for the first time to Lov e wav e and multi-offset surface GPR field data acquired in the Rheinstetten test site, Germany.We propose an input strategy to make the field application more operational.We compare indirect JPI with single inversions and assess our results with direct-push technology (DPT), borehole sample measurements and GPR migration result.The main goal of this paper is to e v aluate the applicability of indirect JPI to solve real-world problems.

Indirect JPI
Indirect JPI is a combination of shallow-seismic and GPR FWIs, where we use the traditional least-squares objective function to quantify the waveform misfit between the observed and synthetic data: and where S and EM are the objective functions of shallow-seismic and GPR FWIs, respecti vel y; the P -wav e v elocity V P , S -wav e v elocity V S and density ρ; m EM are the relative dielectric permittivity ε r and electrical conductivity σ (the magnetic permeability is constant to the value in vacuum and thus is not included).The petrophysical model parameters m P are porosity φ and saturation S w .The seismic objective function S is subject to the seismic petrophysical relation f S , and the GPR objective function EM is subject to the EM petrophysical relation f EM .In this study, the seismic petrophysical relation f S is Gassmann's equations (Gassmann 1951 ), and the EM petrophysical relation f EM is the complex refractive index model (CRIM) and Archie's equation (Archie et al. 1942 ;Birchak et al. 1974 ).
Indirect JPI attempts to integrate shallow-seismic and GPR data via separate objective functions (eqs 1 and 2 ).Like individual FWIs, indirect JPI iterati vel y updates the model parameters by minimizing the objective functions, for example using the preconditioned conjugate-gradient method (Qin 2022 ).We implement indirect JPI in each iteration through the following four steps, which we refer to as indirect petrophysical parametrization (Fig. 1 ).
In step 1, we use Gassmann's equations to calculate the P -wave v elocity V P , S -wav e v elocity V S and density ρ from the petrophysical parameters (Gassmann 1951 ;Abubakar et al. 2012 ):

976
T. Qin, T. Bohlen and Y. Pan where K and μ are the bulk and shear moduli of fluid-saturated rock, respecti vel y; ρ ma , ρ w and ρ a are the densities of the rock matrix, water and air, respecti vel y ( ρ w = 1000 kg m −3 , ρ a = 1.29 kg m −3 ); K ma , K f , K w and K a are the bulk moduli of the rock matrix, pore fluid, water and air, respecti vel y ( K w = 2.17 × 10 9 Pa, K a = 1.49× 10 5 Pa); M is the resulting average modulus and μ ma the shear modulus of the rock matrix; V P ma and V S ma are the P -wave velocity and S -wave velocity of the rock matrix, respecti vel y; β is the Biot's coefficient and φ c the critical porosity.We fix φ c to 0.4, above which the solid becomes a suspension (Nur 1992 ). Thus the range of values for the petrophysical parameters is 0 < φ < φ c and 0 < S w < 1.
In step 2, we use the seismic velocity parametrization in shallowseismic FWI to compute the gradients and update the models of V S and ρ (K öhn et al. 2012 ).The φ model is then transformed from the updated V S and ρ models by using Gassmann's equations as follows: (5) In step 3, we make φ the same size as the GPR models and then use it to calculate the relative dielectric permittivity ε r and electrical conducti vity σ b y CRIM and Archie's equation (Archie et al. 1942 ;Birchak et al. 1974 ;Day-Lewis et al. 2005 ;Bradford et al. 2009 ): where ε rma , ε rw and ε ra are the relative permittivities of the rock matrix, water and air, respectively ( ε rw = 81, ε ra = 1).σ w is the electrical conductivity of the groundwater, m is the cementation exponent and n is the saturation exponent.a is the tortuosity factor, and is used in practice to correct for data errors or poor petrophysical models, for example using Archie's equation when there is surface electrical conductivity (Linde et al. 2006 ;Revil 2013 ).Note that in eq. ( 6), if we fix σ w / a as a constant, dif ferent v alues of a can lead to the same model transformation and thus the same inversion result.
In this study, we simply set a to 1.
In step 4, we use the logarithmic parametrization in GPR FWI to compute the ε r gradient and update the ε r model (Meles et al. 2010 ).We then estimate S w by the well-reconstructed GPR model parameter ( ε r ) through eq. ( 7), which is another form of CRIM (Birchak et al. 1974 ): Finally, we make S w the same size as the seismic models and then use it in step 1 in the next iteration.
The main differences of indirect JPI compared to conventional JPI (Abubakar et al. 2012 ) are summarized as below: (i) We apply non-petrophysical parametrizations (seismic velocity parametrization and logarithmic parametrization) to compute the gradients and update the models of seismic and GPR parameters using their 'natural' model parameters.These parametrizations hav e prov en to be v ery efficient for shallow-seismic and GPR FWIs (Meles et al. 2010 ;K öhn et al. 2012 ).In contrast, conventional JPI uses petrophysical parametrizations where the gradients of porosity and saturation are given by the chain rule and include the gradients of P -wave velocity and electrical conductivity that are in low reliability and may lead to compromised reconstruction.
(ii) We take the sensitivity of data to seismic and GPR parameters into consideration.Previous studies have shown that shallowseismic FWI allows high-quality reconstruction of S -wave velocity (Pan et al. 2019 ), and permittivity is the parameter that can be most ef fecti vel y estimated by GPR FWI (Klotzsche et al. 2019 ).Therefore, w e update S -wa v e v elocity and density (density is required to calculate the shear modulus in eq. 5 ) by shallow-seismic FWI and ignore P -wave velocity, and update permittivity by GPR FWI and ignore electrical conductivity.
(iii) We take into account the sensitivity of seismic and GPR parameters to petrophysical parameters.The sensitivity analysis in Qin et al. ( 2022 ) suggests that the S -wav e v elocity and density are mainly af fected b y porosity, and the permitti vity has a relati vel y strong sensitivity to porosity and saturation (see Fig. 2 ).Thus we transform the S -wave velocity and density into porosity and transform the permittivity into saturation for efficient information exchange.Note that the last panel in Fig. 2 is a correction to figure 1 in Qin et al. ( 2022 ), where we forgot to transpose the model matrix when we plotted the electrical conductivity panel.Ho wever , this error does not affect the conclusions we have drawn in Qin et al. ( 2022 ) and in this paper.
(iv) For these hard-to-recover parameters (the P -wave velocity, density and electrical conductivity), we calculate them by petrophysical relations.Since porosity and saturation contain only the most reliable information from shallow-seismic and GPR data, the reconstruction of the P -wave velocity, density and electrical conductivity can be improved, which in turn helps to estimate petrophysical parameters and other seismic and GPR parameters.Thus, we construct a robust joint inversion framework for multiparameter reconstruction.
(v) We use separate objective functions rather than a combined objective function in the joint inversion to avoid calculating the data weighting matrix and the scaling factor (Heincke et al. 2017 ).In indirect JPI, the contribution of shallow-seismic and GPR data is automaticall y balanced b y the sensiti vity of the geophysical data to petrophysical parameters.This not only ensures that each data does its job but also reduces the trade-offs of multiparameter inversions.Additionally, this requires as little modification as possible from the individual FWIs to indirect JPI, making the programming of the joint inversion much easier.
In the implementation steps 2 and 4, we consider the sensitivity of shallow-seismic and GPR data to seismic and EM parameters and the sensitivity of seismic and EM parameters to petrophysical parameters.The combination of these sensitivities ensures that the information exchange between the two inversions is not disturbed by the weak-sensitive parameters and thus makes the joint inversion robust.

Input strategy
In petrophysical inversion, the seismic and EM petrophysical parameters of the rock matrix ( V P ma , V S ma , ρ ma , ε rma and σ w ) are normally assumed to be known, that is a priori information (For convenience, we treat the groundwater electrical conductivity σ w also as a rock matrix parameter).Ho wever , it is difficult to measure these parameters in the field or in the laboratory.Even more problematic, these parameters are often site-dependent, limiting the application of petrophysical inversion.To address these issues, we propose an input strategy to calculate the rock matrix parameters at the beginning of JPI (Fig. 1 ): with Eqs ( 8) and ( 9) are another form of the petrophysical relations f S and f EM when 0 < φ < φ c and 0 < S w < 1.This input strategy has three advantages: (i) It avoids measurements of rock matrix parameters and makes indirect JPI more operational in solving practical problems.Using this input strategy, the rock matrix parameters ( V P ma , V S ma , ρ ma , ε rma and σ w ) are implicitly included in the initial models ( V P , V S , ρ, ε r , σ , φ and S w , see Fig. 1 ).To build these initial models, many classical methods can be used, such as body-wave refraction tomography ( V P ), dispersion curve inversion ( V S and ρ), velocity analysis ( ε r ), amplitude attenuation estimate ( σ ), measuring the water content of soil samples ( φ and S w ; Xia et al. 1999 ;Annan 2005 ;Booth et al. 2010 ;Boiero & Socco 2014 ).Once the initial models are obtained, we use this input strategy to compute the rock matrix parameters and start the iteration of indirect JPI (Fig. 1 ).
(ii) It makes indirect JPI more robust (see details in Section 3.5 ).Since shallow-seismic and GPR FWIs rely heavily on seismic and GPR initial models ( V P , V S , ρ, ε r and σ ), it can be expected that indirect JPI can still yield similar results as long as the seismic and GPR models are credible.With this input strategy, empirical estimates based on field conditions can even be used for parameters that have no direct role in shallow-seismic and GPR FWIs, such as porosity, saturation, and Archie coefficients.Thus, indirect JPI can be applied like a single FWI with little limitation of petrophysical a priori information.
(iii) It makes indirect JPI also applicable for viscoelastic and dispersive EM media, which are closer to the reality.If the forward solver used to link d syn S and m S is viscoelastic equation (Bohlen 2002 ) or if the forward solver used to link d syn EM and m EM is Maxwell's equation in dispersive media (Bergmann et al. 1998 ), m S and m EM correspond to their values at the reference frequency (Fabien-Ouellet et al. 2017 ;Qin et al. 2023 ) and can be directly used to calculate the rock matrix parameters via eqs ( 8 ) and ( 9).
Note the input strategy also applies to the individual petrophysical inversions (IPIs).For seismic IPI, the rock matrix parameters ( V P ma , V S ma and ρ ma ) can be derived from the initial models ( V P , V S , ρ, φ and S w ) by eqs ( 8) and ( 10 ).For GPR IPI, the rock matrix parameters ( ε rma and σ w ) can be calculated from the initial models ( ε r , σ , φ and S w ) via eq.( 9).

Test site and data acquisition
The data have been acquired at the Rheinstetten test site, Germany, where a V-shaped trench called the Ettlinger Line was excavated in a sedimentary plain covered with gravel and sand from the Rhine river.This trench was refilled with sand a few decades ago and became invisible from the surface at the test site, which is a corner of the glider airfield (Fig. 3 a).The current ground layer is composed of partially saturated soil, and the groundwater table is below 6 m depth (Wittkamp et al. 2019 ).From the outcrop in the forest we found that the trench crosses the test site from the northwest to southeast.We therefore carried out a 2-D investigation with seismic and GPR profiles perpendicular to the Ettlinger Line (Fig. 3 a).
To record the Love wa ves, w e deployed 48 geophones (horizontal crossline component) from −3.5 to 43.5 m in the horizontal direction and used a hammer to blow on a steel beam source in the crossline direction (Figs 3 c and d).We acquired 12 seismograms with a shot spacing of 4 m and a fixed geophone spread.Our GPR data were recorded using a single channel GPR system with antennas of 200 MHz nominal centre frequency (Figs 3 b and c).We deployed the transmitter-receiver orientation in HH mode and acquired 18 radargrams with a transmitter spacing of 2 m.Unlike seismic data acquisition, we used a wide-angle reflection and refraction (WARR) method to acquire the GPR data where we fixed the transmitter and moved the receiver (mounted on a sled) toward or away from the transmitter.To track the receiver coordinates with centimetre-lev el accurac y, we used a real-time kinematic (RTK) positioning with a self-tracking total station (Boniger & Tronicke 2010 ).

Data preprocessing and inversion setup
We pre-process the raw data before using them in the inversion (Table 1 ).
(i) Since we use the 2-D finite-difference time-domain (FDTD) method to simulate seismic and EM wave propagation (Bohlen 2002 ;Irving & Knight 2006 ), we first resample the data to meet the time step requirements of the FDTD methods.This is done by transforming the data to the frequency domain, filling zeros and transforming back to the time domain.After resampling, the time Table 1.Shallow-seismic wave and multi-offset surface GPR data preprocessing steps.

1.
Data resampling in the frequency domain 2.
Interpolation of clipped direct-arri v al amplitudes 3.
DC-shift removal and dew o w 4.
Bad traces removal and offset limitation 6.
Data gridding in the time-offset domain 7.
3-D-to-2-D transformation (vi) After removing bad traces, the data may have irregular trace spacing.Fur ther more, in multi-offset GPR data acquisition, the worker moved the sled at an une ven w alking speed, which is another reason for the different trace spacing.To ensure a balanced illumination in the measurement area, we apply a 2-D spline interpolation in the time-offset domain with regular trace spacing.After data g ridding, seismog rams have 46-48 traces with a 1 m trace spacing, and radargrams have 131-192 traces with a 0.04 m trace spacing (Fig. 3 c).
(vii) The last preprocessing step before inversion is 3-D-to-2-D transformation which can transform the data recorded in the real world (3-D case) to the 2-D case.This step is important because we implement a 2-D inversion in this study.Based on the wave types contributing to the reconstruction, we use a hybrid 3-D-to-2-D transformation of the direct wave and single wave for shallowseismic data (Sch äfer et al. 2014 ), and a 3-D-to-2-D transformation of the reflected wave for GPR data (Forbriger et al. 2014 ).Note the latter may lead to some errors because the surface GPR data are often dominated by the direct waves (air wave and ground wave).
After preprocessing, the data are ready for the 2-D inversion.We then build initial models based on some data features such as dispersion, travel time and attenuation factor (Fig. 4 ).The S -wave velocity and density initial models are an average of dispersive curv e inv ersion results (Xia et al. 2012 ), and the near-ground values of the permittivity and electrical conductivity are estimated from the velocity and amplitude attenuation of the ground wave (Annan 2005 ).We assume that the porosity initial model decreases with depth and saturation is a homogeneous half-space (Fig. 5 ) based on the water content measurement from soil samples (see Section 3.4 ).We use the cementation exponent m = 1.4,and the saturation exponent n = 1.13 in Archie's equation where the low values of m and n are chosen for the possible existence of clay, a typical near-surface sediment (Abubakar et al. 2012 ).To account for the attenuation of the S -wav e v elocity, we use a viscoelastic solver with one relaxation mechanism of 40 Hz relaxation frequency and set the attenuation level Q S ≈ 13.3.The viscoelastic waves and EM waves are modelled by the 2-D FDTD method (Bohlen 2002 ;Irving & Knight 2006 ).We use similar model space but different grid sizes for seismic (0.12 m × 0.12 m) and GPR models (0.04 m × 0.04 m) to simulate them accurately.The convolutional perfectly matched layer (CPML) is included at the model boundaries, except for the free surface of seismic models where an imaging method is applied (Le v ander 1988 ).In GPR models, the air layer of 1 m thickness remains constant during the inversion and thus is not displayed (Fig. 4 b).
We use a multiscale strategy to avoid cycle skipping in the inversion (Bunks et al. 1995 ).We choose five inversion stages to sequentially use data with ever decreasing wavelength.From the first stage to the fifth stage, frequency bands vary from 5 to 20, 35, 45, 60 and 80 Hz for seismic inversion, and frequency bands vary from 5 to 30, 40, 50, 70 and 100 MHz for GPR inversion.At the beginning of each inversion stage, we estimate the source wavelets with a Wiener filter (Groos et al. 2014 ).The abort criterion is that the relative data misfit change is less than 1 per cent, and the maximum number of iterations per stage is 15.In the joint inversion, the program can switch to the next stage if the abort criteria of both individual inversions are satisfied or if the maximum iteration number is reached.In addition, we apply a 1-D Gaussian filter in the horizontal direction to the gradient to suppress the artefacts shorter than the dominant wavelength.

Inversion results
The seismic, GPR and petrophysical models reconstructed by indirect JPI successfully reveal the presence of the Ettlinger Line, shown as a triangle anomaly with low S -wave velocity, low density, high permittivity and high electrical conductivity values (Fig. 4 ), resulting from the high porosity and saturation values (Fig. 5 ).On the one hand, the S -wave velocity result is comparable to that of the 3-D shallow-seismic FWIs of Pan et al. ( 2021 ).On the other hand, the permittivity result is in high agreement with the GPR migration image of Wittkamp et al. ( 2019 ) and Qin ( 2022 ).Due to the constraint of petrophysical relations, the density model also reveals the exact shape of the trench, which is difficult to see from past investigations (Wittkamp et al. 2019 ;Pan et al. 2021 ).For the same reason, the electrical conductivity model has a similar structure to the permittivity model.Note that in the GPR models, the boundaries of the trench become less visible compared to the high permittivity and high electrical conductivity anomalies in the interior.It may result from the high electrical conductivity environment near the surface that degrades the penetration depth of the GPR signal, especially at the bottom of the trench.
Indirect JPI outperforms individual FWIs (Fig. 4 ).Individual FWIs use the same objective functions (eqs 1 and 2 ) but are not subject to petrophysical relations.In the seismic model reconstruction, the S -wave velocity models estimated by all inversions are very comparab le, w hile the density models reconstructed by seismic FWI show a high-density anomaly inside the trench and a low-density anomaly to the left.Based on the latter DPT and borehole measurements in Section 3.4 , these density anomalies could be artefacts caused by the crosstalk from the S -wav e v elocity or by the low sensitivity of Lov e wav e data to density.In comparison, indirect JPI provides significant improvements where there is a low-density triangle anomaly in the middle.When reconstructing the GPR models, individual GPR FWI outlines the shape of the trench, but the model update focuses mainly on the near-surface region (depth shallower than 2 m depth, Fig. 4 b).Indirect JPI can update deeper areas and show clearer interfaces.
We compare indirect JPI with IPIs where seismic IPI follows eq. ( 1), and GPR IPI follows eqs ( 2 ).IPIs use petrophysical parametrization to directly update the petrophysical models and then calculate the geophysical models based on the petrophysical relations.Seismic IPI and seismic FWI reconstruct similar S -wave velocity models, which implies that Gassmann's equations we used are close to the reality (Fig. 4 a).The density model benefits from the petrophysical constraint and shows structures similar to the Swav e v elocity model.Howev er, the resolution of the density model in seismic IPI is lower than in indirect JPI as the latter is improved by higher resolution GPR data.In terms of GPR IPI (Fig. 4 b), the reconstructed models present similar subsurface structures to that of GPR FWI, meaning that CRIM and Archie's equation are also applicable to this test site.Since surface GPR data are dominated by the shor t-wavelength infor mation, the long-wavelength backg round models are difficult to update by GPR FWI or IPI alone.To overcome this dra wback, joint in versions attempt to use the complementary information from Love wave data.In indirect JPI, seismic data successfull y of fers the needed information for GPR FWI through the porosity model, thus allowing reconstruction of the low permittivity and low electrical conductivity background.Hence, compared to individual inversions, indirect JPI using more data can constrain the reconstruction process better.
Indirect JPI provides consistent petrophysical models in seismic and GPR in versions.P etrophysical parametrization used in IPIs  contains weak sensitivity information, for example electrical conductivity gradient in GPR IPI.Consequently, seismic IPI and GPR IPI generate conflicting results at 0.5-1.5 m depth on the right side of the trench (Fig. 5 ).Seismic IPI describes this region as high porosity and high saturation anomalies, while GPR IPI interprets it as low porosity and low saturation anomalies.However, indirect JPI characterizes this region as a water-poor layer consisting of high porosity and low saturation anomalies, which proves to be more reliable because in this case we can match both shallow-seismic and GPR data (see Fig. 6 ).
The seismic and GPR objective functions converge at 35 iterations in indirect JPI.For convenience, we use in the following the relative data misfit, that is dividing the data misfit by that obtained from the initial models (the relative data misfit of the initial models is one).To calculate the data misfit of the initial models, we simulate the synthetic data on the initial models with the source wavelet estimated at the last inversion stage and apply the last frequency bandpass filter to the observed data.The GPR data misfit (the relative data misfit is 0.9208) is higher than seismic data misfit (the relative data misfit is 0.8017) due to four reasons.(1) EM waves attenuate much more than seismic waves at the same propagation distance (see figure A2 4) GPR initial models are quite simple and estimated from the ground wave only, while seismic initial models reference previous studies, which have proven to be useful (Pan et al. 2019 ;Wittkamp et al. 2019 ;Gao et al. 2020 ;Irnaka et al. 2022 ).Overall, the field data fitting could not reach the same good level as the synthetic examples in Qin et al. ( 2022 ).On the one hand, this is due to the complexity of the field data, lower SNR, higher attenuation, and the influence of the test site, equipment and operators; on the other hand, it is a result of data preprocessing errors, simple initial models and limitations of our forward solvers.
The seismogram fitting of seismic IPI (the relative data misfit is 0.7912) is slightly better than that of indirect JPI (Fig. 6 a), and the radargram fitting of GPR IPI (the relative data misfit is 0.9163) and indirect JPI are comparable (Fig. 6 b).Overall, the data misfits of indirect JPI are slightly higher than that of seismic and GPR IPIs (about 1 and 0.4 per cent, respecti vel y) because of the interaction of seismic and GPR data and the additional petrophysical constraint.Ho wever , neither seismic IPI nor GPR IPI can minimize the two data misfits simultaneously.For example, we use the petrophysical models reconstructed by seismic IPI to derive GPR models (eq.6 ) and perform a GPR forward modelling on them; we then use the simulated GPR data to calculate the data misfit and find that it is much higher than that of indirect JPI, and vice versa.From this perspective, indirect JPI is a better way to fit seismogram and radargram at the same time, leading the seismic and GPR model update in an acceptable direction for both.

Assessment
To e v aluate our interpretations of the density model, we compare the inversion results with two independent DPT measurements (see Figs 4 a, 7 a and b).In DPT measurements, a metal pile was hit by a free-falling slide hammer (10 kg) and thus pushed into the ground.By recording the number of hits per 0.1 m depth pushed in, we measured the consolidation degree in the subsurface.The higher the number, the more compact the soil.We have one DPT measurement outside the trench (DPT1) and another one inside the trench (DPT2).These DPT measurements reveal a change from loose topsoil (lower density) to compacted subsoil (higher density) at 1 m depth at the DPT1 position and an interface between refilled sand (lower density) and underlying soil (higher density) at 3 m depth at the DPT2 location, which corresponds to the bottom of the trench (Figs 7 a and b).These findings agree with the density model of seismic IPI and indirect JPI and prove that petrophysical constraint improves the density reconstruction.
To assess the water content given by petrophysical inversions, we drilled two boreholes (BH1 and BH2) close to the DPT locations and collected soil samples every 0.5 m depth (Fig. 5 ).We then measured the gravimetric water content θ g of the soil samples by drying them in the laboratory.The borehole measurements show a change in water content at similar depths to the DPT results.For con venience, we con vert the petrophysical inversion result to the gravimetric water content by θ g = θ/(1 − θ), where θ = θ v / ρ and the volumetric water content θ v = φS w .In GPR IPI, we use the density of seismic IPI to calculate θ g .Note that the uncertainties present in the density may affect the conversion from θ v to θ g .We find that indirect JPI fits well with the BH1 measurement in both high and low water content area (Figs 5 and 7 c).For BH2, seismic IPI overestimates the water content in the trench, and GPR IPI underestimates it and shows result very close to the initial model  (Figs 5 and 7 d).Indirect JPI's result is better than seismic IPI's result, although further improvements are needed.
We compare the indirect JPI results with the GPR migration result along the same profile (Fig. 8 ).After performing a classical GPR processing including bandpass filtering (50-400 MHz), zero-time correction, time-dependent scaling, we apply a Kirchhoff migration routine (Moran et al. 2000 ) to the common-offset data (offset = 0.5 m) using a constant subsurface velocity of 0.1 m ns −1 .The migrated GPR image shows the trench in the centre and the dipping high amplitude reflector on the right-hand side of the trench (Fig. 8 a).The migrated data on the left-hand side of the trench appear more chaotic, with a lower penetration depth.For better comparison, w e overla y the migration image on the indirect JPI results (Figs 8 b and c).We find that the migration image agrees well with the porosity and saturation images in terms of the location and structures of the Ettlinger Line and the right-hand reflector.Therefore, we conclude that indirect JPI is capable of reconstructing the petrophysical models with high reliability, especially for shallow subsurface structures.

Robustness tests
It is well known that the performance of petrophysical inversions tends to depend on a priori information, that is the rock matrix parameters and Archie's coefficients in this study.To e v aluate the inversion robustness, we take the mean structural similarity (MSSIM) index to measure the fidelity of the reconstructed models of indirect JPI using different petrophysical initial models and Archie's coef ficients relati ve to the reference.We adopt the same settings as Boniger & Tronicke ( 2010 ) to compute the MSSIM index.The closer the MSSIM index is to one, the more similar the two compared objects are.Instead of giving the rock matrix parameters explicitly, we compute them from the initial models via the input strategy (eqs 8 and 9 ).Therefore if we fix seismic and GPR initial models, the rock matrix parameters can be changed either by petrophysical initial models (porosity and saturation) or by Archie's coefficients.
We take the indirect JPI shown in Fig. 5 as the reference and other four indirect JPIs using homogeneous petrophysical initial models in half-space for comparison (Fig. 9 ).When we change petrophysical initial models to different values ( φ = 0.2 or 0.3 and S w = 0.5 or 0.7), the reconstructed seismic and GPR models present a high deg ree of str uctural similarity to the reference (see MSSIM index in Fig. 9 ).Compared to the GPR models, the seismic models are reconstr ucted ver y stably (all MSSIM indices are greater than 0.985) because they are influenced only by porosity, whereas the GPR models are affected by both porosity and saturation.The results of φ = 0.3 generally have higher structural similarity than those of φ = 0.2 because the near-ground values of the reference porosity model are close to 0.3 and shallow-seismic and GPR data are more sensitive to shallow than deep zones.With φ = 0.2 and S w = 0.7, we observe the largest differences in permittivity (MSSIM = 0.9355) and electrical conductivity (MSSIM = 0.8552).Ho wever , even in this case, all models can delineate consistent subsurface structures and provide a meaningful geological interpretation.Note that there are significant differences between the petrophysical models reconstructed with different petrophysical initial models.Among them, the reference model best matches the gravimetric water content of the borehole soil samples (Figs 7 c and d) because its initial model allows for the trend of decreasing water content with depth (Fig. 5 ).Therefore, we recommend building the petrophysical initial models based on the borehole data if they are available.If not available, an empirical guess of the petrophysical initial models is acceptable for indirect JPI, but the geological interpretation should be based primarily on the reconstructed seismic and GPR models.
We test indirect JPI with the same initial electrical conductivity model (Fig. 4 b) but different Archie's coefficients, corresponding to dif ferent groundw ater electrical conducti vities (eq. 9 ).The range of m between 0.4 and 2.4 and n between 1.13 and 3.00 gives almost identical results (MSSIM > 0.99, see Fig. 10 ).In other words, the choice of Archie's coefficients has a negligible impact on indirect JPI's performance.This should be attributed to the abandonment of weakly sensitive electrical conductivity information in indirect petrophysical parametrization and to the indirect calculation of the groundwater electrical conductivity in the input strategy (Fig. 1 ).
Ov erall, these e xamples demonstrate that indirect JPI does not need to know accurate a priori information.With suitable seismic and GPR initial models, any petrophysical initial models and Archie's coefficients in a reasonable range can produce similar results.This strength makes indirect JPI a robust and promising technique that could be easily applied to other field environments with as few assumptions about petrophysical relations as possible.

New scientific findings
This is a follo w-up w ork to Qin et al. ( 2022 ) that performed the theoretical study and synthetic test on indirect JPI.  in solving practical problems.To facilitate the field application, we combined indirect JPI with an input strategy that avoids the measurement of the rock matrix parameters, making the algorithm more robust and suitable for complicated environments.Comparison with the three third-party methods showed that the results of indirect JPI are of higher quality than those of FWI and IPI.Note that the robustness tests are different from those in Qin et al. ( 2022 ), where we keep the petrophysical initial models unchanged and perturb the rock matrix parameters (or Archie's coefficients), causing changes in the seismic and GPR initial models and thus having a greater impact on the reconstruction results.In this study, we leave the seismic and GPR initial models fixed, but adopt the input strategy to change the rock matrix parameters by petrophysical initial models or by Archie's coefficients.Shallow-seismic and GPR FWI mainly depend on the seismic and GPR initial models, so indirect JPI in this study suffers less from inaccurate petrophysical a priori information than in the previous study.
Indirect JPI is an algorithm different from conventional JPI.Conventional JPI is a combination of two IPIs in which the petrophysical parameters are reconstructed by the chain rule based on petrophysical relations (Abubakar et al. 2012 ).Conventional JPI merges all the information from two geophysical data and tends to obtain a compromise between the two IPIs.In contrast, indirect JPI integrates two FWIs to update seismic and GPR parameters, and then indirectly transforms the results into petrophysical parameters using petrophysical relations.This is why we call it indirect petrophysical parametrization.The stable performance of the algorithm should be attributed to the rational use of petrophysical sensitivity, non-petrophysical parametrizations and the input strategy (Fig. 1 ).Taking into account the sensiti vity dif ferences of geophysical data to petrophysical parameters ensures a highly efficient exchange of information between two geophysical inversions.Seismic velocity parametrization and logarithmic parametrization fully exploit the reconstruction potential of shallow-seismic and GPR FWI, respecti vel y.The input strategy reduces the difficulty of using petrophysical inversion and makes indirect JPI less dependent on petrophysical initial models and petrophysical relations.The integration of the above three parties makes indirect JPI a robust algorithm that may be applied to the field data.Hence, indirect JPI overcomes the bottleneck of conventional JPI and provides a new way for JPI development.
The moti v ation of this study is to answer the questions about geophysical inversion: Can we use JPI to improve the reconstruction result and avoid conflicting interpretations?Is it possible to easil y appl y JPI to field data?The first question is a basic requirement for all joint inversion studies, and our answer is yes.
The field example revealed that indirect JPI can better use the strengths of shallow-seismic and GPR data and produce results much better than individual inversions.The second question is challenging for petrophysical inversion because the assumed petrophysical relations and a priori information may be site-dependent and hard to obtain.Ho wever , our study indicated that indirect JPI still has a stable performance when there are uncertainties in a priori information.Indirect petrophysical parallelization, in combination with the input strategy, makes the field data application of JPI possible and promising.Therefore, our answer to the second question is yes, and we believe indirect JPI deserves more attention.

Limitations of indirect JPI
The petrophysical parameters we consider in this study are porosity and saturation (water content is their product), which are the objects of most petrophysical inversions.In some cases, it may be necessary to include other petrophysical parameters, for example clay content.Clay content is a key factor affecting seismic and EM properties.For instance, electrical conductivity is controlled by clay content and groundwater electrical conductivity.Although we attempted to use low values of m and n to account for the presence of clay, the adopted Archie's equations may not be the best choice because it w as deri ved from a clay-free matrix (Archie et al. 1942 ).A possible solution is to use petrophysical relations that can consider the effect of clay content on seismic and EM parameters, such as the Voigt-Reuss-Hill (VRH) boundary model and Linde's modification of Archie's equation (Linde et al. 2006 ;Hu et al. 2021 ).Additionally, electrical resistivity data are sensitive to electrical conductivity and may be introduced into the joint inversion to reconstruct the clay content model.In this way, indirect JPI can be extended from two to three geophysical data based on the sensitivity links between shallow-seismic data and porosity, GPR data and saturation, and electrical resistivity data and clay content.
Non-petrophysical parametrizations is another key factor affecting indirect JPI's performance.Seismic velocity parametrization has proven to be very beneficial for seismic multiparameter imaging, mainly due to the similarity in magnitude of seismic velocity and density (K öhn et al. 2012 ).Nevertheless, seismic modulus parametrization may also be worth a try since we actually use the shear modulus to calculate porosity (eq. 5 ).In seismic velocity parametrization, the shear modulus is calculated from the S -wave velocity and density, and may be contaminated by less reliable density information.Whereas, in seismic modulus parametrization, the shear modulus is updated onl y b y the shear modulus gradient without any interferences, and thus may be more reliable.Therefore, indirect JPI with seismic modulus parametrization instead of seismic velocity parametrization may yield better results, which is subject to further validation.
As the first step, we e v aluated our algorithm via Lov e wav e and multi-offset surface GPR field data.As shown in table 2 in Qin et al. ( 2022 ), there are 12 combinations of geometries (crosshole and surface) and wave types ( P -SV, SH, TE and TM) in the 2-D case.Therefore, other combinations need to be tested in the next studies.What we can expect is that indirect JPI with P -SV wave (Ra yleigh wa ve) will be more challenging since more parameters have to be reconstructed, and the crosshole geometry will illuminate better in deeper area than surface geometry.Applying indirect JPI in the 3-D case will be the next step.
The current version of indirect JPI does not include model regularization.Model regularization is an ef fecti ve w ay to stabilize the inversion process and improve the imaging results (Tikhonov & Arsenin 1977 ).For example, the maximum smoothness regularization (Constable et al. 1987 ) may constrain the horizontal roughness and suppress artefacts in the near-surface GPR models (Fig. 4 b); the total variation regularization (Feng et al. 2019 ) may make the trench boundaries and other subsurface interfaces in the seismic models sharper (Fig. 4 a).Therefore, indirect JPI combined with model regularization is worth further investigation.

Implications of the work
The w orkflo w we propose in Fig. 1 can be regarded as a general framework of joint inversion.It is a way to stabilize the reconstruction process and exchange the most reliable information from each data.Our study moti v ates the consideration of sensitivity difference and the use of non-petrophysical parametrization in the joint inversion.In near-surface survey, indirect JPI of shallowseismic and electrical resistivity data can be developed similarly by replacing GPR FWI with electrical resistivity inversion and calculating saturation from electrical conductivity based on Archie's equation (Wagner et al. 2019 ).GPR and electrical resistivity data can be integrated by indirect JPI as well (Domenzain et al. 2020 ), where the frequency difference in the electrical conductivity reconstructed by two inversions should be taken into account (Qin et al. 2023 ).In oil exploration, one can perform indirect JPI of seismic and controlled-source electromagnetic measurements (CSEM) data, namely, replacing GPR FWI with CSEM inversion and computing saturation from electrical conductivity (Abubakar et al. 2012 ).Note that it is possible to calculate porosity from P -wave velocity since seismic FWI uses mainly reflected waves and is also sensitive to P -wave velocity in oil exploration.Fur ther more, JPI is not limited to the traditional framework of using only petrophysical parametrizations.Our study indicates that using non-petrophysical parametrization in JPI can better exploit the reconstruction potential of FWI.Thus, indirect JPI holds great promise as a generalized framework that may be applied to a wide range of joint inversion studies.
The approach presented is not limited to the empirical models gi ven b y Gassmann' s equations, CRIM and Archie' s equation.It is also applicable to other petrophysical relations.For example, the time-averaging equation can be used for the seismic petrophysical relation f S (Wagner et al. 2019 ), and Topp's equation and Ewing's modification of Archie's equation can be taken as the EM petrophysical relation f EM (Topp et al. 1980 ;Ewing & Hunt 2006 ).Some petrophysical relations reveal similar sensitivity trend as Fig. 2 , for example the permittivity increases with porosity and saturation in CRIM and Topp's equations.Indirect JPI using these petrophysical relations is expected to give similar results as in this study.For petrophysical relations with different sensitivity trends, such as Gassmann's equations and the time-averaging equation, indirect JPI ma y ha ve a stable performance but different results.How ever, how to judge whether the adopted relations are valid for the test site is an open question.As an attempt, we suggest to implement FWI and IPI first.If their results are similar (e.g. S -wave velocity in Fig. 4 ), the petrophysical relations are close to the real situation in the test site, and then we can run indirect JPI using these petrophysical relations.If the results of FWI and IPI are far apart, other petrophysical relations should be tested until the most suitable one is found.An alternative is to use JSI to derive possible petrophysical relations from a scatter plot of multiple geophysical parameters (Linde & Doetsch 2010 ), for example S -wave velocity and permitti vity reconstructed b y seismic and GPR FWIs.This deser ves fur ther investigation.
The computational cost of the joint inversion is the sum of two individual inversions.In this paper, we run our algorithm on a 36core computer cluster.Shallow-seismic FWI has 12 seismograms (2500 sampling rate) and model space 412 × 50 grids, simulated with 12 source parallelization and three model domain parallelization (one model decomposed into three subvolumes) (Bohlen 2002 ).GPR FWI has 18 radargrams (2048 sampling rate) and model space 1130 × 175 grids, simulated with 18 source parallelization and two model domain parallelization.We use a subset FWI (SFWI) method to speed up GPR FWI by a factor of five and reduce the memory usage to 20 per cent (Qin 2022 ).As a result, the total computational time of indirect JPI is about nine minutes.The computational cost can be further reduced by applying the source encoding method to shallow-seismic FWI and GPR FWI (Krebs et al. 2009 ;Feng et al. 2023 ).Ho wever , GPR data acquisition requires special attention because each data has an unfixed geometry spread.To make source encoding applicable to GPR data acquisition and to eliminate the crosstalk of uncorrelated wav efields, frequenc ydivision encoding method may be considered (Huang & Schuster 2012 ).

C O N C L U S I O N
In this paper, we applied indirect JPI to shallow-seismic and multioffset surface GPR field data for consistent imaging of the nearsurface targets.Indirect JPI exploits only highly sensitive relations between geophysical and petrophysical parameters.In this case these relations are S -wav e v elocity ∝ porosity and permittivity ∝ saturation.We proposed an input strategy to make indirect JPI more operational and robust in solving practical problems.The application at the Rheinstetten test site showed that this approach not only outperforms FWIs and IPIs in reconstructing seismic and GPR parameters, but also provides more consistent petrophysical models than IPIs.Indirect JPI presented significant improvements in estimating saturation, density, and electrical conductivity, therefore reducing the ambiguity and uncertainty of single geophysical techniques and facilitating the final geological interpretations, such as determining groundwater distribution and facies stratification.This study also suggested that, due to the use of indirect petrophysical parametrization, this approach can ef ficientl y exchange the high confidence information given by each inversion.Thanks to the separate contributions of seismic and GPR data, indirect JPI reduced the reliance on a priori information and exhibited great potential for

Figure 1 .
Figure 1.The w orkflo ws of indirect joint petrophysical inversion (JPI) with an input strategy.The implementation steps of indirect petrophysical parametrization are indicated by numbers one to four.

Figure 2 .
Figure 2. The seismic model parameters ( V P , V S and ρ) and electromagnetic (EM) model parameters ( ε r and σ ) as functions of porosity φ and saturation S w .The rock matrix parameters and Archie's coefficients are the same as Qin et al. ( 2022 ).Note that the last panel (electrical conductivity) is a correction to figure 1 in Qin et al. ( 2022 ).

Figure 3 .
Figure 3. (a) Map of the Rheinstetten test site, where the translucent grey area shows a target trench, the Ettlinger Line, and the black line shows the ground-penetrating radar (GPR) and shallow-seismic surv e y lines.(b) GPR field data acquisition with the receiver mounted on a sled and tracked by a real-time kinematic (RTK) positioning.(c) Observation geometry of 18 multi-offset GPR data and 12 fixed-spread shallow-seismic data, where the translucent grey area shows the trench location.(d) Shallow-seismic data acquisition using a hammer source and fixed geophone locations.

Figure 4 .
Figure 4. (a) Seismic models ( S -wave velocity V S and density ρ) and (b) GPR models (relative permittivity ε r and electrical conductivity σ ).The four columns are the initial models, the reconstructed models of full-waveform inversion (FWI), the reconstructed models of individual petrophysical inversion (IPI), and the reconstructed models of indirect JPI, respecti vel y.In the initial models, the red stars are the sources and the dashed triangle outlines the cross-section of the Ettlinger Line.The borehole histograms overlaid on the density model are the direct-push technology (DPT) results (Figs 7 a and b), where the transparent and translucent areas present the unconsolidated and consolidated soil, respecti vel y.

Figur e 5 .
Figur e 5. Petroph ysical models (porosity φ and saturation S w ) and volumetric water content model θ v .The four columns are the initial models, the reconstructed models of seismic IPI, the reconstructed models of GPR IPI and the reconstructed models of indirect JPI, respecti vel y.The borehole histo grams overlaid on the water content model are the gravimetric water content θ g given by borehole soil samples (Figs 7 c and d), where the transparent and translucent areas exhibit high and low water content, respectively.

Figure 6 .
Figure 6.Comparison of the observed data and the synthetic data corresponding to the results of IPI and indirect JPI.(a) The horizontal velocity seismogram of the 1st shot, shown once every six traces.(b) The horizontal electric field radargram of the 16th GPR transmitter , sho wn once every twenty traces.The rectangular windows display the zoomed waveforms.For better visualization, the data are normalized trace by trace.
in Qin et al. 2022 ), thus GPR inversion mainl y fits near-of fset data, where errors may be greater since the 3-D-to-2-D transformation uses a far-field approximation (Forbriger et al. 2014 ).(2) It is difficult to fit the direct waves in GPR inversions because we use the 3-D-to-2-D transformation of the reflected wav e. (3) Our 2-D solv er cannot simulate the radiation patterns and antenna-ground coupling well.(

Figure 7 .
Figure 7.Comparison of inversion results, DPT, and borehole measurements.The black lines in (a) and (b) are the number of hits measured by DPT, and the black lines in (c) and (d) are the gravimetric water contents ( θ g ) of borehole soil samples.The colourful lines are the inversion results, that is density in (a) and (b) and θ g in (c) and (d).The gray dashed lines mark the interfaces between transparent and translucent areas of the borehole histograms shown in Figs 4 and 5 .

Figure 8 .
Figure 8.Comparison between (a) the migration result of common-offset GPR data, (b) the porosity and (c) saturation results of indirect JPI.The colour scale of (a) is clipped to 20 per cent of the highest amplitude for better visualization.In (b) and (c) we overlay the migration image on the petrophysical models.

Figure 9 .
Figure 9. (a) Seismic models ( S -wave velocity V S and density ρ), (b) GPR models (relative permittivity ε r and electrical conductivity σ ) and (c) petrophysical models (porosity φ and saturation S w ).The five columns are the indirect JPI results with different petrophysical initial models (see the title of each row).These five JPIs have the same seismic and GPR initial models as in Fig. 4 .The referenced petrophysical initial models are the same as in Fig. 5 .The values overlaid on the models in (a) and (b) are the mean structural similarity (MSSIM) index relative the reference.

Figur e 10 .
Figur e 10.Petroph ysical models ( φ and S w ).The five columns are the indirect JPI results with different Archie's coefficients (cementation exponent m and saturation exponent n ).The initial models of five JPIs are the same as the reference in Fig. 9 .The values overlaid on the models are the MSSIM index relative to the reference ( m = 1.4 and n = 1.13).
The major novelty of this work is that we assessed the feasibility of our algorithm Downloaded from https://academic.oup.com/gji/article/237/2/974/7623607 by KIT Library user on 22 April 2024