# Joint inversion of NMR and SIP data to estimate pore size distribution of geomaterials

Joint inversion of NMR and SIP data to estimate pore size distribution of geomaterials Summary There are growing interests in using geophysical tools to characterize the microstructure of geomaterials because of the non-invasive nature and the applicability in field. In these applications, multiple types of geophysical data sets are usually processed separately, which may be inadequate to constrain the key feature of target variables. Therefore, simultaneous processing of multiple data sets could potentially improve the resolution. In this study, we propose a method to estimate pore size distribution by joint inversion of nuclear magnetic resonance (NMR) T2 relaxation and spectral induced polarization (SIP) spectra. The petrophysical relation between NMR T2 relaxation time and SIP relaxation time is incorporated in a nonlinear least squares problem formulation, which is solved using Gauss–Newton method. The joint inversion scheme is applied to a synthetic sample and a Berea sandstone sample. The jointly estimated pore size distributions are very close to the true model and results from other experimental method. Even when the knowledge of the petrophysical models of the sample is incomplete, the joint inversion can still capture the main features of the pore size distribution of the samples, including the general shape and relative peak positions of the distribution curves. It is also found from the numerical example that the surface relaxivity of the sample could be extracted with the joint inversion of NMR and SIP data if the diffusion coefficient of the ions in the electrical double layer is known. Comparing to individual inversions, the joint inversion could improve the resolution of the estimated pore size distribution because of the addition of extra data sets. The proposed approach might constitute a first step towards a comprehensive joint inversion that can extract the full pore geometry information of a geomaterial from NMR and SIP data. Electrical properties, Magnetic properties, Microstructure, Hydrogeophysics, Joint inversion 1 INTRODUCTION The pore space plays a fundamental role in controlling fluid transport in porous media (Dullien 2012). For geomaterials such as soils and sedimentary rocks, the geometric characteristic of the pore space is usually described by porosity, specific surface area, pore size and its distribution. These parameters have been frequently used in estimating hydraulic properties of geomaterials. For instance, several characteristic pore sizes have been found to be correlated with the saturated permeability of sandstones (Johnson et al.1986; Katz & Thompson 1986); the pore size distribution is an important input for predicting water retention curves (Kosugi 1994; Perrier et al.1996) and unsaturated permeability of soils and rocks (Kosugi 1996; Tuli et al.2001; Vervoort & Cattle 2003). The geometric characteristics of pore space can be directly measured using some imaging methods, for example, X-ray imaging technique (see a review in Wildenschild & Sheppard 2013), scanning electron microscope (e.g. Timur et al.1971) and laser scanning confocal microscopy (e.g. Fredrich 1999). In addition, there are also many indirect methods that can be used to estimate the geometric characteristics of the pore space. One example is the mercury injection capillary porosimetry (MICP), and the measurement is usually interpreted in terms of pore throat size rather than pore size as demonstrated in Dullien (1981); another example is the gas adsorption analysis, which has been frequently used to infer the surface area and pore size of rock samples (e.g. Kuila & Prasad 2013). The use of geophysical techniques becomes increasingly attractive in estimating the geometric characteristics of the pore space in geomaterials (Tong et al.2006; Mohnke & Yaramanci 2008; Revil et al.2014), mainly due to its non-invasive nature and high spatiotemporal resolution. Another advantage of using geophysical tools is the applicability to field characterizations. The most promising geophysical tools in practice include nuclear magnetic resonance (NMR) and spectral induced polarization (SIP). The traditional NMR method measures the magnetic response of protons under specific sequences of applied magnetic fields (see Behroozmand et al.2015 for a review). In an NMR relaxation measurement, the signal magnitude is proportional to the absolute number of protons and thus the volume of hydrogen bearing fluid in the sample. The NMR decay curve can be decomposed to form a distribution of relaxation time (Borgia et al.1998), which is related to surface-to-volume ratio (Brownstein & Tarr 1979), thus the effective pore size (and its distribution) of the material can be estimated. The SIP method measures the electrical response of a material under an external alternating electric field, and the measurements can be expressed in terms of complex conductivity (see Kemna et al.2012 for a review). At low frequencies, the variation in complex conductivity of a geomaterial is mainly controlled by the polarization of the electrical double layer (EDL) forming at the solid–liquid interface (Lesmes & Morgan 2001). The measured complex conductivity spectra can be decomposed to obtain a distribution of relaxation time (Florsch et al.2014). Since the relaxation time is related to the geometry of the EDL (Schwarz 1962), the obtained relaxation time distribution thus can give an estimation for the pore size distribution of the material. For simplicity, in this study we use ‘NMR porosimetry’ and ‘SIP porosimetry’ to term the methods of using NMR and SIP data to estimate the pore size distribution of porous media, respectively. The geometric control of pore space on NMR and SIP responses of geomaterials is quite complicated, and there are still many ongoing studies towards a full understanding of the petrophysical relations (see some recent researches in Keating & Falzone 2013; Müller-Petke et al.2015; Niu et al.2016; Volkmann & Klitzsch 2016). In practice, some simplifications are usually made such that the pore size information of the material can be inferred from the NMR and SIP measurements. For instance, the NMR porosimetry assumes (1) that the magnetic relaxation of protons in porous media is in fast diffusion regime and (2) that the effects of bulk fluid diffusion and pore coupling on NMR relaxation are negligible (see a review in Behroozmand et al.2015). However, it has been found that these assumptions may be violated for materials with large pores and/or large values of surface relaxivity (Grunewald & Knight 2009; Dlubac et al.2014). In SIP porosimetry, it is generally assumed that the complex conductivity of geomaterials at frequencies lower than 100 kHz is solely from the EDL polarization and dielectric polarization (Leroy et al.2008; Revil et al.2014). However, other polarization may also contribute to the measured complex conductivity, such as interfacial polarization (Hanai 1960) and membrane polarization (Bücker & Hördt 2013). These approximate petrophysical relations used in NMR and SIP porosimetries may induce significant errors to the estimated pore size distribution. Inversion is needed to decompose the NMR relaxation curve and SIP spectra to obtain the relaxation time distribution. As is well known, most geophysical inversions are underdetermined due to the limited number of measurements; moreover, the inversion is also ill-posed because the measurements are usually contaminated by noise (Tarantola & Valette 1982). To mitigate these problems, joint inversion has become a common approach (e.g. Moorkamp et al.2011) because extra information from other data set could reduce the ambiguity or non-uniqueness of the inversion (Haber & Oldenburg 1997; Linde & Doetsch 2016). This has inspired us to consider the joint use of NMR and SIP data in estimating the pore size distribution of geomaterials. One similar attempt in this research direction is the work of Osterman et al. (2016) where both NMR and SIP measurements were used to estimate the pore geometric parameters and thus the permeability of sandstones. However, in their study different types of data sets are still processed separately, and a joint inversion of NMR and SIP data is still not realized yet. In this study, we explore the idea of joint inversion of NMR and SIP data for the estimation of pore size distributions of geomaterials. We first briefly review theoretical backgrounds of NMR and SIP porosimetries with emphasis on commonly applied petrophysical models. We then introduce a joint inversion approach that can simultaneously process NMR and SIP data. In the joint inversion, one aspect of significant importance is the trade-off between different data sets. We propose a two-step method to balance the weight between NMR and SIP measurements. A numerical experiment is designed to demonstrate the applicability of the proposed inversion method. In the last part of the study, we apply this joint inversion to a Berea sandstone sample and discuss the results as well as the limitations of our proposed approach. 2 THEORETICAL BACKGROUND 2.1 NMR relaxation in porous media In water-saturated geomaterials, NMR relaxation phenomena arises as the result of the interaction of diffusing nuclear spins of protons with their surrounding environment. The NMR relaxation can produce a bulk magnetization, which is governed by the phenomenological Bloch–Torrey equations. Experimentally, the transverse component of the magnetization can be detected using a series of oscillating magnetic field pulses [e.g. the Carr–Purcell–Meiboom–Gill (CPMG) pulse sequence] and the measured signal is in the form of a voltage decay curve, Exy(t). The solution to the diffusion problem can be expressed as the sum of a number of exponential decays (Brownstein & Tarr 1979) and thus the NMR signal can be expressed as   $${E_{xy}}\ \left( t \right) = {E_0}\ \mathop \sum \nolimits_i^{} {f_i}{e^{ - \frac{t}{{T_2^i}}}}$$ (1)where Exy is the transverse component of the NMR signal, E0 is the initial signal magnitude, fi is the relative intensities of the magnetic field relaxing with the relaxation time of $$T_2^i$$. The initial magnitude E0 is directly proportional to the number of protons and thus the amount of water in the material and therefore it can be used to estimate the water content of the material. Assuming a uniform magnetic field, the relaxation with time constant T2i only involves two processes (Brownstein & Tarr 1979) with   $$\frac{1}{{T_2^i}} = \frac{1}{{{T_{2B}}}}\ + \frac{1}{{T_{2S}^i}}$$ (2)where relaxation time T2B is related to the bulk fluid diffusion and $$T_{2S}^i$$ is related to surface relaxation. The diffusion rate of protons in bulk water in geomaterials is small, ∼0.3 s−1; the relaxation can be greatly enhanced at the mineral surface due to the presence of paramagnetic impurities such as Mn2+ and Fe3+ (Foley et al.1996). Compared with bulk diffusion, the surface relaxation could dominate the relaxation process in a saturated geomaterial such that $$\frac{1}{{{T_{2B}}}}$$ in eq. (2) could be dropped. Assuming no coupling between pores (e.g. see Grunewald & Knight 2009), a relaxation with time constant $$T_2^i$$ in eq. (1) can be related to a pore by   $$\frac{1}{{T_2^i}} = \frac{1}{{T_{2S}^i}}\ = \frac{1}{{\frac{{{r_i}}}{{\alpha {\rho _s}}} + \frac{{r_i^2}}{{2\alpha D}}}}\$$ (3)where ri is the effective size of the pore, α describes the shape of the pore (α = 3 for spheres), and D is the diffusion coefficient of water. The parameter ρs, termed as surface relaxivity, qualifies the efficiency of the surface relaxation. In practice, the term $$\frac{{{r^2_i}}}{{2\alpha D}}$$ in eq. (3) is usually ignored (i.e. assuming in fast diffusion regime), and the relaxation time can be related to the pore size by   $$T_2^i = \frac{{{r_i}}}{{\alpha {\rho _s}}}\ .$$ (4) Thus, the relative intensities fi of the relaxation time $$T_2^i$$ can be interpreted as pore size distribution. It should be noted that the bulk diffusion term in eq. (2) could be significant if the material is characterized with large pore space and low magnetic susceptibility (e.g. see Grunewald & Knight 2011). In that case, eq. (4) might be insufficient in describing the T2 relaxation time, and thus the bulk diffusion rate should also be included. 2.2 SIP and EDL polarization in porous media The SIP measures the complex conductivity response σ* of a porous medium at frequencies lower than 100 kHz. While the real part of the complex conductivity σ΄ quantifies the material's ability to conduct an electric current, the imaginary part σ″ reflects the storage capability of electrical energy from the applied electric field. The two parameters are not independent; they are related to each other by the Kramers–Kronig relation (e.g. see Volkmann & Klitzsch 2015). In this study, our discussion will be focused on the imaginary conductivity. The low-frequency imaginary conductivity of geomaterials is usually interpreted in the framework of EDL polarization. In geomaterials, the mineral surface is usually negatively charged because of isomorphous substitutions (McBride 1989) and/or the presence of hydroxyl group (Zhuravlev 2000). When the material is saturated with an electrolyte (e.g. NaCl solution), some ions (mostly cations) will accumulate at the liquid-solid interface to balance the surface charge, forming the EDL. Consider a spherical pore space saturated with NaCl solution. After applying an electric field, the ions in the EDL will redistribute to gain electric dipole moment. For the spherical EDL, the induced complex conductance follows Debye relaxation (Schwarz 1962),   $${{C}}_i^* = \frac{{i\omega \tau }}{{1 + i\omega \tau }}\ {{\rm{\Sigma }}_S},$$ (5)where ΣS is the surface conductance related to Stern layer. The relaxation time constant τ is related to the radius of the curvature of the EDL rc (see Appendix  A for details) by   $$\tau \ = \frac{{{r_c}^2}}{{2{D^ + }}}\$$ (6)where D+ is the diffusion coefficient of the ions in the EDL. Note that eq. (5) does not consider the ion exchange between EDL and bulk solution. For colloids, the EDL is coating the particles, and therefore, rc is equal the particle radius rd as shown in Schwarz (1962). For a spherical pore embedded in rock matrix, the radius of the curvature of the EDL rc is equal to the pore radius r, and therefore the relaxation time τ is   $$\tau \ = \frac{{{r^2}}}{{2{D^ + }}}\ .$$ (7) In geomaterials, the shape of the EDL seems to be controlled by the pore space and thus eq. (7) applies. The correlation between SIP relaxation time and pore size have been experimentally observed for a number of geological materials (Binley et al.2005; Kruschwitz et al.2010; Titov et al.2010; Revil et al.2012). For porous media with different pore sizes, the contribution from each pore can be summed up to form the overall response (e.g. see Leroy et al.2008),   $${{C}}_{{\rm{EDL}}}^* = \mathop \sum \nolimits_i^{} {f_i}\frac{{i\omega {\tau _i}}}{{1 + i\omega {\tau _i}}}\ {{\rm{\Sigma }}_S},$$ (8)where $${{C}}_{{\rm{EDL}}}^*$$ is the complex conductance induced by EDL polarization, and fi is the relative weighting of the pore with size ri, which corresponds to the relaxation time τi. For geological materials, the (complex) conductance induced by EDL (in S) is usually treated as a perturbation of the pore fluid conductivity (in S m−1) by considering Joule dissipation (Johnson et al.1986; Niu et al.2016). The increased complex conductivity of the pore fluid is sometimes termed as specific surface (complex) conductivity $${\rm{\sigma }}_{ss}^*$$, and it has the following form,   $${\rm{\sigma }}_{ss}^* = \frac{{2\left( {{{C}}_{{\rm{EDL}}}^* + C_S^0} \right)}}{{\rm{\Lambda }}}\ ,$$ (9)where $$C_S^0$$ is the DC surface conductance (e.g. see Revil & Glover 1998) and Λ is a characteristic pore size (Johnson et al.1986), which can be related to surface area-to-pore volume Spor and formation factor F by   $$\frac{2}{{\rm{\Lambda }}} = m{S_{{\rm{por}}}}$$ (10)where m is the cementation factor. After applying Archie's law, the complex conductivity of the material can be obtained as   $${{\rm{\sigma }}^*} = \frac{1}{F}\ \left( {{{\rm{\sigma }}_f} + {\rm{\sigma }}_{ss}^*} \right)$$ (11)where σf is the fluid conductivity. Inserting eqs (9) and (10) into eq. (11) yields   $${{\rm{\sigma }}^*} = \frac{{{{\rm{\sigma }}_f}}}{F}\ + \frac{{m{S_{{\rm{por}}}}}}{F}\left( {\mathop \sum \nolimits_i^{} {f_i}\frac{{i\omega {\tau _i}}}{{1 + i\omega {\tau _i}}}{{\rm{\Sigma }}_S} + C_S^0} \right).$$ (12) In eq. (12), the first term in the right-hand side (a real number) is associated with ionic conduction occurring in the bulk fluid in the pore space. The real part of the second term related to the surface conduction occurring in the EDL (Revil & Glover 1998) and the imaginary part quantifies the EDL polarization. The parameters ΣS and and $$C_S^0$$ can be expressed in terms of fluid chemistry and electrochemical properties of the EDL (e.g. see Revil 2012). In this study, we will not address the physical origin of the petrophysical parameters in eq. (12) such as F, m, ΣS and $$C_S^0$$; instead, we simplify eq. (12) by using several widely used lumped parameters. If the permittivity induced by dielectric polarization εinf is considered, the above equation is reformulated as,   $${{\rm{\sigma }}^*} = {\sigma _\infty }\ - \left( {{\sigma _\infty } - {\sigma _0}} \right)\mathop \sum \nolimits_i^{} \frac{{{f_i}}}{{1 + i\omega {\tau _i}}} + i\omega {\varepsilon _{{\rm{inf}}}},$$ (13)where σ∞ and σ0 are the real conductivities of the sample at low- and high-frequency limits, respectively and their relations to fluid chemistry and EDL properties can be found, for example, in Revil (2012). Since both σ∞ and σ0 are optimized during the inversion, to a certain extent, the effects of surface charge density, fluid salinity, and pH are also considered in interpreting the complex conductivity data. 2.3 Petrophysical simplifications in NMR and SIP porosimetries In this subsection, we summarize the petrophysical assumptions made in NMR and SIP porosimetries. In NMR porosimetry, it is assumed that the effect of diffusion on NMR relaxation in a magnetic gradient is negligible . The magnetic gradient arises from the magnetic susceptibility contrast between grain minerals and pore fluid (Cotts et al.1989). This assumption is usually valid if the NMR T2 relaxation is measured with CPMG pulse sequence with short echo times (Kleinberg & Horsfield 1990). Another assumption in NMR porosimetry is that the pores are isolated and there is no pore coupling effect (D’Orazio et al.1989). For most consolidated sandstones, the narrow pore throat could restrict the molecular diffusion between two interconnected pore spaces during the time of an NMR experiment, and therefore, this assumption is approximately valid (Borgia et al.1996). However, for samples with well-connected heterogeneous pores such as soil aggregates and analogous synthetic silica gels, the assumption might break down (Grunewald & Knight 2009). In addition, the NMR porosimetry assumes that T2 relaxation is controlled by surface relaxation in the fast-diffusion regime, in which the diffusion-controlled exchange between surface and bulk water is fast enough such that all the water molecules interact with the surface during the NMR measurement (Bryar et al.2000). This fast diffusion regime assumption may not be appropriate for materials with high surface relaxivity and/or large pores (e.g. Keating & Falzone 2013). In SIP porosimetry, the measurements are usually interpreted with the EDL polarization model. The most important assumption is that the pore size r is equal to the radius of EDL curvature rc so eq. (7) holds. For isolated pores, this assumption is clearly valid. For geomaterials with interconnected pores, although the correlation between r and τ has been experimentally observed (Revil et al.2012), it remains still unclear whether the relation is quadratic (Scott & Barker 2003; Titov et al.2010). Second, eq. (7) also assumes there is no ion exchange between EDL and bulk fluid (Lyklema et al.1983; Chelidze & Gueguen 1999; Niu & Revil 2016). As a matter of fact, the ion exchange could speed up the ion redistribution in the EDL and thus decrease the relaxation time τ as shown in Lyklema et al. (1983). In SIP porosimetry, the interfacial polarization and dielectric effects may also be considered. However, the theory combing EDL polarization, interfacial polarization, and dielectric polarization still significantly underestimates the measured permittivity of sandstones (Lesmes & Morgan 2001; Kremer et al.2016). Lastly, although not stated explicitly, SIP porosimetry applies the isolated pore model, which assumes the ion movement is limited to a single pore or EDL. That is, no ion exchange occurs between two adjacent pores during a SIP measurement. 3 JOINT INVERSION APPROACH As discussed, NMR and SIP porosimetries involve a number of assumptions about the petrophysical relations of a geomaterial. Since it is difficult to justify these assumptions for a real geomaterial, the pore size information obtained separately from NMR and SIP data may not be accurate. To reduce the uncertainties and thus increase the resolution, in this study we suggest to jointly process the NMR and SIP measurements for the estimation of the pore size distribution of geomaterials. The joint inversion scheme is illustrated in Fig. 1, and each step of the joint inversion is discussed in detail below. Figure 1. View largeDownload slide The joint inversion scheme to estimate pore size distribution from NMR and SIP measurements. Figure 1. View largeDownload slide The joint inversion scheme to estimate pore size distribution from NMR and SIP measurements. 3.1 Petrophysical relations In practice, there are two types of approaches for joint inversion of two geophysical measurements: the petrophysical approach (e.g. Jardani & Revil 2009) and structural approach (e.g. Gallardo & Meju 2003). The petrophysical approach assumes that different geophysical measurements are functions of the same properties of a material. Thus, multiple geophysical data can be integrated in the inversion. For instance, it assumes both resistivity and seismic velocity of a rock are controlled by its porosity and saturation, and thus resistivity and seismic data can be jointly inverted (Gao et al.2012). Alternatively, the structural approach assumes different geophysical methods yield structurally similar geophysical images of the subsurface. That is, the spatial and directional variations of different geophysical responses are similar. Mathematically, the structural difference between two images can be quantified using the cross-gradient function, and it can be combined with data misfit and smoothness constraints to form the objective function in the joint inversion (Gallardo & Meju 2003; Gallardo & Meju 2004). In this study, we apply the petrophysical approach in the joint inversion of NMR and SIP measurements. We assume both NMR T2 relaxation time and SIP relaxation time τ are related to the pore size r, that is, eqs (4) and (7) are valid. This implies that, for samples with similar surface properties (e.g. D+ and ρs), the relaxation time τ is proportional to squared T2. We verify the relation between τ and T2 using Sherwood sandstone samples from Osterman et al. (2016). The T2 squared and τ of these samples are cross-plotted in Fig. 2 where T2 is the peak NMR relaxation time and τ is determined from the characteristic frequency chosen using the approach in Revil et al. (2015). Despite an outliner point, a strong linear correlation is found between T22 and τ. The experimental evidence lends additional support to the analytical relation between T2 and τ as presented in the previous section. Figure 2. View largeDownload slide The relation between squared NMR T2 relaxation time and SIP relaxation time τ of Sherwood sandstone samples. The data are from Osterman et al. (2016). The datum with empty triangle is an outlier point. Figure 2. View largeDownload slide The relation between squared NMR T2 relaxation time and SIP relaxation time τ of Sherwood sandstone samples. The data are from Osterman et al. (2016). The datum with empty triangle is an outlier point. 3.2 Model parameters The common relation to pore size r makes it possible to jointly invert the pore size distribution [e.g. the relative intensities/weighting fi in eqs (1) and (13)] from NMR and SIP measurements. In addition to fi, other parameters related to the petrophysical models may also be inverted during the inversion, including σ∞, σ0, εinf, and E0 (e.g. Florsch et al.2014). The surface properties D+ and ρs are also involved in the inversion, but they cannot be recovered in separate inversions because of the non-uniqueness in $$\frac{r}{{3{\rho _s}}}$$ and $$\frac{{{r^2}}}{{2{D^ + }}}$$, as nicely demonstrated by Mohnke (2014) for ρs. As shown in the previous section and Fig. 2, relaxation time T2 and τ are related to each other for a given material through the pore size. This implies that the surface properties D+ and ρs are not independent, and therefore it is possible to recover one parameter once the other is known. In fact, previous research has shown, with extra data sets, the joint inversion can properly estimate the material's ρs (Mohnke 2014). In this study, we assume D+ is known as suggested by Revil (2013). The surface relaxivity ρs in the inversion can be either assumed as a constant (if we have a prior information on it) or updated iteratively. We define a model vector m, which contains all the unknowns to be recovered from the joint inversion. The model vector m can be expressed as   $${{\bf m}}\ = \left( {\begin{array}{@{}*{1}{c}@{}} {{\bf x}}\\ {\ln {\sigma _0}}\\ {\begin{array}{@{}*{1}{c}@{}} {\begin{array}{@{}*{1}{c}@{}} {\ln {\sigma _\infty }}\\ {\ln {\varepsilon _{{\rm{inf}}}}} \end{array}}\\ {\ln {E_0}}\\ {\ln {\rho _s}} \end{array}} \end{array}} \right)\ ,$$ (14)where x is a vector containing the relative weighting fi for the considered pore size ri in the form of xi = ln fi. The use of log-transformed pore size distribution (relative weighting) and other parameters in eq. (14) is to ensure positivity of the unknowns. Moreover, since the relative weighting fi of some pores are considerably smaller than the others, the log-transformation of fi ensures their relative changes are also important during the inversion (Weigand & Kemna 2016). 3.3 Objective function The objective function Ψ consists of data misfit Ψd and a model regularization term Ψm (Aster et al.2013),   $${\rm{\Psi \ }} = {{\rm{\Psi }}_d}\ + \lambda {{\rm{\Psi }}_m},$$ (15)where the damping parameter λ controls the relative contribution of the regularization term Ψm on the objective function. The regularization is needed to stabilize the inversion process to obtain a unique model (Titterington 1985). In NMR and SIP inversions, Tikhonov regularization is commonly used (Florsch et al.2014; Weigand & Kemna 2016). In this study, we also employ the Tikhonov regularization and the l2 regularization term Ψm can be expressed as   $${{\rm{\Psi }}_m} = \left\| {{{{\bf W}}_{{\bf m}}}\left( {{{\bf m}} - {{{\bf m}}_{{\rm{ref}}}}} \right)} \right\|_2^2$$ (16)where ${{{\bf W}}_{{\bf m}}} = ( \scriptsize{\begin{array}{@{}*{2}{c}@{}} {{\bf L}}&0\\ 0&{{\bf I}} \end{array}} ){\boldsymbol{\ }}$ is a regularization matrix and m is a reference model that may contain a priori information. While the discrete operator L in Wm applies smoothness constraints to the pore size distribution x in the model vector m, the identify matrix I forces the petrophysical unknowns (i.e. σ∞, σ0, εinf, E0 and ρs) to be close to reference values defined in mref (a similar treatment can be found in Florsch et al.2014). In this study, the operator L is chosen as the second-order row difference of an identity matrix (i.e. eq. 19 in Weigand & Kemna 2016). The minimization of the objective function Ψ can be realized using an iterative Gauss–Newton scheme such that the model update can be determined (e.g. Park & Van 1991) from   \begin{eqnarray} \left[ {{{{\bf J}}^T}{{\bf W}}_d^T{{{\bf W}}_d}{{\bf J}} + {\rm{\lambda }}{{\bf W}}_m^T{{{\bf W}}_{{m}}}} \right]\Delta {{{\bf m}}_k}& =& {{{\bf J}}^T}\ {{\bf W}}_d^T{{{\bf W}}_d}\left( {{{\bf d}} - {{\bf G}}\left( {{{{\bf m}}_k}} \right)} \right)\nonumber\\ && - {\rm{\lambda }}{{\bf W}}_m^T{{{\bf W}}_{{m}}}\left( {{{{\bf m}}_k} - {{{\bf m}}_{{\rm{ref}}}}} \right) \end{eqnarray} (17)where J is the sensitivity matrix corresponding to the forward modelling function G, which includes both eqs (1) and (13), Wd is the data weighting matrix, Δmk is the model update for the model mk at the iteration step k, and d is the data vector. In the joint inversion, the vector d is expressed as   $${{\bf d}}\ = \left( {\begin{array}{@{}*{1}{c}@{}} {{{\bf d}}_{{\rm{SIP}}}^{{\rm{Re}}}}\\ {{{\bf d}}_{{\rm{SIP}}}^{{\rm{Im}}}}\\ {{{{\bf d}}_{{\rm{NMR}}}}} \end{array}} \right)\$$ (18)where vectors $${{\bf d}}_{{\rm{SIP}}}^{{\rm{Re}}}$$ and $${{\bf d}}_{{\rm{SIP}}}^{{\rm{Im}}}$$ (with the size of m × 1) contain the real part and imaginary part of the frequency dependent complex conductivity, respectively, and the vector dNMR contains the measured NMR T2 relaxation time measurements. Since data points in an NMR experiment is usually at the order of thousands, it will be reduced in the inversion so that the dimension of dNMR (with the size of n × 1) is comparable to $${{\bf d}}_{{\rm{SIP}}}^{{\rm{Re}}}$$ and $${{\bf d}}_{{\rm{SIP}}}^{{\rm{Im}}}$$. The data weighting Wd matrix contains the errors related to the data vector, and it can be expressed as a diagonal matrix,   $${{{\bf W}}_d} = \left( {\begin{array}{c@{\quad}c@{\quad}c} {\beta {\boldsymbol \epsilon} _{{\rm{SIP}}}^{{\rm{Re}}}}&0&0\\ 0&{\beta {\boldsymbol \epsilon} _{{\rm{SIP}}}^{{\rm{Im}}}}&0\\ 0&0&{{{\boldsymbol \epsilon} _{{\rm{NMR}}}}} \end{array}} \right),$$ (19)where $${\boldsymbol \epsilon} _{{\rm{SIP}}}^{{\rm{Re}}} = \ {\rm{diag}}( {\frac{1}{{{\epsilon} _j^{Re}}}} )$$, $${\boldsymbol \epsilon} _{{\rm{SIP}}}^{{\rm{Im}}} = \ {\rm{diag}}( {\frac{1}{{\epsilon _j^{Im}}}} )$$, $${\epsilon _{{\rm{NMR}}}} = \ {\rm{diag}}( {\frac{1}{{{\epsilon _k}}}} )$$, and β is a trade-off parameter balancing the weight between SIP and NMR measurements. The parameters $$\epsilon _j^{\rm Re}$$ and $$\epsilon _j^{\rm Im}$$ are the errors related to the real part and imaginary part of the jth complex conductivity measurements, and εk is the error related to the kth reduced NMR measurement. In practice, these errors can be estimated from the standard deviation of the measurement noise (e.g. Pidlisecky et al.2007). Starting with a reference model mref, the model vector can be updated iteratively until the root-mean-square value of the misfit is smaller than a pre-defined value. At the step k, the model update can be calculated from eq. (17), and thus updated model mk+1 can be determined as   $${{{\bf m}}_{k + 1}} = {{{\bf m}}_k}\ + \chi \Delta {{{\bf m}}_k}$$ (20)where the step parameter χ belongs to [0, 1]. The restriction of χ ≤ 1 prevents the model from overshooting due to non-linearity, and it can be determined using a line search. For details, one can refer to Weigand & Kemna (2016). 3.4 Sensitivity calculation The sensitivity matrix (also known as Jacobian matrix) J in eq. (17) contains the partial derivatives of the model responses with respect to the model parameters, expressed as   \begin{eqnarray} {{\bf J}} \!=\!\! \left( {\begin{array}{@{}*{8}{c}@{}} {\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {f_1}}}}& \cdots &{\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {f_L}}}}&{\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {\sigma _0}}}}&{\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {\sigma _\infty }}}}&{\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {\varepsilon _{{\rm{inf}}}}}}}&{\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {E_0}}}}&{\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {\rho _s}}}}\\ \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {f_1}}}}& \cdots &{\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {f_L}}}}&{\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {\sigma _0}}}}&{\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {\sigma _\infty }}}}&{\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {\varepsilon _{{\rm{inf}}}}}}}&{\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {E_0}}}}&{\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {\rho _s}}}}\\ {\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {f_1}}}}& \cdots &{\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {f_L}}}}&{\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {\sigma _0}}}}&{\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {\sigma _\infty }}}}&{\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {\varepsilon _{{\rm{inf}}}}}}}&{\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {E_0}}}}&{\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {\rho _s}}}}\\ \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {f_1}}}}& \cdots &{\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {f_L}}}}&{\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {\sigma _0}}}}&{\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {\sigma _\infty }}}}&{\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {\varepsilon _{{\rm{inf}}}}}}}&{\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {E_0}}}}&{\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {\rho _s}}}}\\ {\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {f_1}}}}& \cdots &{\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {f_L}}}}&{\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {\sigma _0}}}}&{\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {\sigma _\infty }}}}&{\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {\varepsilon _{{\rm{inf}}}}}}}&{\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {E_0}}}}&{\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {\rho _s}}}}\\ \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {f_1}}}}& \cdots &{\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {f_L}}}}&{\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {\sigma _0}}}}&{\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {\sigma _\infty }}}}&{\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {\varepsilon _{{\rm{inf}}}}}}}&{\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {E_0}}}}&{\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {\rho _s}}}} \end{array}} \right),\nonumber\\ \end{eqnarray} (21)where σ'(ωi) and σ''(ωi) are the real and imaginary conductivities of the ith calculated SIP data (i = 1, 2, … m) and E(tj) is the jth calculated NMR data (j = 1, 2, … n). The sensitivity can be calculated by several methods (McGillivray & Oldenburg 1990). In this study, we analytically determine each element of the sensitivity matrix J in eq. (21). For most of the elements in Jacobian matrix J, they are dependent on the model parameters x. To ensure a correct update direction, in the inversion the matrix J is calculated at each iteration step. If one of the petrophysical model parameters is known (e.g. the surface relaxivity ρs is independently measured or estimated), we can set the related sensitivity components (e.g. last column in J in eq. 21) as zero so that the parameter will not be updated in the inversion. 3.5 Choice of λ and β The regularization parameter λ is a trade-off between the data misfit and model roughness. For separate inversion, the L-curve method (Hansen 2000) has been frequently used to determine λ. In the joint inversion, the data weighting parameter β balances the relative weighting between SIP and NMR data sets. Similar parameters can also be found in other joint inversions (e.g. Gao et al.2012). When more controlling parameters are involved in the inversion, it is always better to search these two parameters simultaneously using the generalized L-curve method as proposed in Murat et al. (2002). However, this method is computationally costly, and in practice the selection of multiple controlling parameters is mainly based on trial-and-error (e.g. Gallardo & Meju 2004) or simple criteria (e.g. Commer & Newman 2009). In this study, we suggest a two-step method for the determination of the regularization parameter λ and data weighting parameter β. The method is illustrated in Fig. 3 using the synthetic example presented in Section 4. In this two-step method, the regularization parameter λ is first fixed, for example, as an intermediate value. Then, the joint inversion is carried out with β varying over a broad range, for example, from 10−3 to 103 in Fig. 3(a). By plotting the data misfit and β, a reasonable trade-off can then be chosen based on the shape of the curve. For instance, in Fig. 3(a) the optimal β is chosen to ensure the data misfit reaches a minimum. The second step is to determine the regularization parameter λ. Similarly, the joint inversion is conducted with the optimal β and a varying λ. The parameter λ and the related data misfit can be plotted to form a curve where the minimum can be chosen, for example, as shown in Fig. 3(b). We chose the associated λ as the optimal regularization parameter. Other methods such as the conventional L-curve method and its updated version (Hansen et al.2007) can also be used here to determine the optimal regularization parameter. Figure 3. View largeDownload slide The determination of the data weighting parameter β and regularization parameter λ in the synthetic example presented in Section 4. The circle symbols indicate the positions of the optimal β and λ. Figure 3. View largeDownload slide The determination of the data weighting parameter β and regularization parameter λ in the synthetic example presented in Section 4. The circle symbols indicate the positions of the optimal β and λ. 4 A SYNTHETIC SAMPLE In this section, a synthetic porous sample with a bimodal pore size distribution (see Fig. 4a) is used to demonstrate the applicability of the proposed joint inversion approach. The density of the pore size is modelled by a combination of two lognormal distributions. The related mean pore sizes are 20 μm and 2 μm; the standard deviations are 0.3 and 0.3; and the relative weighting is 0.4 and 0.6. Numerical experiments are conducted to produce the NMR responses [using eqs (1) and (4)] and SIP responses [using eqs (7) and (13)] of the sample. The NMR and SIP measurements are then jointly inverted to estimate the pore size distribution. Two cases were tested and are presented here. In Case 1, the numerical experiments and the inversion use the same petrophysical models. In Case 2, it is assumed we lack information on the petrophysical models. Thus, the inversion uses simplified petrophysical models, which are different from those used in the numerical experiments. Figure 4. View largeDownload slide The synthetic example in Case 1: (a) the pore size distribution, (b) the ‘measured’ NMR T2 relaxation, (c) the ‘measured’ real conductivity spectra and (d) the ‘measured’ imaginary conductivity spectra. Figure 4. View largeDownload slide The synthetic example in Case 1: (a) the pore size distribution, (b) the ‘measured’ NMR T2 relaxation, (c) the ‘measured’ real conductivity spectra and (d) the ‘measured’ imaginary conductivity spectra. 4.1 Case 1: inversion with exact petrophysical models In Case 1, the NMR T2 relaxation was numerically calculated using eqs (1) and (4) and the ‘measured’ decay curve is shown in Fig. 4(b). In the NMR experiment, the initial signal E0 is set as 0.5 μV; the surface relaxivity ρs is set as 10 μm s−1; the NMR signal was recorded for 5 s and the total number of the measurements are 3000. The random Gaussian errors were added to the simulated decay curve and the simulated noise level is 1 per cent. The complex conductivity of the sample is numerically calculated using eqs (7) and (13), and the ‘measured’ real and imaginary conductivities are shown in Figs 4(c) and (d), respectively. In SIP experiment, the following petrophysical parameters were used: σ0 = 0.045 S m−1, σ∞ = 0.05 S m−1, D+ = 1.32 × 10−9 m2 s−1, εinf = 100ε0 (where ε0 is the vacuum dielectric permittivity, 8.85 × 10−12 F m−1). The use of 100ε0 as the high frequency permittivity εinf is consistent with experimental measurements (e.g. see Jougnot et al.2010; Kremer et al.2016). The origin of this abnormally high value (for comparison, the permittivity of pure water is ∼79ε0) is still unclear and it could be caused by Maxwell-Wagner polarization, grain shape influence, inductive electromagnetic coupling or electrode influence (Kremer et al.2016). The applied minimum and maximum frequencies are 10−3 Hz and 100 kHz, respectively, and the total number of SIP measurements is 50. Similarly, random Gaussian errors are added to the calculation to simulate the measurement noise. The noise level is 1 per cent for real conductivity and 35 per cent for imaginary conductivity. Three individual calculations were conducted for each measured frequency. The mean value was used as the measured conductivity and the standard deviation was used as the measurement error as shown in Figs 4(c) and (d). The proposed joint inversion approach was used to estimate the pore size distribution of the sample based on the NMR and SIP measurements in Fig. 4. The proposed two-step method was used to determine the regularization parameter λ and data weighting parameter β. As shown in Figs 3(a) and (b), the reasonable trade-off was chosen when the data misfit reaches a minimum, and the determined optimal λ and β are 9771 and 153, respectively. In the joint inversion, the surface relaxivity was taken as the true value, that is, 10 μm s−1; and therefore, it was not updated during the inversion. The joint-inversion result is shown in Fig. 5(a). Separate NMR and SIP inversions were also carried out with the same regularization parameter, and the results are shown in Figs 5(b) and (c). The true density (normalized value) of the pore size is also plotted in Fig. 5 as well as discrepancies between the inverted and true densities. Figure 5. View largeDownload slide The inversion results of Case 1: (a) joint inversion, (b) separate NMR inversion, (c) separate SIP inversion and (d) discrepancies between inverted and true distributions of the pore size. Figure 5. View largeDownload slide The inversion results of Case 1: (a) joint inversion, (b) separate NMR inversion, (c) separate SIP inversion and (d) discrepancies between inverted and true distributions of the pore size. As shown in Fig. 5(a), the joint inversion successfully recovers the bimodal shape of the pore size distribution. The estimated positions of the two peaks are very close to the true locations although the related densities are slightly smaller than the true values. In contrast, the results from separate inversions have very low resolutions. The NMR inversion, to a certain extent, recovers the bimodal nature of the pore size but it does not correctly locate the dominant pore size; the SIP inversion renders a pore size distribution with only one peak. Comparing with separate inversions, the joint inversion significantly improves the resolution of the estimated pore size distribution. Clearly, this improvement is due to the addition of extra measurements in the inversion. We also tested whether the joint inversion can correctly estimate the surface relaxivity of the sample by allowing ρs to update during the inversion. The initial value of ρs was set as 1 μm s−1and other parameters were unchanged. The inverted pore size distribution is very close to the one presented in Fig. 5(a). The evolution of ρs during the inversion is plotted in Fig. 6. It shows that although the initial ρs is relatively small, after several iterations, it approaches to the true value. The final inverted surface relaxivity is 9.47 μm s−1, which is very close to the true value 10 μm s−1. Figure 6. View largeDownload slide The surface relaxivity evolution during the joint inversion in Case 1. Figure 6. View largeDownload slide The surface relaxivity evolution during the joint inversion in Case 1. 4.2 Case 2: inversion with inaccurate petrophysical models In Case 2, we use eqs (1) and (3) to simulate the NMR response of the sample so that the NMR relaxation is not necessarily within the fast diffusion regime. In SIP experiments, a Cole–Cole model is used to describe the unknown polarization that might be responsible for the anomalously high permittivity at intermediate frequencies (e.g. Lesmes & Morgan 2001; Kremer et al.2016). Thus, the complex conductivity is calculated by combing eqs (7) and (13) and the added Cole–Cole model. The NMR and SIP responses calculated using these ‘realistic’ petrophysical models are shown in Fig. 7. In the calculation, same model parameters were used here except that the surface relaxivity is changed to 100 μm s−1. For comparison, NMR and SIP responses calculated using the simplified petrophysical models (i.e. those used in Case 1) were also presented. As shown in Fig. 7, the ‘realistic’ and simplified petrophysical models give different NMR and SIP responses. The former relaxes the fast diffusion condition, and thus the related NMR relaxation rate is apparently slower (Fig. 7a). Moreover, ‘realistic’ petrophysical models give relatively large imaginary conductivities at frequencies larger 100 Hz (Fig. 7b), and the spectrum resembles experimental measurements shown in Lesmes & Morgan (2001). Figure 7. View largeDownload slide The measured (a) NMR and (b) SIP responses of the synthetic sample in Case 2. Compared with simplified petrophysical models, the realistic ones relax the fast diffusion condition and consider the unknown polarization responsible for anomalously high imaginary conductivity at intermediate frequencies. Figure 7. View largeDownload slide The measured (a) NMR and (b) SIP responses of the synthetic sample in Case 2. Compared with simplified petrophysical models, the realistic ones relax the fast diffusion condition and consider the unknown polarization responsible for anomalously high imaginary conductivity at intermediate frequencies. The NMR and SIP responses (simulated from ‘realistic’ petrophysical models) are inverted to estimate the pore size distribution jointly or separately. In the inversions, it is assumed that we do not know the exact petrophysical relations and the simplified petrophysical models [eqs (1), (4), (7), and (13)] are used. The estimated pore size distributions from joint and separate inversions are shown in Fig. 8. The joint inversion still recovers the bimodal nature of the pore size and two peak locations are also well defined. However, the estimated densities are significantly lower than the true values. In addition, the inversion gives a very broad range for the pore size distribution. Compared to Fig. 5(a), the joint inversion results in Fig. 8(a) differ notably from the true model. This is attributed to the use of the simplified petrophysical models in the inversion. Again, we compare the results of joint inversion to those from separate inversions. It appears the separate inversion, neither SIP nor NMR inversion, can capture the bimodal nature of the pore size. The joint inversion still improves the quality of the estimated pore size distribution despite that we lack enough information on the exact petrophysical models of the sample. Figure 8. View largeDownload slide The inversion results of Case 2: (a) joint inversion, (b) separate NMR inversion, (c) separate SIP inversion and (d) discrepancies between inverted and true distributions of the pore size. Figure 8. View largeDownload slide The inversion results of Case 2: (a) joint inversion, (b) separate NMR inversion, (c) separate SIP inversion and (d) discrepancies between inverted and true distributions of the pore size. 5 APPLICATION TO BEREA SANDSTONE In this section, we apply the proposed joint inversion approach to a Berea sandstone. The estimated pore size distribution is compared with those from separate inversions. 5.1 Sample and measurements The well-characterized Berea sandstone is chosen in this study because both NMR and SIP data are available. The sample has a porosity ϕ of 20.2 per cent and the specific surface area Sm is determined with nitrogen adsorption as 0.5 m2 g−1. Previous research shows quartz is the most abundant mineral in the sandstone and it also contains a significant amount of clays (e.g. Knight & Nur 1987). The intrinsic formation factor F of the sample is 16.3 as determined from multi-salinity measurement. The sample was saturated with NaCl solution (σw = 0.043 S m−1) for NMR and SIP measurement. The NMR T2 relaxation of the sample was measured using a 2 MHz Magritek Rock Core AnalyzerTM with CPMG pulse sequence and the echo spacing is 100 μs. The complex conductivity of the sample was measured in the frequency range between 10−3 Hz and 10 kHz using the SIP-LAB-II instrument (Radic). The NMR and SIP responses of the sample are shown in Fig. 9. Figure 9. View largeDownload slide The measured and fitted (a) NMR and (b) SIP responses (imaginary conductivity) of the Berea sandstone sample. The fitted curves are from joint and separate inversions. Figure 9. View largeDownload slide The measured and fitted (a) NMR and (b) SIP responses (imaginary conductivity) of the Berea sandstone sample. The fitted curves are from joint and separate inversions. 5.2 Inversion results The NMR and SIP data of the sample in Fig. 9 are jointly inverted to estimate the pore size distribution of Berea sandstone. The data weighting and regularization parameters β and λ are determined using the proposed two-step method as shown in Fig. 10. In the inversion, the diffusion coefficient of ions in the EDL D+ was taken as 3.2 × 10−12 m2 s−1 as suggested by Revil (2012), considering that the sample has a significant amount of clays; the value may need to be adjusted by trial-and-error in cases where the calculated and measured peak frequencies in SIP spectra do not coincide. The surface relaxivity is allowed to change in the inversion and the initial value is set as ρs = 5 μm s−1. However, we find ρs does not change during the inversion. The reason could be that the initial value was too close to the true surface relaxivity of the sample. The recovered petrophysical parameters from joint inversion are as follows: σ∞ = 0.0028 S m−1, σ0 = 0.0026 S m−1, εinf = 128ε0, E0 = 7 μV, ρs = 4.995 μm s−1. For comparison, separate inversions are also conducted. In NMR inversion, the surface relaxivity ρs is kept as 5 μm s−1; in SIP inversion, the diffusion coefficient D+ is taken as 3.2 × 10−12 m2 s−1. The regularization parameters used in separate inversions are individually determined using the traditional L-curve method so they may be different from that used in the joint inversion. The estimated pore size distributions from separate and joint inversions are shown in Fig. 11(a). The NMR and SIP responses calculated from the inversion results are also plotted in Fig. 9 for comparison. Figure 10. View largeDownload slide The determination of the data weighting parameter β and regularization parameter λ for the Berea sandstone sample. The circle symbols indicate the positions of the optimal β and λ. Figure 10. View largeDownload slide The determination of the data weighting parameter β and regularization parameter λ for the Berea sandstone sample. The circle symbols indicate the positions of the optimal β and λ. Figure 11. View largeDownload slide The estimated pore size distribution of the Berea sandstone using (a) SIP and NMR porosimetries and (b) mercury injection capillary pressure (MICP). The MICP-determined pore (throat) size distribution is from the Berea sandstone sample in Shi et al. (2011), which has a similar porosity (ϕ = 0.19) with the one used in this study (ϕ = 0.2). Figure 11. View largeDownload slide The estimated pore size distribution of the Berea sandstone using (a) SIP and NMR porosimetries and (b) mercury injection capillary pressure (MICP). The MICP-determined pore (throat) size distribution is from the Berea sandstone sample in Shi et al. (2011), which has a similar porosity (ϕ = 0.19) with the one used in this study (ϕ = 0.2). The microstructure of Berea sandstone has been extensively studied using different methods, including MICP (Shi et al.2011) and X-ray CT scan (Dong & Blunt 2009). Here, we use the sample in Shi et al. (2011) as an example to show the main features of pore (throat) size of Berea sandstone. The sample has a porosity of ∼19 per cent, which is very close to the sample in our study (∼20 per cent). The pore (throat) size distribution measured from MICP (Shi et al.2011) is plotted in Fig. 11(b). One profound feature of the sample is the narrow range of the pore size, which centres on ∼10 μm. In contrast, the pore size distribution estimated from separate SIP inversion is quite broad, covering more than 2 orders. Comparing to MICP results, the SIP inversion significantly overestimates the contribution of small pores (r < 3 μm) although the imaginary conductivity spectrum is well reproduced base on this broad pore size distribution (Fig. 9b). As shown in Fig. 11, the pore size estimated from joint inversion is consistent with the MICP results; both the shape and the relative position of the curve are well recovered, demonstrating the applicability of the joint inversion approach. The NMR response calculated from joint inversion fits the measurement perfectly. The simulated imaginary conductivity spectra agree with the measurements only at low (<0.1 Hz) and high frequencies (>10 kHz). At intermediate frequencies, the simulated σ″ is much smaller than the measurements. The discrepancy may be explained by the presence of other polarization, which is responsible for the measured quadrature conductivity at intermediate frequencies (Kremer et al.2016). The estimated pore size distribution from NMR inversion is very similar to that of joint inversion, indicating that more weights are put on NMR data during the joint inversion. That implies, for this sample, the NMR data contain more information on the pore size than the SIP measurements. 5.3 Discussion In this subsection, we briefly discuss limitations of the proposed joint inversion approach and possible directions for further studies. First, we point out that the Stern layer polarization model (i.e. eq. 5) may only be valid for relatively large pores (>∼0.1 μm) because it assumes the EDL thickness is much smaller than the radius of the curvature (see Schwarz 1962 or Appendix  A in this study). Therefore, for geomaterials characterized by small pore size (<∼0.1 μm) such as mudstones or clayey soils, the Stern layer polarization model may not be able to describe the measured SIP spectrum. In that case, the joint inversion of SIP and NMR data with our approach may not improve the estimation of pore size distribution if compared with separate NMR porosimetry. In addition, since both NMR and SIP porosimetries assume isolated pores, the joint inversion approach may not work properly for geomaterials with well-connected heterogeneous pores such as loosely packed soils. Despite the fact that the Stern layer polarization model has been frequently used in practice for various applications (e.g. Osterman et al. 2016; Leroy et al. 2017; Revil et al. 2017), it is still unclear whether it alone is sufficient to explain the measured complex conductivity of geomaterials over a broad frequency range. Thus, it is necessary to assess the applicability of the Stern layer polarization model for a given sample, and the proposed joint inversion approach is well suited for this purpose. Indeed, the simulated SIP response in Fig. 9(b) based on joint inversion of NMR and SIP data has indicated that Stern layer polarization may not explain the high imaginary conductivity of Berea sandstone at ∼100 Hz. In addition to EDL polarization, membrane polarization (Marshall & Madden 1959) is also commonly used to explain the induced polarization of porous media (e.g. Bairlein et al. 2016; Hördt et al. 2016). Different from Stern layer polarization, which relates the relaxation time to pore/particle size, membrane polarization suggests that the relaxation time may also be controlled by the geometry of the pore throat (Bücker & Hördt 2013). This conclusion was also supported by experimental evidence (e.g. Scott & Barker 2003). Considering that pore throat is usually smaller than pore space in geomaterials, membrane polarization could potentially explain the abnormal high imaginary conductivity at intermediate frequencies reported in Lesmes & Morgan (2001) and Kremer et al. (2016). In this perspective, the joint inversion incorporating membrane polarization would give additional information on the pore geometry of a geomaterial. Therefore, our current work might constitute a first step towards a comprehensive joint inversion that can extract the full pore geometry information from NMR and SIP data. 6 CONCLUSIONS In this study, we propose to jointly invert NMR and SIP measurements for the estimation of the pore size distribution of geomaterials. The joint inversion is based on the fact that both NMR T2 relaxation and low-frequency induced polarization of geomaterials are influenced by geometric characteristics of the pore space. In the joint inversion, the NMR T2 relaxation time and the SIP relaxation time τ are linked through their petrophysical models. The Tikhonov regularization is used to constrain the smoothness of the pore size distribution. The trade-offs between two data sets and the balance between data misfit and model regularization are determined using a two-step method. The inversion problem is solved iteratively using the Gauss–Newton method and the related sensitivity matrix at each iteration is determined analytically. Numerical examples have shown that the joint inversion could recover the main features of the pore size of the synthetic sample even if petrophysical models used in the inversion are inaccurate. Comparing to separate inversions, the addition of extra data set can improve the resolution of the estimated pore size distribution. In addition, the results from synthetic example indicate that the surface relaxivity could be properly estimated by the joint inversion when the diffusion coefficient of ions in the EDL is known. The proposed method has been applied to a Berea sandstone sample and the estimated pore size distribution is fairly consistent with that estimated from mercury injection capillary pressure measurements, demonstrating the applicability and usefulness of the proposed method. Acknowledgements We thank the editor Mark Everett for handling our submission, and Norbert Klitzsch along with an anonymous reviewer for their valuable review comments. We thank Manika Prasad at the Department of Petroleum Engineering of Colorado School of Mines for providing the Berea sandstone data sets. REFERENCES Aster R.C., Borchers B., Thurber C.H., 2013. Tikhonov regularization, in Parameter Estimation and Inverse Problems , 2nd edn, pp. 93– 127, Academic Press. Google Scholar CrossRef Search ADS   Bairlein K., Bücker M., Hördt A., Hinze B. 2016. Temperature dependence of spectral induced polarization data: experimental results and membrane polarization theory, Geophys. J. Int. , 205( 1), 440– 453. https://doi.org/10.1093/gji/ggw027 Google Scholar CrossRef Search ADS   Behroozmand A.A., Keating K., Auken E., 2015. A review of the principles and applications of the NMR technique for near-surface characterization, Surv. Geophys. , 36( 1), 27– 85. https://doi.org/10.1007/s10712-014-9304-0 Google Scholar CrossRef Search ADS   Binley A., Slater L.D., Fukes M., Cassiani G., 2005. Relationship between spectral induced polarization and hydraulic properties of saturated and unsaturated sandstone, Water Resour. Res. , 41, W12417, doi:10.1029/2005WR004202. https://doi.org/10.1029/2005WR004202 Borgia G.C., Bortolotti V., Brancolini A., Brown R.J.S., Fantazzini P., 1996. Developments in core analysis by NMR measurements, Magn. Reson. Imaging , 14( 7–8), 751– 760. https://doi.org/10.1016/S0730-725X(96)00160-9 Google Scholar CrossRef Search ADS PubMed  Borgia G.C., Brown R.J.S., Fantazzini P., 1998. Uniform-penalty inversion of multiexponential decay data, J. Magn. Reson. , 132( 1), 65– 77. https://doi.org/10.1006/jmre.1998.1387 Google Scholar CrossRef Search ADS PubMed  Brownstein K.R., Tarr C.E., 1979. Importance of classical diffusion in NMR studies of water in biological cells, Phys. Rev. A , 19( 6), 2446– 2453. https://doi.org/10.1103/PhysRevA.19.2446 Google Scholar CrossRef Search ADS   Bryar T.R., Daughney C.J., Knight R.J., 2000. Paramagnetic effects of iron(III) species on nuclear magnetic relaxation of fluid protons in porous media, J. Magn. Reson. , 142( 1), 74– 85. https://doi.org/10.1006/jmre.1999.1917 Google Scholar CrossRef Search ADS PubMed  Bücker M., Hördt A., 2013. Analytical modelling of membrane polarization with explicit parametrization of pore radii and the electrical double layer, Geophys. J. Int. , 194( 2), 804– 813. https://doi.org/10.1093/gji/ggt136 Google Scholar CrossRef Search ADS   Chelidze T.L., Gueguen Y., 1999. Electrical spectroscopy of porous rocks: a review—I. Theoretical models, Geophys. J. Int. , 137( 1), 1– 15. https://doi.org/10.1046/j.1365-246x.1999.00799.x Google Scholar CrossRef Search ADS   Commer M., Newman G.A., 2009. Three-dimensional controlled-source electromagnetic and magnetotelluric joint inversion, Geophys. J. Int. , 178 ( 3), 1305– 1316. https://doi.org/10.1111/j.1365-246X.2009.04216.x Google Scholar CrossRef Search ADS   Cotts R.M., Hoch M.J.R., Sun T., Markert J.T., 1989. Pulsed field gradient stimulated echo methods for improved NMR diffusion measurements in heterogeneous systems, J. Magn. Reson. (1969) , 83( 2), 252– 266. https://doi.org/10.1016/0022-2364(89)90189-3 Google Scholar CrossRef Search ADS   Dlubac K., Knight R., Keating K., 2014. A numerical study of the relationship between NMR relaxation and permeability in sands and gravels, Near Surf. Geophys. , 12( 2), 219– 230. https://doi.org/10.3997/1873-0604.2013042 Dong H., Blunt M.J., 2009. Pore-network extraction from micro-computerized-tomography images, Phys. Rev. E , 80( 3), doi:10.1103/PhysRevE.80.036307. https://doi.org/10.1103/PhysRevE.80.036307 D’Orazio F., Tarczon J.C., Halperin W.P., Eguchi K., Mizusaki T., 1989. Application of nuclear magnetic resonance pore structure analysis to porous silica glass, J. Appl. Phys. , 65( 2), 742– 751. https://doi.org/10.1063/1.343088 Google Scholar CrossRef Search ADS   Dullien F.A., 2012. Porous Media: Fluid Transport and Pore Structure , Academic press. Dullien F.A.L., 1981. Wood's metal porosimetry and its relation to mercury porosimetry, Powder Technol. , 29( 1), 109– 116. https://doi.org/10.1016/0032-5910(81)85009-7 Google Scholar CrossRef Search ADS   Florsch N., Revil A., Camerlynck C., 2014. Inversion of generalized relaxation time distributions with optimized damping parameter, J. Appl. Geophys. , 109, 119– 132. https://doi.org/10.1016/j.jappgeo.2014.07.013 Google Scholar CrossRef Search ADS   Foley I., Farooqui S.A., Kleinberg R.L., 1996. Effect of paramagnetic ions on NMR relaxation of fluids at solid surfaces, J. Magn. Reson. A , 123( 1), 95– 104. https://doi.org/10.1006/jmra.1996.0218 Google Scholar CrossRef Search ADS PubMed  Fredrich J.T., 1999. 3D imaging of porous media using laser scanning confocal microscopy with application to microscale transport processes, Phys. Chem. Earth, Part A: Solid Earth Geod. , 24( 7), 551– 561. https://doi.org/10.1016/S1464-1895(99)00079-4 Google Scholar CrossRef Search ADS   Gallardo L.A., Meju M.A., 2003. Characterization of heterogeneous near-surface materials by joint 2D inversion of dc resistivity and seismic data, Geophys. Res. Lett. , 30( 13), 1658, doi:10.1029/2003GL017370. https://doi.org/10.1029/2003GL017370 Gallardo L.A., Meju M.A., 2004. Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints, J. geophys. Res. , 109, B03311, doi:10.1029/2003JB002716. https://doi.org/10.1029/2003JB002716 Gao G., Abubakar A., Habashy T.M., 2012. Joint petrophysical inversion of electromagnetic and full-waveform seismic data, Geophysics , 77( 3), WA3– WA18. https://doi.org/10.1190/geo2011-0157.1 Google Scholar CrossRef Search ADS   Grunewald E., Knight R., 2009. A laboratory study of NMR relaxation times and pore coupling in heterogeneous media, Geophysics , 74( 6), E215– E221. https://doi.org/10.1190/1.3223712 Google Scholar CrossRef Search ADS   Grunewald E., Knight R., 2011. The effect of pore size and magnetic susceptibility on the surface NMR relaxation parameter T2*, Near Surf. Geophys. , 9( 2), 169– 178. Haber E., Oldenburg D., 1997. Joint inversion: a structural approach, Inverse Probl. , 13( 1), 63– 77. https://doi.org/10.1088/0266-5611/13/1/006 Google Scholar CrossRef Search ADS   Hanai T., 1960. Theory of the dielectric dispersion due to the interfacial polarization and its application to emulsions, Kolloid-Zeitschrift , 171( 1), 23– 31. https://doi.org/10.1007/BF01520320 Google Scholar CrossRef Search ADS   Hansen P.C., 2000. The L-curve and its use in the numerical treatment of inverse problems, in Computational Inverse Problems in Electrocardiology , pp. 119– 142, ed. Johnston P., WIT Press. Hansen P.C., Jensen T.K., Rodriguez G., 2007. An adaptive pruning algorithm for the discrete L-curve criterion, J. Comput. Appl. Math. , 198( 2), 483– 492. https://doi.org/10.1016/j.cam.2005.09.026 Google Scholar CrossRef Search ADS   Hördt A., Bairlein K., Bielefeld A., Bücker M., Kuhn E., Nordsiek S., Stebner H. 2016. The dependence of induced polarization on fluid salinity and pH, studied with an extended model of membrane polarization, J. Appl. Geophys. , 135, 408– 417. https://doi.org/10.1016/j.jappgeo.2016.02.007 Google Scholar CrossRef Search ADS   Jardani A., Revil A., 2009. Stochastic joint inversion of temperature and self-potential data, Geophys. J. Int. , 179( 1), 640– 654. https://doi.org/10.1111/j.1365-246X.2009.04295.x Google Scholar CrossRef Search ADS   Johnson D.L., Koplik J., Schwartz L.M., 1986. New pore-size parameter characterizing transport in porous media, Phys. Rev. Lett. , 57( 20), 2564– 2567. https://doi.org/10.1103/PhysRevLett.57.2564 Google Scholar CrossRef Search ADS PubMed  Jougnot D., Ghorbani A., Revil A., Leroy P., Cosenza P. 2010. Spectral induced polarization of partially saturated clay-rocks: a mechanistic approach, Geophys. J. Int. , 180( 1), 210– 224. https://doi.org/10.1111/j.1365-246X.2009.04426.x Google Scholar CrossRef Search ADS   Katz A.J., Thompson A.H., 1986. Quantitative prediction of permeability in porous rock, Phys. Rev. B , 34( 11), 8179– 8181. https://doi.org/10.1103/PhysRevB.34.8179 Google Scholar CrossRef Search ADS   Keating K., Falzone S., 2013. Relating nuclear magnetic resonance relaxation time distributions to void-size distributions for unconsolidated sand packs, Geophysics , 78( 6), D461– D472. https://doi.org/10.1190/geo2012-0461.1 Google Scholar CrossRef Search ADS   Kemna A. et al.  , 2012. An overview of the spectral induced polarization method for near-surface applications, Near Surf. Geophys. , 10( 6), 453– 468. https://doi.org/10.3997/1873-0604.2012027 Kleinberg R.L., Horsfield M.A., 1990. Transverse relaxation processes in porous sedimentary rock, J. Magn. Reson. (1969) , 88( 1), 9– 19. https://doi.org/10.1016/0022-2364(90)90104-H Google Scholar CrossRef Search ADS   Knight R.J., Nur A., 1987. The dielectric constant of sandstones, 60 kHz to 4 MHz, Geophysics , 52( 5), 644– 654. https://doi.org/10.1190/1.1442332 Google Scholar CrossRef Search ADS   Kosugi K., 1996. Lognormal distribution model for unsaturated soil hydraulic properties, Water Resour. Res. , 32( 9), 2697– 2703. https://doi.org/10.1029/96WR01776 Google Scholar CrossRef Search ADS   Kosugi K.I., 1994. Three-parameter lognormal distribution model for soil water retention, Water Resour. Res. , 30( 4), 891– 901. https://doi.org/10.1029/93WR02931 Google Scholar CrossRef Search ADS   Kremer T., Schmutz M., Leroy P., Agrinier P., Maineult A., 2016. Modelling the spectral induced polarization response of water-saturated sands in the intermediate frequency range (102–105 Hz) using mechanistic and empirical approaches, Geophys. J. Int. , 207( 2), 1303– 1312. https://doi.org/10.1093/gji/ggw334 Google Scholar CrossRef Search ADS   Kruschwitz S., Binley A., Lesmes D., Elshenawy A., 2010. Textural controls on low-frequency electrical spectra of porous media, Geophysics , 75( 4), WA113– WA123. https://doi.org/10.1190/1.3479835 Google Scholar CrossRef Search ADS   Kuila U., Prasad M., 2013. Specific surface area and pore-size distribution in clays and shales, Geophys. Prospect. , 61( 2), 341– 362. https://doi.org/10.1111/1365-2478.12028 Google Scholar CrossRef Search ADS   Leroy P., Revil A., Kemna A., Cosenza P., Ghorbani A., 2008. Complex conductivity of water-saturated packs of glass beads, J. Colloid Interface Sci. , 321( 1), 103– 117. https://doi.org/10.1016/j.jcis.2007.12.031 Google Scholar CrossRef Search ADS PubMed  Leroy P., Li S., Jougnot D., Revil A., Wu Y., 2017. Modelling the evolution of complex conductivity during calcite precipitation on glass beads. Geophys. J. Int. , 209( 1), 123– 140. Lesmes D.P., Morgan F.D., 2001. Dielectric spectroscopy of sedimentary rocks, J. geophys. Res. , 106( B7), 13 329–13 346. https://doi.org/10.1029/2000JB900402 Google Scholar CrossRef Search ADS   Linde N., Doetsch J., 2016. Joint inversion in hydrogeophysics and near-surface geophysics, in Integrated Imaging of the Earth , ch. 7, pp. 119– 135, eds Moorkamp M., Lelievre P., Linde N., Khan A., John Wiley & Sons, Inc. Google Scholar CrossRef Search ADS   Lyklema J., Dukhin S., Shilov V., 1983. The relaxation of the double layer around colloidal particles and the low-frequency dielectric dispersion, J. Electroanal. Chem. Interfacial Electrochem. , 143( 1–2), 1– 21. https://doi.org/10.1016/S0022-0728(83)80251-4 Google Scholar CrossRef Search ADS   Marshall D.J., Madden T.R. 1959. Induced polarization, a study of its causes, Geophysics , 24( 4), 790– 816. https://doi.org/10.1190/1.1438659 Google Scholar CrossRef Search ADS   McBride M., 1989. Surface chemistry of soil minerals, in Minerals in Soil Environments , pp. 3588, eds Dixon, J.B. & Weed, S.B., Soil Science Society of America. McGillivray P.R., Oldenburg D.W., 1990. Methods for calculating Frechet derivatives and sensitivities for the non-linear inverse problem: a comparative study, Geophys. Prospect. , 38, 499– 524. https://doi.org/10.1111/j.1365-2478.1990.tb01859.x Google Scholar CrossRef Search ADS   Mohnke O., 2014. Jointly deriving NMR surface relaxivity and pore size distributions by NMR relaxation experiments on partially desaturated rocks, Water Resour. Res. , 50, 5309– 5321. https://doi.org/10.1002/2014WR015282 Google Scholar CrossRef Search ADS   Mohnke O., Yaramanci U., 2008. Pore size distributions and hydraulic conductivities of rocks derived from Magnetic Resonance Sounding relaxation data using multi-exponential decay time inversion, J. Appl. Geophys. , 66, 73– 81. https://doi.org/10.1016/j.jappgeo.2008.05.002 Google Scholar CrossRef Search ADS   Moorkamp M., Heincke B., Jegen M., Roberts A.W., Hobbs R.W., 2011. A framework for 3-D joint inversion of MT, gravity and seismic refraction data, Geophys. J. Int. , 184, 477– 493. https://doi.org/10.1111/j.1365-246X.2010.04856.x Google Scholar CrossRef Search ADS   Müller-Petke M., Dlugosch R., Lehmann-Horn J., Ronczka M., 2015. Nuclear magnetic resonance average pore-size estimations outside the fast-diffusion regime, Geophysics , 80, D195– D206. https://doi.org/10.1190/geo2014-0167.1 Google Scholar CrossRef Search ADS   Murat B., Misha E.K., Eric L.M., 2002. Efficient determination of multiple regularization parameters in a generalized L-curve framework, Inverse Probl. , 18, 1161– 1183. https://doi.org/10.1088/0266-5611/18/4/314 Google Scholar CrossRef Search ADS   Niu Q., Revil A., 2016. Connecting complex conductivity spectra to mercury porosimetry of sedimentary rocks, Geophysics , 81, E17– E32. https://doi.org/10.1190/geo2015-0072.1 Google Scholar CrossRef Search ADS   Niu Q., Prasad M., Revil A., Saidian M., 2016. Textural control on the quadrature conductivity of porous media, Geophysics , 81, E297– E309. https://doi.org/10.1190/geo2015-0715.1 Google Scholar CrossRef Search ADS   Osterman G., Keating K., Binley A., Slater L., 2016. A laboratory study to estimate pore geometric parameters of sandstones using complex conductivity and nuclear magnetic resonance for permeability prediction, Water Resour. Res. , 52, 4321– 4337. https://doi.org/10.1002/2015WR018472 Google Scholar CrossRef Search ADS   Park S.K., Van G.P., 1991. Inversion of pole–pole data for 3-D resistivity structure beneath arrays of electrodes, Geophysics , 56, 951– 960. https://doi.org/10.1190/1.1443128 Google Scholar CrossRef Search ADS   Perrier E., Rieu M., Sposito G., de Marsily G., 1996. Models of the water retention curve for soils with a fractal pore size distribution, Water Resour. Res. , 32, 3025– 3031. https://doi.org/10.1029/96WR01779 Google Scholar CrossRef Search ADS   Pidlisecky A., Haber E., Knight R., 2007. RESINVM3D: a 3D resistivity inversion package, Geophysics , 72, H1– H10. https://doi.org/10.1190/1.2402499 Google Scholar CrossRef Search ADS   Revil A., 2012. Spectral induced polarization of shaly sands: influence of the electrical double layer, Water Resour. Res. , 48, W02517, doi:10.1029/2011WR011260. https://doi.org/10.1029/2011WR011260 Revil A., 2013. Effective conductivity and permittivity of unsaturated porous materials in the frequency range 1 mHz–1 GHz, Water Resour. Res. , 49, 306– 327. https://doi.org/10.1029/2012WR012700 Google Scholar CrossRef Search ADS PubMed  Revil A., Glover P.W.J., 1998. Nature of surface electrical conductivity in natural sands, sandstones, and clays, Geophys. Res. Lett. , 25, 691– 694. https://doi.org/10.1029/98GL00296 Google Scholar CrossRef Search ADS   Revil A., Koch K., Holliger K., 2012. Is it the grain size or the characteristic pore size that controls the induced polarization relaxation time of clean sands and sandstones?, Water Resour. Res. , 48, W05602, doi:10.1029/2011WR011561. https://doi.org/10.1029/2011WR011561 Revil A., Florsch N., Camerlynck C., 2014. Spectral induced polarization porosimetry, Geophys. J. Int. , 198, 1016– 1033. https://doi.org/10.1093/gji/ggu180 Google Scholar CrossRef Search ADS   Revil A., Binley A., Mejus L., Kessouri P., 2015. Predicting permeability from the characteristic relaxation time and intrinsic formation factor of complex conductivity spectra, Water Resour. Res. , 51, 6672– 6700. https://doi.org/10.1002/2015WR017074 Google Scholar CrossRef Search ADS   Revil A. et al.  , 2017. Complex conductivity of soils, Water Resour. Res. , 53, 7121– 7147. https://doi.org/10.1002/2017WR020655 Google Scholar CrossRef Search ADS   Schwarz G., 1962. A theory of the low-frequency dielectric dispersion of colloidal particles in electrolyte solution, J. Phys. Chem. , 66, 2636– 2642. https://doi.org/10.1021/j100818a067 Google Scholar CrossRef Search ADS   Scott J.B.T., Barker R.D., 2003. Determining pore-throat size in Permo-Triassic sandstones from low-frequency electrical spectroscopy, Geophys. Res. Lett. , 30(9), 1450, doi:10.1029/2003GL016951. https://doi.org/10.1029/2003GL016951 Shi J.-Q., Xue Z., Durucan S., 2011. Supercritical CO2 core flooding and imbibition in Berea sandstone—CT imaging and numerical simulation, Energy Procedia , 4, 5001– 5008. https://doi.org/10.1016/j.egypro.2011.02.471 Google Scholar CrossRef Search ADS   Tarantola A., Valette B., 1982. Generalized nonlinear inverse problems solved using the least squares criterion, Rev. Geophys. , 20, 219– 232. https://doi.org/10.1029/RG020i002p00219 Google Scholar CrossRef Search ADS   Timur A., Hempkins W.B., Weinbrandt R.M., 1971. Scanning electron microscope study of pore systems in rocks, J. geophys. Res. , 76, 4932– 4948. https://doi.org/10.1029/JB076i020p04932 Google Scholar CrossRef Search ADS   Titov K., Tarasov A., Ilyin Y., Seleznev N., Boyd A., 2010. Relationships between induced polarization relaxation time and hydraulic properties of sandstone, Geophys. J. Int. , 180, 1095– 1106. https://doi.org/10.1111/j.1365-246X.2009.04465.x Google Scholar CrossRef Search ADS   Titterington D., 1985. General structure of regularization procedures in image reconstruction, Astron. Astrophys. , 144, 381– 387. Tong M., Li L., Wang W., Jiang Y., 2006. Determining capillary-pressure curve, pore-size distribution, and permeability from induced polarization of shaley sand, Geophysics , 71, N33– N40. https://doi.org/10.1190/1.2195989 Google Scholar CrossRef Search ADS   Tuli A., Kosugi K., Hopmans J.W., 2001. Simultaneous scaling of soil water retention and unsaturated hydraulic conductivity functions assuming lognormal pore-size distribution, Adv. Water Resour. , 24, 677– 688. https://doi.org/10.1016/S0309-1708(00)00070-1 Google Scholar CrossRef Search ADS   Vervoort R., Cattle S., 2003. Linking hydraulic conductivity and tortuosity parameters to pore space geometry and pore-size distribution, J. Hydrol. , 272, 36– 49. https://doi.org/10.1016/S0022-1694(02)00253-6 Google Scholar CrossRef Search ADS   Volkmann J., Klitzsch N., 2015. Wideband impedance spectroscopy from 1mHz to 10 MHz by combination of four- and two-electrode methods, J. Appl. Geophys. , 114, 191– 201. https://doi.org/10.1016/j.jappgeo.2015.01.012 Google Scholar CrossRef Search ADS   Volkmann J., Klitzsch N., 2016. Evaluation of low frequency polarization models using well characterized sintered porous glass samples, J. Appl. Geophys. , 124, 39– 53. https://doi.org/10.1016/j.jappgeo.2015.11.011 Google Scholar CrossRef Search ADS   Weigand M., Kemna A., 2016. Debye decomposition of time-lapse spectral induced polarisation data, Comput. Geosci. , 86, 34– 45. https://doi.org/10.1016/j.cageo.2015.09.021 Google Scholar CrossRef Search ADS   Wildenschild D., Sheppard A.P., 2013. X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems, Adv. Water Resour. , 51, 217– 246. https://doi.org/10.1016/j.advwatres.2012.07.018 Google Scholar CrossRef Search ADS   Zhuravlev L.T., 2000. The surface chemistry of amorphous silica. Zhuravlev model, Colloids Surf. A , 173, 1– 38. https://doi.org/10.1016/S0927-7757(00)00556-2 Google Scholar CrossRef Search ADS   APPENDIX A: RELAXATION TIME OF EDL POLARIZATION Consider a spherical surface Γ separating two domains Ωi and Ωe (Fig. A1). The charge density of the surface is QS, and the radius of the curvature of the surface is rc. Apply an oscillating electric field E with its magnitude in the form of E = E0eiωt where ω is the angular frequency. Define spherical coordinates r and θ with r = 0 at the centre of domain Ωi and θ = 0 in the direction of E. Figure A1. View largeDownload slide A spherical surface Γ separating two domains Ωi and Ωe. The radius of the curvature of the surface is rc. Figure A1. View largeDownload slide A spherical surface Γ separating two domains Ωi and Ωe. The radius of the curvature of the surface is rc. The potential distribution in the domains φi and φe follows Laplace equation,   $${\nabla ^2}\ {\varphi ^i} = \ 0 \left( {r < {r_c}} \right)$$ (A1)and   $${\nabla ^2}\ {\varphi ^e} = \ 0 \left( {r > {r_c}} \right).$$ (A2) Define the potential at the surface Γ as φS. At least, eqs (A1) and (A2) satisfy the following boundary conditions: φe → −Ercos θ as r → ∞, φi is finite as r → 0, and $${\varphi _S} = \mathop {\lim }\limits_{r = {r_c}} {\varphi ^e}{\rm{\ }} = \mathop {\lim }\limits_{r = {r_c}} {\varphi ^i}{\rm{\ }}$$. The solutions to eqs (A1) and (A2) are well known and can be expressed as   $${\varphi ^i} = \sum \left( {A_n^i{r^n} + B_n^i\frac{1}{{{r^{n + 1}}}}} \right)\ {P_n}\left( {\cos \theta } \right)$$ (A3)and   $${\varphi ^e} = \sum \left( {A_n^e{r^n} + B_n^e\frac{1}{{{r^{n + 1}}}}} \right)\ {P_n}\left( {\cos \theta } \right)$$ (A4)where $$A_n^i$$, $$B_n^i$$, $$A_n^e$$, $$B_n^e$$ are unknowns in the coefficients for the nth term of Legendre polynomials Pn(cos θ). Applying boundary condition (1), eq. (A4) can be simplified as   $${\varphi ^e} = \left( {A_1^er + B_1^e\frac{1}{{{r^2}}}} \right)\ \cos \theta .$$ (A5)Then, the potential at the surface can be obtained as   $${\varphi _S} = \left( {A_1^e{r_c} + B_1^e\frac{1}{{{r_c}^2}}} \right)\ \cos \theta .$$ (A6) If we expand φS using Legendre polynomials, φS has the following form:   $${\varphi _S} = \sum {\beta _n}{P_n}\left( {\cos \theta } \right),$$ (A7)where βn is the associated coefficient. Comparing eqs (A6) and (A7), we find that only the first term (n = 1) in eq. (A7) is non-zero. This result will be used later to determine the relaxation time of the change in QS in the surface Γ. Assume the charge on Γ (a thin EDL) only moves tangentially. The movement of charge on Γ is due to two driving forces: (1) electric field and (2) ion concentration gradient. Considering the continuity, the resulting time derivative of the surface charge density can be written as   $$\frac{{\partial {Q_S}}}{{\partial t}} = {e_0}{\rm{\ }}\nabla \cdot \left( {{j_E} + {j_D}} \right),$$ (A8)where e0 is the charge carried by an ion in the EDL, and jE and jD are surface fluxes of ions induced by electric field and ion diffusion, respectively. Consider the change in QS is ΔQS, which is much less than QS. Also, consider the expressions of jE and jD. Then, eq. (A8) can be expressed as   \begin{eqnarray}i\omega \ \Delta {Q_S} \!=\! \frac{{\mu kT}}{{{r_c}^2}}\ \frac{1}{{\sin \theta }}\frac{{\rm{d}}}{{{\rm{d}}\theta }}\left( {\sin \theta \frac{{{\rm{d}}\Delta {Q_S}}}{{{\rm{d}}\theta }} \!+\! \frac{{{e_0}{Q_S}}}{{kT}}\sin \theta \frac{{{\rm{d}}{\varphi _S}}}{{{\rm{d}}\theta }}} \right),\end{eqnarray} (A9)where μ is the mobility of ions in EDL, k is Boltzmann constant and T is the absolute temperature. We also expand the surface charge variation ΔQS using Legendre polynomials as   $$\Delta {Q_S} = \sum {\alpha _n}{P_n}\left( {\cos \theta } \right).$$ (A10) Insert eqs (A7) and (A10) into eq. (A9), and after some derivative calculations, eq. (A9) becomes   \begin{eqnarray}i\omega \ \sum {\alpha _n}{P_n} &=& \frac{{\mu kT}}{{{r_c}^2}}\ \sum {\alpha _n}\left( {\frac{{\cos \theta }}{{\sin \theta }}P_n^{\prime} + P_n^{\prime\prime}} \right) \nonumber\\ &&+ \frac{{\mu {e_0}{Q_S}}}{{{r_c}^2}}\sum {\beta _n}\left( {\frac{{\cos \theta }}{{\sin \theta }}P_n^{\prime} + P_n^{\prime\prime}} \right).\end{eqnarray} (A11) Considering the following equality,   $$P_n^{\prime\prime} + \frac{{\cos \theta }}{{\sin \theta }}P_n^\prime + n\left( {n + 1} \right)\ {P_n} = \ 0,$$ (A12)eq. (A11) can be rewritten as   \begin{eqnarray}\sum {\alpha _n}\left( {i\omega \!+\! \frac{{\mu kT}}{{{r_c}^2}}n\left( {n + 1} \right)} \right){P_n} \!=\! \ \sum {\beta _n}\frac{{\mu {e_0}{Q_S}}}{{{r_c}^2}}n\left( {n + 1} \right){P_n}.\nonumber\\ \end{eqnarray} (A13) Because of the orthogonality of Legendre polynomials, the following relation is valid for all the terms,   $${\alpha _n} = \ - \frac{1}{{1 + i\omega \frac{{{r_c}^2}}{{n\left( {n + 1} \right)\mu kT}}}}\frac{{{e_0}{Q_S}}}{{kT}}{\beta _n}.$$ (A14) As discussed, the coefficient βnfor the surface potential φS is non-zero only when n = 1, we have the following relation between ΔQS and φS:   $$\Delta {Q_S} = \ - \frac{1}{{1 + i\omega \tau }}\frac{{{e_0}{Q_S}}}{{kT}}{\varphi _S}$$ (A15) with   $$\tau \ = \frac{{{r_c}^2}}{{2\mu kT}}\ = \frac{{{r_c}^2}}{{2{D^ + }}}\$$ (A16)where D+ is the diffusion coefficient of the ions in the EDL. Eq. (A16) states that the surface charge density follows a simple Debye relaxation with relaxation time τ, which is controlled by the radius of the curvature of the surface Γ. In the derivation, the electrical properties of domains Ωi and Ωe were not specified, and therefore, eq. (A15) holds for both following cases: (1) domain Ωi is solid particle and domain Ωe is electrolyte; and (2) domain Ωi is electrolyte and domain Ωe is solid particle. In the former case, the radius of the curvature of the surface Γ is equal to the particle radius; for the latter, rc is equal to the pore size. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Joint inversion of NMR and SIP data to estimate pore size distribution of geomaterials

, Volume 212 (3) – Mar 1, 2018
15 pages

/lp/ou_press/joint-inversion-of-nmr-and-sip-data-to-estimate-pore-size-distribution-DY67tBjzjO
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx501
Publisher site
See Article on Publisher Site

### Abstract

Summary There are growing interests in using geophysical tools to characterize the microstructure of geomaterials because of the non-invasive nature and the applicability in field. In these applications, multiple types of geophysical data sets are usually processed separately, which may be inadequate to constrain the key feature of target variables. Therefore, simultaneous processing of multiple data sets could potentially improve the resolution. In this study, we propose a method to estimate pore size distribution by joint inversion of nuclear magnetic resonance (NMR) T2 relaxation and spectral induced polarization (SIP) spectra. The petrophysical relation between NMR T2 relaxation time and SIP relaxation time is incorporated in a nonlinear least squares problem formulation, which is solved using Gauss–Newton method. The joint inversion scheme is applied to a synthetic sample and a Berea sandstone sample. The jointly estimated pore size distributions are very close to the true model and results from other experimental method. Even when the knowledge of the petrophysical models of the sample is incomplete, the joint inversion can still capture the main features of the pore size distribution of the samples, including the general shape and relative peak positions of the distribution curves. It is also found from the numerical example that the surface relaxivity of the sample could be extracted with the joint inversion of NMR and SIP data if the diffusion coefficient of the ions in the electrical double layer is known. Comparing to individual inversions, the joint inversion could improve the resolution of the estimated pore size distribution because of the addition of extra data sets. The proposed approach might constitute a first step towards a comprehensive joint inversion that can extract the full pore geometry information of a geomaterial from NMR and SIP data. Electrical properties, Magnetic properties, Microstructure, Hydrogeophysics, Joint inversion 1 INTRODUCTION The pore space plays a fundamental role in controlling fluid transport in porous media (Dullien 2012). For geomaterials such as soils and sedimentary rocks, the geometric characteristic of the pore space is usually described by porosity, specific surface area, pore size and its distribution. These parameters have been frequently used in estimating hydraulic properties of geomaterials. For instance, several characteristic pore sizes have been found to be correlated with the saturated permeability of sandstones (Johnson et al.1986; Katz & Thompson 1986); the pore size distribution is an important input for predicting water retention curves (Kosugi 1994; Perrier et al.1996) and unsaturated permeability of soils and rocks (Kosugi 1996; Tuli et al.2001; Vervoort & Cattle 2003). The geometric characteristics of pore space can be directly measured using some imaging methods, for example, X-ray imaging technique (see a review in Wildenschild & Sheppard 2013), scanning electron microscope (e.g. Timur et al.1971) and laser scanning confocal microscopy (e.g. Fredrich 1999). In addition, there are also many indirect methods that can be used to estimate the geometric characteristics of the pore space. One example is the mercury injection capillary porosimetry (MICP), and the measurement is usually interpreted in terms of pore throat size rather than pore size as demonstrated in Dullien (1981); another example is the gas adsorption analysis, which has been frequently used to infer the surface area and pore size of rock samples (e.g. Kuila & Prasad 2013). The use of geophysical techniques becomes increasingly attractive in estimating the geometric characteristics of the pore space in geomaterials (Tong et al.2006; Mohnke & Yaramanci 2008; Revil et al.2014), mainly due to its non-invasive nature and high spatiotemporal resolution. Another advantage of using geophysical tools is the applicability to field characterizations. The most promising geophysical tools in practice include nuclear magnetic resonance (NMR) and spectral induced polarization (SIP). The traditional NMR method measures the magnetic response of protons under specific sequences of applied magnetic fields (see Behroozmand et al.2015 for a review). In an NMR relaxation measurement, the signal magnitude is proportional to the absolute number of protons and thus the volume of hydrogen bearing fluid in the sample. The NMR decay curve can be decomposed to form a distribution of relaxation time (Borgia et al.1998), which is related to surface-to-volume ratio (Brownstein & Tarr 1979), thus the effective pore size (and its distribution) of the material can be estimated. The SIP method measures the electrical response of a material under an external alternating electric field, and the measurements can be expressed in terms of complex conductivity (see Kemna et al.2012 for a review). At low frequencies, the variation in complex conductivity of a geomaterial is mainly controlled by the polarization of the electrical double layer (EDL) forming at the solid–liquid interface (Lesmes & Morgan 2001). The measured complex conductivity spectra can be decomposed to obtain a distribution of relaxation time (Florsch et al.2014). Since the relaxation time is related to the geometry of the EDL (Schwarz 1962), the obtained relaxation time distribution thus can give an estimation for the pore size distribution of the material. For simplicity, in this study we use ‘NMR porosimetry’ and ‘SIP porosimetry’ to term the methods of using NMR and SIP data to estimate the pore size distribution of porous media, respectively. The geometric control of pore space on NMR and SIP responses of geomaterials is quite complicated, and there are still many ongoing studies towards a full understanding of the petrophysical relations (see some recent researches in Keating & Falzone 2013; Müller-Petke et al.2015; Niu et al.2016; Volkmann & Klitzsch 2016). In practice, some simplifications are usually made such that the pore size information of the material can be inferred from the NMR and SIP measurements. For instance, the NMR porosimetry assumes (1) that the magnetic relaxation of protons in porous media is in fast diffusion regime and (2) that the effects of bulk fluid diffusion and pore coupling on NMR relaxation are negligible (see a review in Behroozmand et al.2015). However, it has been found that these assumptions may be violated for materials with large pores and/or large values of surface relaxivity (Grunewald & Knight 2009; Dlubac et al.2014). In SIP porosimetry, it is generally assumed that the complex conductivity of geomaterials at frequencies lower than 100 kHz is solely from the EDL polarization and dielectric polarization (Leroy et al.2008; Revil et al.2014). However, other polarization may also contribute to the measured complex conductivity, such as interfacial polarization (Hanai 1960) and membrane polarization (Bücker & Hördt 2013). These approximate petrophysical relations used in NMR and SIP porosimetries may induce significant errors to the estimated pore size distribution. Inversion is needed to decompose the NMR relaxation curve and SIP spectra to obtain the relaxation time distribution. As is well known, most geophysical inversions are underdetermined due to the limited number of measurements; moreover, the inversion is also ill-posed because the measurements are usually contaminated by noise (Tarantola & Valette 1982). To mitigate these problems, joint inversion has become a common approach (e.g. Moorkamp et al.2011) because extra information from other data set could reduce the ambiguity or non-uniqueness of the inversion (Haber & Oldenburg 1997; Linde & Doetsch 2016). This has inspired us to consider the joint use of NMR and SIP data in estimating the pore size distribution of geomaterials. One similar attempt in this research direction is the work of Osterman et al. (2016) where both NMR and SIP measurements were used to estimate the pore geometric parameters and thus the permeability of sandstones. However, in their study different types of data sets are still processed separately, and a joint inversion of NMR and SIP data is still not realized yet. In this study, we explore the idea of joint inversion of NMR and SIP data for the estimation of pore size distributions of geomaterials. We first briefly review theoretical backgrounds of NMR and SIP porosimetries with emphasis on commonly applied petrophysical models. We then introduce a joint inversion approach that can simultaneously process NMR and SIP data. In the joint inversion, one aspect of significant importance is the trade-off between different data sets. We propose a two-step method to balance the weight between NMR and SIP measurements. A numerical experiment is designed to demonstrate the applicability of the proposed inversion method. In the last part of the study, we apply this joint inversion to a Berea sandstone sample and discuss the results as well as the limitations of our proposed approach. 2 THEORETICAL BACKGROUND 2.1 NMR relaxation in porous media In water-saturated geomaterials, NMR relaxation phenomena arises as the result of the interaction of diffusing nuclear spins of protons with their surrounding environment. The NMR relaxation can produce a bulk magnetization, which is governed by the phenomenological Bloch–Torrey equations. Experimentally, the transverse component of the magnetization can be detected using a series of oscillating magnetic field pulses [e.g. the Carr–Purcell–Meiboom–Gill (CPMG) pulse sequence] and the measured signal is in the form of a voltage decay curve, Exy(t). The solution to the diffusion problem can be expressed as the sum of a number of exponential decays (Brownstein & Tarr 1979) and thus the NMR signal can be expressed as   $${E_{xy}}\ \left( t \right) = {E_0}\ \mathop \sum \nolimits_i^{} {f_i}{e^{ - \frac{t}{{T_2^i}}}}$$ (1)where Exy is the transverse component of the NMR signal, E0 is the initial signal magnitude, fi is the relative intensities of the magnetic field relaxing with the relaxation time of $$T_2^i$$. The initial magnitude E0 is directly proportional to the number of protons and thus the amount of water in the material and therefore it can be used to estimate the water content of the material. Assuming a uniform magnetic field, the relaxation with time constant T2i only involves two processes (Brownstein & Tarr 1979) with   $$\frac{1}{{T_2^i}} = \frac{1}{{{T_{2B}}}}\ + \frac{1}{{T_{2S}^i}}$$ (2)where relaxation time T2B is related to the bulk fluid diffusion and $$T_{2S}^i$$ is related to surface relaxation. The diffusion rate of protons in bulk water in geomaterials is small, ∼0.3 s−1; the relaxation can be greatly enhanced at the mineral surface due to the presence of paramagnetic impurities such as Mn2+ and Fe3+ (Foley et al.1996). Compared with bulk diffusion, the surface relaxation could dominate the relaxation process in a saturated geomaterial such that $$\frac{1}{{{T_{2B}}}}$$ in eq. (2) could be dropped. Assuming no coupling between pores (e.g. see Grunewald & Knight 2009), a relaxation with time constant $$T_2^i$$ in eq. (1) can be related to a pore by   $$\frac{1}{{T_2^i}} = \frac{1}{{T_{2S}^i}}\ = \frac{1}{{\frac{{{r_i}}}{{\alpha {\rho _s}}} + \frac{{r_i^2}}{{2\alpha D}}}}\$$ (3)where ri is the effective size of the pore, α describes the shape of the pore (α = 3 for spheres), and D is the diffusion coefficient of water. The parameter ρs, termed as surface relaxivity, qualifies the efficiency of the surface relaxation. In practice, the term $$\frac{{{r^2_i}}}{{2\alpha D}}$$ in eq. (3) is usually ignored (i.e. assuming in fast diffusion regime), and the relaxation time can be related to the pore size by   $$T_2^i = \frac{{{r_i}}}{{\alpha {\rho _s}}}\ .$$ (4) Thus, the relative intensities fi of the relaxation time $$T_2^i$$ can be interpreted as pore size distribution. It should be noted that the bulk diffusion term in eq. (2) could be significant if the material is characterized with large pore space and low magnetic susceptibility (e.g. see Grunewald & Knight 2011). In that case, eq. (4) might be insufficient in describing the T2 relaxation time, and thus the bulk diffusion rate should also be included. 2.2 SIP and EDL polarization in porous media The SIP measures the complex conductivity response σ* of a porous medium at frequencies lower than 100 kHz. While the real part of the complex conductivity σ΄ quantifies the material's ability to conduct an electric current, the imaginary part σ″ reflects the storage capability of electrical energy from the applied electric field. The two parameters are not independent; they are related to each other by the Kramers–Kronig relation (e.g. see Volkmann & Klitzsch 2015). In this study, our discussion will be focused on the imaginary conductivity. The low-frequency imaginary conductivity of geomaterials is usually interpreted in the framework of EDL polarization. In geomaterials, the mineral surface is usually negatively charged because of isomorphous substitutions (McBride 1989) and/or the presence of hydroxyl group (Zhuravlev 2000). When the material is saturated with an electrolyte (e.g. NaCl solution), some ions (mostly cations) will accumulate at the liquid-solid interface to balance the surface charge, forming the EDL. Consider a spherical pore space saturated with NaCl solution. After applying an electric field, the ions in the EDL will redistribute to gain electric dipole moment. For the spherical EDL, the induced complex conductance follows Debye relaxation (Schwarz 1962),   $${{C}}_i^* = \frac{{i\omega \tau }}{{1 + i\omega \tau }}\ {{\rm{\Sigma }}_S},$$ (5)where ΣS is the surface conductance related to Stern layer. The relaxation time constant τ is related to the radius of the curvature of the EDL rc (see Appendix  A for details) by   $$\tau \ = \frac{{{r_c}^2}}{{2{D^ + }}}\$$ (6)where D+ is the diffusion coefficient of the ions in the EDL. Note that eq. (5) does not consider the ion exchange between EDL and bulk solution. For colloids, the EDL is coating the particles, and therefore, rc is equal the particle radius rd as shown in Schwarz (1962). For a spherical pore embedded in rock matrix, the radius of the curvature of the EDL rc is equal to the pore radius r, and therefore the relaxation time τ is   $$\tau \ = \frac{{{r^2}}}{{2{D^ + }}}\ .$$ (7) In geomaterials, the shape of the EDL seems to be controlled by the pore space and thus eq. (7) applies. The correlation between SIP relaxation time and pore size have been experimentally observed for a number of geological materials (Binley et al.2005; Kruschwitz et al.2010; Titov et al.2010; Revil et al.2012). For porous media with different pore sizes, the contribution from each pore can be summed up to form the overall response (e.g. see Leroy et al.2008),   $${{C}}_{{\rm{EDL}}}^* = \mathop \sum \nolimits_i^{} {f_i}\frac{{i\omega {\tau _i}}}{{1 + i\omega {\tau _i}}}\ {{\rm{\Sigma }}_S},$$ (8)where $${{C}}_{{\rm{EDL}}}^*$$ is the complex conductance induced by EDL polarization, and fi is the relative weighting of the pore with size ri, which corresponds to the relaxation time τi. For geological materials, the (complex) conductance induced by EDL (in S) is usually treated as a perturbation of the pore fluid conductivity (in S m−1) by considering Joule dissipation (Johnson et al.1986; Niu et al.2016). The increased complex conductivity of the pore fluid is sometimes termed as specific surface (complex) conductivity $${\rm{\sigma }}_{ss}^*$$, and it has the following form,   $${\rm{\sigma }}_{ss}^* = \frac{{2\left( {{{C}}_{{\rm{EDL}}}^* + C_S^0} \right)}}{{\rm{\Lambda }}}\ ,$$ (9)where $$C_S^0$$ is the DC surface conductance (e.g. see Revil & Glover 1998) and Λ is a characteristic pore size (Johnson et al.1986), which can be related to surface area-to-pore volume Spor and formation factor F by   $$\frac{2}{{\rm{\Lambda }}} = m{S_{{\rm{por}}}}$$ (10)where m is the cementation factor. After applying Archie's law, the complex conductivity of the material can be obtained as   $${{\rm{\sigma }}^*} = \frac{1}{F}\ \left( {{{\rm{\sigma }}_f} + {\rm{\sigma }}_{ss}^*} \right)$$ (11)where σf is the fluid conductivity. Inserting eqs (9) and (10) into eq. (11) yields   $${{\rm{\sigma }}^*} = \frac{{{{\rm{\sigma }}_f}}}{F}\ + \frac{{m{S_{{\rm{por}}}}}}{F}\left( {\mathop \sum \nolimits_i^{} {f_i}\frac{{i\omega {\tau _i}}}{{1 + i\omega {\tau _i}}}{{\rm{\Sigma }}_S} + C_S^0} \right).$$ (12) In eq. (12), the first term in the right-hand side (a real number) is associated with ionic conduction occurring in the bulk fluid in the pore space. The real part of the second term related to the surface conduction occurring in the EDL (Revil & Glover 1998) and the imaginary part quantifies the EDL polarization. The parameters ΣS and and $$C_S^0$$ can be expressed in terms of fluid chemistry and electrochemical properties of the EDL (e.g. see Revil 2012). In this study, we will not address the physical origin of the petrophysical parameters in eq. (12) such as F, m, ΣS and $$C_S^0$$; instead, we simplify eq. (12) by using several widely used lumped parameters. If the permittivity induced by dielectric polarization εinf is considered, the above equation is reformulated as,   $${{\rm{\sigma }}^*} = {\sigma _\infty }\ - \left( {{\sigma _\infty } - {\sigma _0}} \right)\mathop \sum \nolimits_i^{} \frac{{{f_i}}}{{1 + i\omega {\tau _i}}} + i\omega {\varepsilon _{{\rm{inf}}}},$$ (13)where σ∞ and σ0 are the real conductivities of the sample at low- and high-frequency limits, respectively and their relations to fluid chemistry and EDL properties can be found, for example, in Revil (2012). Since both σ∞ and σ0 are optimized during the inversion, to a certain extent, the effects of surface charge density, fluid salinity, and pH are also considered in interpreting the complex conductivity data. 2.3 Petrophysical simplifications in NMR and SIP porosimetries In this subsection, we summarize the petrophysical assumptions made in NMR and SIP porosimetries. In NMR porosimetry, it is assumed that the effect of diffusion on NMR relaxation in a magnetic gradient is negligible . The magnetic gradient arises from the magnetic susceptibility contrast between grain minerals and pore fluid (Cotts et al.1989). This assumption is usually valid if the NMR T2 relaxation is measured with CPMG pulse sequence with short echo times (Kleinberg & Horsfield 1990). Another assumption in NMR porosimetry is that the pores are isolated and there is no pore coupling effect (D’Orazio et al.1989). For most consolidated sandstones, the narrow pore throat could restrict the molecular diffusion between two interconnected pore spaces during the time of an NMR experiment, and therefore, this assumption is approximately valid (Borgia et al.1996). However, for samples with well-connected heterogeneous pores such as soil aggregates and analogous synthetic silica gels, the assumption might break down (Grunewald & Knight 2009). In addition, the NMR porosimetry assumes that T2 relaxation is controlled by surface relaxation in the fast-diffusion regime, in which the diffusion-controlled exchange between surface and bulk water is fast enough such that all the water molecules interact with the surface during the NMR measurement (Bryar et al.2000). This fast diffusion regime assumption may not be appropriate for materials with high surface relaxivity and/or large pores (e.g. Keating & Falzone 2013). In SIP porosimetry, the measurements are usually interpreted with the EDL polarization model. The most important assumption is that the pore size r is equal to the radius of EDL curvature rc so eq. (7) holds. For isolated pores, this assumption is clearly valid. For geomaterials with interconnected pores, although the correlation between r and τ has been experimentally observed (Revil et al.2012), it remains still unclear whether the relation is quadratic (Scott & Barker 2003; Titov et al.2010). Second, eq. (7) also assumes there is no ion exchange between EDL and bulk fluid (Lyklema et al.1983; Chelidze & Gueguen 1999; Niu & Revil 2016). As a matter of fact, the ion exchange could speed up the ion redistribution in the EDL and thus decrease the relaxation time τ as shown in Lyklema et al. (1983). In SIP porosimetry, the interfacial polarization and dielectric effects may also be considered. However, the theory combing EDL polarization, interfacial polarization, and dielectric polarization still significantly underestimates the measured permittivity of sandstones (Lesmes & Morgan 2001; Kremer et al.2016). Lastly, although not stated explicitly, SIP porosimetry applies the isolated pore model, which assumes the ion movement is limited to a single pore or EDL. That is, no ion exchange occurs between two adjacent pores during a SIP measurement. 3 JOINT INVERSION APPROACH As discussed, NMR and SIP porosimetries involve a number of assumptions about the petrophysical relations of a geomaterial. Since it is difficult to justify these assumptions for a real geomaterial, the pore size information obtained separately from NMR and SIP data may not be accurate. To reduce the uncertainties and thus increase the resolution, in this study we suggest to jointly process the NMR and SIP measurements for the estimation of the pore size distribution of geomaterials. The joint inversion scheme is illustrated in Fig. 1, and each step of the joint inversion is discussed in detail below. Figure 1. View largeDownload slide The joint inversion scheme to estimate pore size distribution from NMR and SIP measurements. Figure 1. View largeDownload slide The joint inversion scheme to estimate pore size distribution from NMR and SIP measurements. 3.1 Petrophysical relations In practice, there are two types of approaches for joint inversion of two geophysical measurements: the petrophysical approach (e.g. Jardani & Revil 2009) and structural approach (e.g. Gallardo & Meju 2003). The petrophysical approach assumes that different geophysical measurements are functions of the same properties of a material. Thus, multiple geophysical data can be integrated in the inversion. For instance, it assumes both resistivity and seismic velocity of a rock are controlled by its porosity and saturation, and thus resistivity and seismic data can be jointly inverted (Gao et al.2012). Alternatively, the structural approach assumes different geophysical methods yield structurally similar geophysical images of the subsurface. That is, the spatial and directional variations of different geophysical responses are similar. Mathematically, the structural difference between two images can be quantified using the cross-gradient function, and it can be combined with data misfit and smoothness constraints to form the objective function in the joint inversion (Gallardo & Meju 2003; Gallardo & Meju 2004). In this study, we apply the petrophysical approach in the joint inversion of NMR and SIP measurements. We assume both NMR T2 relaxation time and SIP relaxation time τ are related to the pore size r, that is, eqs (4) and (7) are valid. This implies that, for samples with similar surface properties (e.g. D+ and ρs), the relaxation time τ is proportional to squared T2. We verify the relation between τ and T2 using Sherwood sandstone samples from Osterman et al. (2016). The T2 squared and τ of these samples are cross-plotted in Fig. 2 where T2 is the peak NMR relaxation time and τ is determined from the characteristic frequency chosen using the approach in Revil et al. (2015). Despite an outliner point, a strong linear correlation is found between T22 and τ. The experimental evidence lends additional support to the analytical relation between T2 and τ as presented in the previous section. Figure 2. View largeDownload slide The relation between squared NMR T2 relaxation time and SIP relaxation time τ of Sherwood sandstone samples. The data are from Osterman et al. (2016). The datum with empty triangle is an outlier point. Figure 2. View largeDownload slide The relation between squared NMR T2 relaxation time and SIP relaxation time τ of Sherwood sandstone samples. The data are from Osterman et al. (2016). The datum with empty triangle is an outlier point. 3.2 Model parameters The common relation to pore size r makes it possible to jointly invert the pore size distribution [e.g. the relative intensities/weighting fi in eqs (1) and (13)] from NMR and SIP measurements. In addition to fi, other parameters related to the petrophysical models may also be inverted during the inversion, including σ∞, σ0, εinf, and E0 (e.g. Florsch et al.2014). The surface properties D+ and ρs are also involved in the inversion, but they cannot be recovered in separate inversions because of the non-uniqueness in $$\frac{r}{{3{\rho _s}}}$$ and $$\frac{{{r^2}}}{{2{D^ + }}}$$, as nicely demonstrated by Mohnke (2014) for ρs. As shown in the previous section and Fig. 2, relaxation time T2 and τ are related to each other for a given material through the pore size. This implies that the surface properties D+ and ρs are not independent, and therefore it is possible to recover one parameter once the other is known. In fact, previous research has shown, with extra data sets, the joint inversion can properly estimate the material's ρs (Mohnke 2014). In this study, we assume D+ is known as suggested by Revil (2013). The surface relaxivity ρs in the inversion can be either assumed as a constant (if we have a prior information on it) or updated iteratively. We define a model vector m, which contains all the unknowns to be recovered from the joint inversion. The model vector m can be expressed as   $${{\bf m}}\ = \left( {\begin{array}{@{}*{1}{c}@{}} {{\bf x}}\\ {\ln {\sigma _0}}\\ {\begin{array}{@{}*{1}{c}@{}} {\begin{array}{@{}*{1}{c}@{}} {\ln {\sigma _\infty }}\\ {\ln {\varepsilon _{{\rm{inf}}}}} \end{array}}\\ {\ln {E_0}}\\ {\ln {\rho _s}} \end{array}} \end{array}} \right)\ ,$$ (14)where x is a vector containing the relative weighting fi for the considered pore size ri in the form of xi = ln fi. The use of log-transformed pore size distribution (relative weighting) and other parameters in eq. (14) is to ensure positivity of the unknowns. Moreover, since the relative weighting fi of some pores are considerably smaller than the others, the log-transformation of fi ensures their relative changes are also important during the inversion (Weigand & Kemna 2016). 3.3 Objective function The objective function Ψ consists of data misfit Ψd and a model regularization term Ψm (Aster et al.2013),   $${\rm{\Psi \ }} = {{\rm{\Psi }}_d}\ + \lambda {{\rm{\Psi }}_m},$$ (15)where the damping parameter λ controls the relative contribution of the regularization term Ψm on the objective function. The regularization is needed to stabilize the inversion process to obtain a unique model (Titterington 1985). In NMR and SIP inversions, Tikhonov regularization is commonly used (Florsch et al.2014; Weigand & Kemna 2016). In this study, we also employ the Tikhonov regularization and the l2 regularization term Ψm can be expressed as   $${{\rm{\Psi }}_m} = \left\| {{{{\bf W}}_{{\bf m}}}\left( {{{\bf m}} - {{{\bf m}}_{{\rm{ref}}}}} \right)} \right\|_2^2$$ (16)where ${{{\bf W}}_{{\bf m}}} = ( \scriptsize{\begin{array}{@{}*{2}{c}@{}} {{\bf L}}&0\\ 0&{{\bf I}} \end{array}} ){\boldsymbol{\ }}$ is a regularization matrix and m is a reference model that may contain a priori information. While the discrete operator L in Wm applies smoothness constraints to the pore size distribution x in the model vector m, the identify matrix I forces the petrophysical unknowns (i.e. σ∞, σ0, εinf, E0 and ρs) to be close to reference values defined in mref (a similar treatment can be found in Florsch et al.2014). In this study, the operator L is chosen as the second-order row difference of an identity matrix (i.e. eq. 19 in Weigand & Kemna 2016). The minimization of the objective function Ψ can be realized using an iterative Gauss–Newton scheme such that the model update can be determined (e.g. Park & Van 1991) from   \begin{eqnarray} \left[ {{{{\bf J}}^T}{{\bf W}}_d^T{{{\bf W}}_d}{{\bf J}} + {\rm{\lambda }}{{\bf W}}_m^T{{{\bf W}}_{{m}}}} \right]\Delta {{{\bf m}}_k}& =& {{{\bf J}}^T}\ {{\bf W}}_d^T{{{\bf W}}_d}\left( {{{\bf d}} - {{\bf G}}\left( {{{{\bf m}}_k}} \right)} \right)\nonumber\\ && - {\rm{\lambda }}{{\bf W}}_m^T{{{\bf W}}_{{m}}}\left( {{{{\bf m}}_k} - {{{\bf m}}_{{\rm{ref}}}}} \right) \end{eqnarray} (17)where J is the sensitivity matrix corresponding to the forward modelling function G, which includes both eqs (1) and (13), Wd is the data weighting matrix, Δmk is the model update for the model mk at the iteration step k, and d is the data vector. In the joint inversion, the vector d is expressed as   $${{\bf d}}\ = \left( {\begin{array}{@{}*{1}{c}@{}} {{{\bf d}}_{{\rm{SIP}}}^{{\rm{Re}}}}\\ {{{\bf d}}_{{\rm{SIP}}}^{{\rm{Im}}}}\\ {{{{\bf d}}_{{\rm{NMR}}}}} \end{array}} \right)\$$ (18)where vectors $${{\bf d}}_{{\rm{SIP}}}^{{\rm{Re}}}$$ and $${{\bf d}}_{{\rm{SIP}}}^{{\rm{Im}}}$$ (with the size of m × 1) contain the real part and imaginary part of the frequency dependent complex conductivity, respectively, and the vector dNMR contains the measured NMR T2 relaxation time measurements. Since data points in an NMR experiment is usually at the order of thousands, it will be reduced in the inversion so that the dimension of dNMR (with the size of n × 1) is comparable to $${{\bf d}}_{{\rm{SIP}}}^{{\rm{Re}}}$$ and $${{\bf d}}_{{\rm{SIP}}}^{{\rm{Im}}}$$. The data weighting Wd matrix contains the errors related to the data vector, and it can be expressed as a diagonal matrix,   $${{{\bf W}}_d} = \left( {\begin{array}{c@{\quad}c@{\quad}c} {\beta {\boldsymbol \epsilon} _{{\rm{SIP}}}^{{\rm{Re}}}}&0&0\\ 0&{\beta {\boldsymbol \epsilon} _{{\rm{SIP}}}^{{\rm{Im}}}}&0\\ 0&0&{{{\boldsymbol \epsilon} _{{\rm{NMR}}}}} \end{array}} \right),$$ (19)where $${\boldsymbol \epsilon} _{{\rm{SIP}}}^{{\rm{Re}}} = \ {\rm{diag}}( {\frac{1}{{{\epsilon} _j^{Re}}}} )$$, $${\boldsymbol \epsilon} _{{\rm{SIP}}}^{{\rm{Im}}} = \ {\rm{diag}}( {\frac{1}{{\epsilon _j^{Im}}}} )$$, $${\epsilon _{{\rm{NMR}}}} = \ {\rm{diag}}( {\frac{1}{{{\epsilon _k}}}} )$$, and β is a trade-off parameter balancing the weight between SIP and NMR measurements. The parameters $$\epsilon _j^{\rm Re}$$ and $$\epsilon _j^{\rm Im}$$ are the errors related to the real part and imaginary part of the jth complex conductivity measurements, and εk is the error related to the kth reduced NMR measurement. In practice, these errors can be estimated from the standard deviation of the measurement noise (e.g. Pidlisecky et al.2007). Starting with a reference model mref, the model vector can be updated iteratively until the root-mean-square value of the misfit is smaller than a pre-defined value. At the step k, the model update can be calculated from eq. (17), and thus updated model mk+1 can be determined as   $${{{\bf m}}_{k + 1}} = {{{\bf m}}_k}\ + \chi \Delta {{{\bf m}}_k}$$ (20)where the step parameter χ belongs to [0, 1]. The restriction of χ ≤ 1 prevents the model from overshooting due to non-linearity, and it can be determined using a line search. For details, one can refer to Weigand & Kemna (2016). 3.4 Sensitivity calculation The sensitivity matrix (also known as Jacobian matrix) J in eq. (17) contains the partial derivatives of the model responses with respect to the model parameters, expressed as   \begin{eqnarray} {{\bf J}} \!=\!\! \left( {\begin{array}{@{}*{8}{c}@{}} {\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {f_1}}}}& \cdots &{\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {f_L}}}}&{\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {\sigma _0}}}}&{\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {\sigma _\infty }}}}&{\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {\varepsilon _{{\rm{inf}}}}}}}&{\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {E_0}}}}&{\frac{{\partial \sigma '\left( {{\omega _1}} \right)}}{{\partial \ln {\rho _s}}}}\\ \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {f_1}}}}& \cdots &{\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {f_L}}}}&{\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {\sigma _0}}}}&{\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {\sigma _\infty }}}}&{\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {\varepsilon _{{\rm{inf}}}}}}}&{\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {E_0}}}}&{\frac{{\partial \sigma '\left( {{\omega _m}} \right)}}{{\partial \ln {\rho _s}}}}\\ {\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {f_1}}}}& \cdots &{\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {f_L}}}}&{\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {\sigma _0}}}}&{\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {\sigma _\infty }}}}&{\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {\varepsilon _{{\rm{inf}}}}}}}&{\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {E_0}}}}&{\frac{{\partial \sigma ''\left( {{\omega _1}} \right)}}{{\partial \ln {\rho _s}}}}\\ \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {f_1}}}}& \cdots &{\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {f_L}}}}&{\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {\sigma _0}}}}&{\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {\sigma _\infty }}}}&{\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {\varepsilon _{{\rm{inf}}}}}}}&{\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {E_0}}}}&{\frac{{\partial \sigma ''\left( {{\omega _m}} \right)}}{{\partial \ln {\rho _s}}}}\\ {\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {f_1}}}}& \cdots &{\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {f_L}}}}&{\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {\sigma _0}}}}&{\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {\sigma _\infty }}}}&{\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {\varepsilon _{{\rm{inf}}}}}}}&{\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {E_0}}}}&{\frac{{\partial E\left( {{t_1}} \right)}}{{\partial \ln {\rho _s}}}}\\ \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {f_1}}}}& \cdots &{\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {f_L}}}}&{\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {\sigma _0}}}}&{\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {\sigma _\infty }}}}&{\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {\varepsilon _{{\rm{inf}}}}}}}&{\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {E_0}}}}&{\frac{{\partial E\left( {{t_n}} \right)}}{{\partial \ln {\rho _s}}}} \end{array}} \right),\nonumber\\ \end{eqnarray} (21)where σ'(ωi) and σ''(ωi) are the real and imaginary conductivities of the ith calculated SIP data (i = 1, 2, … m) and E(tj) is the jth calculated NMR data (j = 1, 2, … n). The sensitivity can be calculated by several methods (McGillivray & Oldenburg 1990). In this study, we analytically determine each element of the sensitivity matrix J in eq. (21). For most of the elements in Jacobian matrix J, they are dependent on the model parameters x. To ensure a correct update direction, in the inversion the matrix J is calculated at each iteration step. If one of the petrophysical model parameters is known (e.g. the surface relaxivity ρs is independently measured or estimated), we can set the related sensitivity components (e.g. last column in J in eq. 21) as zero so that the parameter will not be updated in the inversion. 3.5 Choice of λ and β The regularization parameter λ is a trade-off between the data misfit and model roughness. For separate inversion, the L-curve method (Hansen 2000) has been frequently used to determine λ. In the joint inversion, the data weighting parameter β balances the relative weighting between SIP and NMR data sets. Similar parameters can also be found in other joint inversions (e.g. Gao et al.2012). When more controlling parameters are involved in the inversion, it is always better to search these two parameters simultaneously using the generalized L-curve method as proposed in Murat et al. (2002). However, this method is computationally costly, and in practice the selection of multiple controlling parameters is mainly based on trial-and-error (e.g. Gallardo & Meju 2004) or simple criteria (e.g. Commer & Newman 2009). In this study, we suggest a two-step method for the determination of the regularization parameter λ and data weighting parameter β. The method is illustrated in Fig. 3 using the synthetic example presented in Section 4. In this two-step method, the regularization parameter λ is first fixed, for example, as an intermediate value. Then, the joint inversion is carried out with β varying over a broad range, for example, from 10−3 to 103 in Fig. 3(a). By plotting the data misfit and β, a reasonable trade-off can then be chosen based on the shape of the curve. For instance, in Fig. 3(a) the optimal β is chosen to ensure the data misfit reaches a minimum. The second step is to determine the regularization parameter λ. Similarly, the joint inversion is conducted with the optimal β and a varying λ. The parameter λ and the related data misfit can be plotted to form a curve where the minimum can be chosen, for example, as shown in Fig. 3(b). We chose the associated λ as the optimal regularization parameter. Other methods such as the conventional L-curve method and its updated version (Hansen et al.2007) can also be used here to determine the optimal regularization parameter. Figure 3. View largeDownload slide The determination of the data weighting parameter β and regularization parameter λ in the synthetic example presented in Section 4. The circle symbols indicate the positions of the optimal β and λ. Figure 3. View largeDownload slide The determination of the data weighting parameter β and regularization parameter λ in the synthetic example presented in Section 4. The circle symbols indicate the positions of the optimal β and λ. 4 A SYNTHETIC SAMPLE In this section, a synthetic porous sample with a bimodal pore size distribution (see Fig. 4a) is used to demonstrate the applicability of the proposed joint inversion approach. The density of the pore size is modelled by a combination of two lognormal distributions. The related mean pore sizes are 20 μm and 2 μm; the standard deviations are 0.3 and 0.3; and the relative weighting is 0.4 and 0.6. Numerical experiments are conducted to produce the NMR responses [using eqs (1) and (4)] and SIP responses [using eqs (7) and (13)] of the sample. The NMR and SIP measurements are then jointly inverted to estimate the pore size distribution. Two cases were tested and are presented here. In Case 1, the numerical experiments and the inversion use the same petrophysical models. In Case 2, it is assumed we lack information on the petrophysical models. Thus, the inversion uses simplified petrophysical models, which are different from those used in the numerical experiments. Figure 4. View largeDownload slide The synthetic example in Case 1: (a) the pore size distribution, (b) the ‘measured’ NMR T2 relaxation, (c) the ‘measured’ real conductivity spectra and (d) the ‘measured’ imaginary conductivity spectra. Figure 4. View largeDownload slide The synthetic example in Case 1: (a) the pore size distribution, (b) the ‘measured’ NMR T2 relaxation, (c) the ‘measured’ real conductivity spectra and (d) the ‘measured’ imaginary conductivity spectra. 4.1 Case 1: inversion with exact petrophysical models In Case 1, the NMR T2 relaxation was numerically calculated using eqs (1) and (4) and the ‘measured’ decay curve is shown in Fig. 4(b). In the NMR experiment, the initial signal E0 is set as 0.5 μV; the surface relaxivity ρs is set as 10 μm s−1; the NMR signal was recorded for 5 s and the total number of the measurements are 3000. The random Gaussian errors were added to the simulated decay curve and the simulated noise level is 1 per cent. The complex conductivity of the sample is numerically calculated using eqs (7) and (13), and the ‘measured’ real and imaginary conductivities are shown in Figs 4(c) and (d), respectively. In SIP experiment, the following petrophysical parameters were used: σ0 = 0.045 S m−1, σ∞ = 0.05 S m−1, D+ = 1.32 × 10−9 m2 s−1, εinf = 100ε0 (where ε0 is the vacuum dielectric permittivity, 8.85 × 10−12 F m−1). The use of 100ε0 as the high frequency permittivity εinf is consistent with experimental measurements (e.g. see Jougnot et al.2010; Kremer et al.2016). The origin of this abnormally high value (for comparison, the permittivity of pure water is ∼79ε0) is still unclear and it could be caused by Maxwell-Wagner polarization, grain shape influence, inductive electromagnetic coupling or electrode influence (Kremer et al.2016). The applied minimum and maximum frequencies are 10−3 Hz and 100 kHz, respectively, and the total number of SIP measurements is 50. Similarly, random Gaussian errors are added to the calculation to simulate the measurement noise. The noise level is 1 per cent for real conductivity and 35 per cent for imaginary conductivity. Three individual calculations were conducted for each measured frequency. The mean value was used as the measured conductivity and the standard deviation was used as the measurement error as shown in Figs 4(c) and (d). The proposed joint inversion approach was used to estimate the pore size distribution of the sample based on the NMR and SIP measurements in Fig. 4. The proposed two-step method was used to determine the regularization parameter λ and data weighting parameter β. As shown in Figs 3(a) and (b), the reasonable trade-off was chosen when the data misfit reaches a minimum, and the determined optimal λ and β are 9771 and 153, respectively. In the joint inversion, the surface relaxivity was taken as the true value, that is, 10 μm s−1; and therefore, it was not updated during the inversion. The joint-inversion result is shown in Fig. 5(a). Separate NMR and SIP inversions were also carried out with the same regularization parameter, and the results are shown in Figs 5(b) and (c). The true density (normalized value) of the pore size is also plotted in Fig. 5 as well as discrepancies between the inverted and true densities. Figure 5. View largeDownload slide The inversion results of Case 1: (a) joint inversion, (b) separate NMR inversion, (c) separate SIP inversion and (d) discrepancies between inverted and true distributions of the pore size. Figure 5. View largeDownload slide The inversion results of Case 1: (a) joint inversion, (b) separate NMR inversion, (c) separate SIP inversion and (d) discrepancies between inverted and true distributions of the pore size. As shown in Fig. 5(a), the joint inversion successfully recovers the bimodal shape of the pore size distribution. The estimated positions of the two peaks are very close to the true locations although the related densities are slightly smaller than the true values. In contrast, the results from separate inversions have very low resolutions. The NMR inversion, to a certain extent, recovers the bimodal nature of the pore size but it does not correctly locate the dominant pore size; the SIP inversion renders a pore size distribution with only one peak. Comparing with separate inversions, the joint inversion significantly improves the resolution of the estimated pore size distribution. Clearly, this improvement is due to the addition of extra measurements in the inversion. We also tested whether the joint inversion can correctly estimate the surface relaxivity of the sample by allowing ρs to update during the inversion. The initial value of ρs was set as 1 μm s−1and other parameters were unchanged. The inverted pore size distribution is very close to the one presented in Fig. 5(a). The evolution of ρs during the inversion is plotted in Fig. 6. It shows that although the initial ρs is relatively small, after several iterations, it approaches to the true value. The final inverted surface relaxivity is 9.47 μm s−1, which is very close to the true value 10 μm s−1. Figure 6. View largeDownload slide The surface relaxivity evolution during the joint inversion in Case 1. Figure 6. View largeDownload slide The surface relaxivity evolution during the joint inversion in Case 1. 4.2 Case 2: inversion with inaccurate petrophysical models In Case 2, we use eqs (1) and (3) to simulate the NMR response of the sample so that the NMR relaxation is not necessarily within the fast diffusion regime. In SIP experiments, a Cole–Cole model is used to describe the unknown polarization that might be responsible for the anomalously high permittivity at intermediate frequencies (e.g. Lesmes & Morgan 2001; Kremer et al.2016). Thus, the complex conductivity is calculated by combing eqs (7) and (13) and the added Cole–Cole model. The NMR and SIP responses calculated using these ‘realistic’ petrophysical models are shown in Fig. 7. In the calculation, same model parameters were used here except that the surface relaxivity is changed to 100 μm s−1. For comparison, NMR and SIP responses calculated using the simplified petrophysical models (i.e. those used in Case 1) were also presented. As shown in Fig. 7, the ‘realistic’ and simplified petrophysical models give different NMR and SIP responses. The former relaxes the fast diffusion condition, and thus the related NMR relaxation rate is apparently slower (Fig. 7a). Moreover, ‘realistic’ petrophysical models give relatively large imaginary conductivities at frequencies larger 100 Hz (Fig. 7b), and the spectrum resembles experimental measurements shown in Lesmes & Morgan (2001). Figure 7. View largeDownload slide The measured (a) NMR and (b) SIP responses of the synthetic sample in Case 2. Compared with simplified petrophysical models, the realistic ones relax the fast diffusion condition and consider the unknown polarization responsible for anomalously high imaginary conductivity at intermediate frequencies. Figure 7. View largeDownload slide The measured (a) NMR and (b) SIP responses of the synthetic sample in Case 2. Compared with simplified petrophysical models, the realistic ones relax the fast diffusion condition and consider the unknown polarization responsible for anomalously high imaginary conductivity at intermediate frequencies. The NMR and SIP responses (simulated from ‘realistic’ petrophysical models) are inverted to estimate the pore size distribution jointly or separately. In the inversions, it is assumed that we do not know the exact petrophysical relations and the simplified petrophysical models [eqs (1), (4), (7), and (13)] are used. The estimated pore size distributions from joint and separate inversions are shown in Fig. 8. The joint inversion still recovers the bimodal nature of the pore size and two peak locations are also well defined. However, the estimated densities are significantly lower than the true values. In addition, the inversion gives a very broad range for the pore size distribution. Compared to Fig. 5(a), the joint inversion results in Fig. 8(a) differ notably from the true model. This is attributed to the use of the simplified petrophysical models in the inversion. Again, we compare the results of joint inversion to those from separate inversions. It appears the separate inversion, neither SIP nor NMR inversion, can capture the bimodal nature of the pore size. The joint inversion still improves the quality of the estimated pore size distribution despite that we lack enough information on the exact petrophysical models of the sample. Figure 8. View largeDownload slide The inversion results of Case 2: (a) joint inversion, (b) separate NMR inversion, (c) separate SIP inversion and (d) discrepancies between inverted and true distributions of the pore size. Figure 8. View largeDownload slide The inversion results of Case 2: (a) joint inversion, (b) separate NMR inversion, (c) separate SIP inversion and (d) discrepancies between inverted and true distributions of the pore size. 5 APPLICATION TO BEREA SANDSTONE In this section, we apply the proposed joint inversion approach to a Berea sandstone. The estimated pore size distribution is compared with those from separate inversions. 5.1 Sample and measurements The well-characterized Berea sandstone is chosen in this study because both NMR and SIP data are available. The sample has a porosity ϕ of 20.2 per cent and the specific surface area Sm is determined with nitrogen adsorption as 0.5 m2 g−1. Previous research shows quartz is the most abundant mineral in the sandstone and it also contains a significant amount of clays (e.g. Knight & Nur 1987). The intrinsic formation factor F of the sample is 16.3 as determined from multi-salinity measurement. The sample was saturated with NaCl solution (σw = 0.043 S m−1) for NMR and SIP measurement. The NMR T2 relaxation of the sample was measured using a 2 MHz Magritek Rock Core AnalyzerTM with CPMG pulse sequence and the echo spacing is 100 μs. The complex conductivity of the sample was measured in the frequency range between 10−3 Hz and 10 kHz using the SIP-LAB-II instrument (Radic). The NMR and SIP responses of the sample are shown in Fig. 9. Figure 9. View largeDownload slide The measured and fitted (a) NMR and (b) SIP responses (imaginary conductivity) of the Berea sandstone sample. The fitted curves are from joint and separate inversions. Figure 9. View largeDownload slide The measured and fitted (a) NMR and (b) SIP responses (imaginary conductivity) of the Berea sandstone sample. The fitted curves are from joint and separate inversions. 5.2 Inversion results The NMR and SIP data of the sample in Fig. 9 are jointly inverted to estimate the pore size distribution of Berea sandstone. The data weighting and regularization parameters β and λ are determined using the proposed two-step method as shown in Fig. 10. In the inversion, the diffusion coefficient of ions in the EDL D+ was taken as 3.2 × 10−12 m2 s−1 as suggested by Revil (2012), considering that the sample has a significant amount of clays; the value may need to be adjusted by trial-and-error in cases where the calculated and measured peak frequencies in SIP spectra do not coincide. The surface relaxivity is allowed to change in the inversion and the initial value is set as ρs = 5 μm s−1. However, we find ρs does not change during the inversion. The reason could be that the initial value was too close to the true surface relaxivity of the sample. The recovered petrophysical parameters from joint inversion are as follows: σ∞ = 0.0028 S m−1, σ0 = 0.0026 S m−1, εinf = 128ε0, E0 = 7 μV, ρs = 4.995 μm s−1. For comparison, separate inversions are also conducted. In NMR inversion, the surface relaxivity ρs is kept as 5 μm s−1; in SIP inversion, the diffusion coefficient D+ is taken as 3.2 × 10−12 m2 s−1. The regularization parameters used in separate inversions are individually determined using the traditional L-curve method so they may be different from that used in the joint inversion. The estimated pore size distributions from separate and joint inversions are shown in Fig. 11(a). The NMR and SIP responses calculated from the inversion results are also plotted in Fig. 9 for comparison. Figure 10. View largeDownload slide The determination of the data weighting parameter β and regularization parameter λ for the Berea sandstone sample. The circle symbols indicate the positions of the optimal β and λ. Figure 10. View largeDownload slide The determination of the data weighting parameter β and regularization parameter λ for the Berea sandstone sample. The circle symbols indicate the positions of the optimal β and λ. Figure 11. View largeDownload slide The estimated pore size distribution of the Berea sandstone using (a) SIP and NMR porosimetries and (b) mercury injection capillary pressure (MICP). The MICP-determined pore (throat) size distribution is from the Berea sandstone sample in Shi et al. (2011), which has a similar porosity (ϕ = 0.19) with the one used in this study (ϕ = 0.2). Figure 11. View largeDownload slide The estimated pore size distribution of the Berea sandstone using (a) SIP and NMR porosimetries and (b) mercury injection capillary pressure (MICP). The MICP-determined pore (throat) size distribution is from the Berea sandstone sample in Shi et al. (2011), which has a similar porosity (ϕ = 0.19) with the one used in this study (ϕ = 0.2). The microstructure of Berea sandstone has been extensively studied using different methods, including MICP (Shi et al.2011) and X-ray CT scan (Dong & Blunt 2009). Here, we use the sample in Shi et al. (2011) as an example to show the main features of pore (throat) size of Berea sandstone. The sample has a porosity of ∼19 per cent, which is very close to the sample in our study (∼20 per cent). The pore (throat) size distribution measured from MICP (Shi et al.2011) is plotted in Fig. 11(b). One profound feature of the sample is the narrow range of the pore size, which centres on ∼10 μm. In contrast, the pore size distribution estimated from separate SIP inversion is quite broad, covering more than 2 orders. Comparing to MICP results, the SIP inversion significantly overestimates the contribution of small pores (r < 3 μm) although the imaginary conductivity spectrum is well reproduced base on this broad pore size distribution (Fig. 9b). As shown in Fig. 11, the pore size estimated from joint inversion is consistent with the MICP results; both the shape and the relative position of the curve are well recovered, demonstrating the applicability of the joint inversion approach. The NMR response calculated from joint inversion fits the measurement perfectly. The simulated imaginary conductivity spectra agree with the measurements only at low (<0.1 Hz) and high frequencies (>10 kHz). At intermediate frequencies, the simulated σ″ is much smaller than the measurements. The discrepancy may be explained by the presence of other polarization, which is responsible for the measured quadrature conductivity at intermediate frequencies (Kremer et al.2016). The estimated pore size distribution from NMR inversion is very similar to that of joint inversion, indicating that more weights are put on NMR data during the joint inversion. That implies, for this sample, the NMR data contain more information on the pore size than the SIP measurements. 5.3 Discussion In this subsection, we briefly discuss limitations of the proposed joint inversion approach and possible directions for further studies. First, we point out that the Stern layer polarization model (i.e. eq. 5) may only be valid for relatively large pores (>∼0.1 μm) because it assumes the EDL thickness is much smaller than the radius of the curvature (see Schwarz 1962 or Appendix  A in this study). Therefore, for geomaterials characterized by small pore size (<∼0.1 μm) such as mudstones or clayey soils, the Stern layer polarization model may not be able to describe the measured SIP spectrum. In that case, the joint inversion of SIP and NMR data with our approach may not improve the estimation of pore size distribution if compared with separate NMR porosimetry. In addition, since both NMR and SIP porosimetries assume isolated pores, the joint inversion approach may not work properly for geomaterials with well-connected heterogeneous pores such as loosely packed soils. Despite the fact that the Stern layer polarization model has been frequently used in practice for various applications (e.g. Osterman et al. 2016; Leroy et al. 2017; Revil et al. 2017), it is still unclear whether it alone is sufficient to explain the measured complex conductivity of geomaterials over a broad frequency range. Thus, it is necessary to assess the applicability of the Stern layer polarization model for a given sample, and the proposed joint inversion approach is well suited for this purpose. Indeed, the simulated SIP response in Fig. 9(b) based on joint inversion of NMR and SIP data has indicated that Stern layer polarization may not explain the high imaginary conductivity of Berea sandstone at ∼100 Hz. In addition to EDL polarization, membrane polarization (Marshall & Madden 1959) is also commonly used to explain the induced polarization of porous media (e.g. Bairlein et al. 2016; Hördt et al. 2016). Different from Stern layer polarization, which relates the relaxation time to pore/particle size, membrane polarization suggests that the relaxation time may also be controlled by the geometry of the pore throat (Bücker & Hördt 2013). This conclusion was also supported by experimental evidence (e.g. Scott & Barker 2003). Considering that pore throat is usually smaller than pore space in geomaterials, membrane polarization could potentially explain the abnormal high imaginary conductivity at intermediate frequencies reported in Lesmes & Morgan (2001) and Kremer et al. (2016). In this perspective, the joint inversion incorporating membrane polarization would give additional information on the pore geometry of a geomaterial. Therefore, our current work might constitute a first step towards a comprehensive joint inversion that can extract the full pore geometry information from NMR and SIP data. 6 CONCLUSIONS In this study, we propose to jointly invert NMR and SIP measurements for the estimation of the pore size distribution of geomaterials. The joint inversion is based on the fact that both NMR T2 relaxation and low-frequency induced polarization of geomaterials are influenced by geometric characteristics of the pore space. In the joint inversion, the NMR T2 relaxation time and the SIP relaxation time τ are linked through their petrophysical models. The Tikhonov regularization is used to constrain the smoothness of the pore size distribution. The trade-offs between two data sets and the balance between data misfit and model regularization are determined using a two-step method. The inversion problem is solved iteratively using the Gauss–Newton method and the related sensitivity matrix at each iteration is determined analytically. Numerical examples have shown that the joint inversion could recover the main features of the pore size of the synthetic sample even if petrophysical models used in the inversion are inaccurate. Comparing to separate inversions, the addition of extra data set can improve the resolution of the estimated pore size distribution. In addition, the results from synthetic example indicate that the surface relaxivity could be properly estimated by the joint inversion when the diffusion coefficient of ions in the EDL is known. The proposed method has been applied to a Berea sandstone sample and the estimated pore size distribution is fairly consistent with that estimated from mercury injection capillary pressure measurements, demonstrating the applicability and usefulness of the proposed method. Acknowledgements We thank the editor Mark Everett for handling our submission, and Norbert Klitzsch along with an anonymous reviewer for their valuable review comments. We thank Manika Prasad at the Department of Petroleum Engineering of Colorado School of Mines for providing the Berea sandstone data sets. REFERENCES Aster R.C., Borchers B., Thurber C.H., 2013. Tikhonov regularization, in Parameter Estimation and Inverse Problems , 2nd edn, pp. 93– 127, Academic Press. Google Scholar CrossRef Search ADS   Bairlein K., Bücker M., Hördt A., Hinze B. 2016. Temperature dependence of spectral induced polarization data: experimental results and membrane polarization theory, Geophys. J. Int. , 205( 1), 440– 453. https://doi.org/10.1093/gji/ggw027 Google Scholar CrossRef Search ADS   Behroozmand A.A., Keating K., Auken E., 2015. A review of the principles and applications of the NMR technique for near-surface characterization, Surv. Geophys. , 36( 1), 27– 85. https://doi.org/10.1007/s10712-014-9304-0 Google Scholar CrossRef Search ADS   Binley A., Slater L.D., Fukes M., Cassiani G., 2005. Relationship between spectral induced polarization and hydraulic properties of saturated and unsaturated sandstone, Water Resour. Res. , 41, W12417, doi:10.1029/2005WR004202. https://doi.org/10.1029/2005WR004202 Borgia G.C., Bortolotti V., Brancolini A., Brown R.J.S., Fantazzini P., 1996. Developments in core analysis by NMR measurements, Magn. Reson. Imaging , 14( 7–8), 751– 760. https://doi.org/10.1016/S0730-725X(96)00160-9 Google Scholar CrossRef Search ADS PubMed  Borgia G.C., Brown R.J.S., Fantazzini P., 1998. Uniform-penalty inversion of multiexponential decay data, J. Magn. Reson. , 132( 1), 65– 77. https://doi.org/10.1006/jmre.1998.1387 Google Scholar CrossRef Search ADS PubMed  Brownstein K.R., Tarr C.E., 1979. Importance of classical diffusion in NMR studies of water in biological cells, Phys. Rev. A , 19( 6), 2446– 2453. https://doi.org/10.1103/PhysRevA.19.2446 Google Scholar CrossRef Search ADS   Bryar T.R., Daughney C.J., Knight R.J., 2000. Paramagnetic effects of iron(III) species on nuclear magnetic relaxation of fluid protons in porous media, J. Magn. Reson. , 142( 1), 74– 85. https://doi.org/10.1006/jmre.1999.1917 Google Scholar CrossRef Search ADS PubMed  Bücker M., Hördt A., 2013. Analytical modelling of membrane polarization with explicit parametrization of pore radii and the electrical double layer, Geophys. J. Int. , 194( 2), 804– 813. https://doi.org/10.1093/gji/ggt136 Google Scholar CrossRef Search ADS   Chelidze T.L., Gueguen Y., 1999. Electrical spectroscopy of porous rocks: a review—I. Theoretical models, Geophys. J. Int. , 137( 1), 1– 15. https://doi.org/10.1046/j.1365-246x.1999.00799.x Google Scholar CrossRef Search ADS   Commer M., Newman G.A., 2009. Three-dimensional controlled-source electromagnetic and magnetotelluric joint inversion, Geophys. J. Int. , 178 ( 3), 1305– 1316. https://doi.org/10.1111/j.1365-246X.2009.04216.x Google Scholar CrossRef Search ADS   Cotts R.M., Hoch M.J.R., Sun T., Markert J.T., 1989. Pulsed field gradient stimulated echo methods for improved NMR diffusion measurements in heterogeneous systems, J. Magn. Reson. (1969) , 83( 2), 252– 266. https://doi.org/10.1016/0022-2364(89)90189-3 Google Scholar CrossRef Search ADS   Dlubac K., Knight R., Keating K., 2014. A numerical study of the relationship between NMR relaxation and permeability in sands and gravels, Near Surf. Geophys. , 12( 2), 219– 230. https://doi.org/10.3997/1873-0604.2013042 Dong H., Blunt M.J., 2009. Pore-network extraction from micro-computerized-tomography images, Phys. Rev. E , 80( 3), doi:10.1103/PhysRevE.80.036307. https://doi.org/10.1103/PhysRevE.80.036307 D’Orazio F., Tarczon J.C., Halperin W.P., Eguchi K., Mizusaki T., 1989. Application of nuclear magnetic resonance pore structure analysis to porous silica glass, J. Appl. Phys. , 65( 2), 742– 751. https://doi.org/10.1063/1.343088 Google Scholar CrossRef Search ADS   Dullien F.A., 2012. Porous Media: Fluid Transport and Pore Structure , Academic press. Dullien F.A.L., 1981. Wood's metal porosimetry and its relation to mercury porosimetry, Powder Technol. , 29( 1), 109– 116. https://doi.org/10.1016/0032-5910(81)85009-7 Google Scholar CrossRef Search ADS   Florsch N., Revil A., Camerlynck C., 2014. Inversion of generalized relaxation time distributions with optimized damping parameter, J. Appl. Geophys. , 109, 119– 132. https://doi.org/10.1016/j.jappgeo.2014.07.013 Google Scholar CrossRef Search ADS   Foley I., Farooqui S.A., Kleinberg R.L., 1996. Effect of paramagnetic ions on NMR relaxation of fluids at solid surfaces, J. Magn. Reson. A , 123( 1), 95– 104. https://doi.org/10.1006/jmra.1996.0218 Google Scholar CrossRef Search ADS PubMed  Fredrich J.T., 1999. 3D imaging of porous media using laser scanning confocal microscopy with application to microscale transport processes, Phys. Chem. Earth, Part A: Solid Earth Geod. , 24( 7), 551– 561. https://doi.org/10.1016/S1464-1895(99)00079-4 Google Scholar CrossRef Search ADS   Gallardo L.A., Meju M.A., 2003. Characterization of heterogeneous near-surface materials by joint 2D inversion of dc resistivity and seismic data, Geophys. Res. Lett. , 30( 13), 1658, doi:10.1029/2003GL017370. https://doi.org/10.1029/2003GL017370 Gallardo L.A., Meju M.A., 2004. Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints, J. geophys. Res. , 109, B03311, doi:10.1029/2003JB002716. https://doi.org/10.1029/2003JB002716 Gao G., Abubakar A., Habashy T.M., 2012. Joint petrophysical inversion of electromagnetic and full-waveform seismic data, Geophysics , 77( 3), WA3– WA18. https://doi.org/10.1190/geo2011-0157.1 Google Scholar CrossRef Search ADS   Grunewald E., Knight R., 2009. A laboratory study of NMR relaxation times and pore coupling in heterogeneous media, Geophysics , 74( 6), E215– E221. https://doi.org/10.1190/1.3223712 Google Scholar CrossRef Search ADS   Grunewald E., Knight R., 2011. The effect of pore size and magnetic susceptibility on the surface NMR relaxation parameter T2*, Near Surf. Geophys. , 9( 2), 169– 178. Haber E., Oldenburg D., 1997. Joint inversion: a structural approach, Inverse Probl. , 13( 1), 63– 77. https://doi.org/10.1088/0266-5611/13/1/006 Google Scholar CrossRef Search ADS   Hanai T., 1960. Theory of the dielectric dispersion due to the interfacial polarization and its application to emulsions, Kolloid-Zeitschrift , 171( 1), 23– 31. https://doi.org/10.1007/BF01520320 Google Scholar CrossRef Search ADS   Hansen P.C., 2000. The L-curve and its use in the numerical treatment of inverse problems, in Computational Inverse Problems in Electrocardiology , pp. 119– 142, ed. Johnston P., WIT Press. Hansen P.C., Jensen T.K., Rodriguez G., 2007. An adaptive pruning algorithm for the discrete L-curve criterion, J. Comput. Appl. Math. , 198( 2), 483– 492. https://doi.org/10.1016/j.cam.2005.09.026 Google Scholar CrossRef Search ADS   Hördt A., Bairlein K., Bielefeld A., Bücker M., Kuhn E., Nordsiek S., Stebner H. 2016. The dependence of induced polarization on fluid salinity and pH, studied with an extended model of membrane polarization, J. Appl. Geophys. , 135, 408– 417. https://doi.org/10.1016/j.jappgeo.2016.02.007 Google Scholar CrossRef Search ADS   Jardani A., Revil A., 2009. Stochastic joint inversion of temperature and self-potential data, Geophys. J. Int. , 179( 1), 640– 654. https://doi.org/10.1111/j.1365-246X.2009.04295.x Google Scholar CrossRef Search ADS   Johnson D.L., Koplik J., Schwartz L.M., 1986. New pore-size parameter characterizing transport in porous media, Phys. Rev. Lett. , 57( 20), 2564– 2567. https://doi.org/10.1103/PhysRevLett.57.2564 Google Scholar CrossRef Search ADS PubMed  Jougnot D., Ghorbani A., Revil A., Leroy P., Cosenza P. 2010. Spectral induced polarization of partially saturated clay-rocks: a mechanistic approach, Geophys. J. Int. , 180( 1), 210– 224. https://doi.org/10.1111/j.1365-246X.2009.04426.x Google Scholar CrossRef Search ADS   Katz A.J., Thompson A.H., 1986. Quantitative prediction of permeability in porous rock, Phys. Rev. B , 34( 11), 8179– 8181. https://doi.org/10.1103/PhysRevB.34.8179 Google Scholar CrossRef Search ADS   Keating K., Falzone S., 2013. Relating nuclear magnetic resonance relaxation time distributions to void-size distributions for unconsolidated sand packs, Geophysics , 78( 6), D461– D472. https://doi.org/10.1190/geo2012-0461.1 Google Scholar CrossRef Search ADS   Kemna A. et al.  , 2012. An overview of the spectral induced polarization method for near-surface applications, Near Surf. Geophys. , 10( 6), 453– 468. https://doi.org/10.3997/1873-0604.2012027 Kleinberg R.L., Horsfield M.A., 1990. Transverse relaxation processes in porous sedimentary rock, J. Magn. Reson. (1969) , 88( 1), 9– 19. https://doi.org/10.1016/0022-2364(90)90104-H Google Scholar CrossRef Search ADS   Knight R.J., Nur A., 1987. The dielectric constant of sandstones, 60 kHz to 4 MHz, Geophysics , 52( 5), 644– 654. https://doi.org/10.1190/1.1442332 Google Scholar CrossRef Search ADS   Kosugi K., 1996. Lognormal distribution model for unsaturated soil hydraulic properties, Water Resour. Res. , 32( 9), 2697– 2703. https://doi.org/10.1029/96WR01776 Google Scholar CrossRef Search ADS   Kosugi K.I., 1994. Three-parameter lognormal distribution model for soil water retention, Water Resour. Res. , 30( 4), 891– 901. https://doi.org/10.1029/93WR02931 Google Scholar CrossRef Search ADS   Kremer T., Schmutz M., Leroy P., Agrinier P., Maineult A., 2016. Modelling the spectral induced polarization response of water-saturated sands in the intermediate frequency range (102–105 Hz) using mechanistic and empirical approaches, Geophys. J. Int. , 207( 2), 1303– 1312. https://doi.org/10.1093/gji/ggw334 Google Scholar CrossRef Search ADS   Kruschwitz S., Binley A., Lesmes D., Elshenawy A., 2010. Textural controls on low-frequency electrical spectra of porous media, Geophysics , 75( 4), WA113– WA123. https://doi.org/10.1190/1.3479835 Google Scholar CrossRef Search ADS   Kuila U., Prasad M., 2013. Specific surface area and pore-size distribution in clays and shales, Geophys. Prospect. , 61( 2), 341– 362. https://doi.org/10.1111/1365-2478.12028 Google Scholar CrossRef Search ADS   Leroy P., Revil A., Kemna A., Cosenza P., Ghorbani A., 2008. Complex conductivity of water-saturated packs of glass beads, J. Colloid Interface Sci. , 321( 1), 103– 117. https://doi.org/10.1016/j.jcis.2007.12.031 Google Scholar CrossRef Search ADS PubMed  Leroy P., Li S., Jougnot D., Revil A., Wu Y., 2017. Modelling the evolution of complex conductivity during calcite precipitation on glass beads. Geophys. J. Int. , 209( 1), 123– 140. Lesmes D.P., Morgan F.D., 2001. Dielectric spectroscopy of sedimentary rocks, J. geophys. Res. , 106( B7), 13 329–13 346. https://doi.org/10.1029/2000JB900402 Google Scholar CrossRef Search ADS   Linde N., Doetsch J., 2016. Joint inversion in hydrogeophysics and near-surface geophysics, in Integrated Imaging of the Earth , ch. 7, pp. 119– 135, eds Moorkamp M., Lelievre P., Linde N., Khan A., John Wiley & Sons, Inc. Google Scholar CrossRef Search ADS   Lyklema J., Dukhin S., Shilov V., 1983. The relaxation of the double layer around colloidal particles and the low-frequency dielectric dispersion, J. Electroanal. Chem. Interfacial Electrochem. , 143( 1–2), 1– 21. https://doi.org/10.1016/S0022-0728(83)80251-4 Google Scholar CrossRef Search ADS   Marshall D.J., Madden T.R. 1959. Induced polarization, a study of its causes, Geophysics , 24( 4), 790– 816. https://doi.org/10.1190/1.1438659 Google Scholar CrossRef Search ADS   McBride M., 1989. Surface chemistry of soil minerals, in Minerals in Soil Environments , pp. 3588, eds Dixon, J.B. & Weed, S.B., Soil Science Society of America. McGillivray P.R., Oldenburg D.W., 1990. Methods for calculating Frechet derivatives and sensitivities for the non-linear inverse problem: a comparative study, Geophys. Prospect. , 38, 499– 524. https://doi.org/10.1111/j.1365-2478.1990.tb01859.x Google Scholar CrossRef Search ADS   Mohnke O., 2014. Jointly deriving NMR surface relaxivity and pore size distributions by NMR relaxation experiments on partially desaturated rocks, Water Resour. Res. , 50, 5309– 5321. https://doi.org/10.1002/2014WR015282 Google Scholar CrossRef Search ADS   Mohnke O., Yaramanci U., 2008. Pore size distributions and hydraulic conductivities of rocks derived from Magnetic Resonance Sounding relaxation data using multi-exponential decay time inversion, J. Appl. Geophys. , 66, 73– 81. https://doi.org/10.1016/j.jappgeo.2008.05.002 Google Scholar CrossRef Search ADS   Moorkamp M., Heincke B., Jegen M., Roberts A.W., Hobbs R.W., 2011. A framework for 3-D joint inversion of MT, gravity and seismic refraction data, Geophys. J. Int. , 184, 477– 493. https://doi.org/10.1111/j.1365-246X.2010.04856.x Google Scholar CrossRef Search ADS   Müller-Petke M., Dlugosch R., Lehmann-Horn J., Ronczka M., 2015. Nuclear magnetic resonance average pore-size estimations outside the fast-diffusion regime, Geophysics , 80, D195– D206. https://doi.org/10.1190/geo2014-0167.1 Google Scholar CrossRef Search ADS   Murat B., Misha E.K., Eric L.M., 2002. Efficient determination of multiple regularization parameters in a generalized L-curve framework, Inverse Probl. , 18, 1161– 1183. https://doi.org/10.1088/0266-5611/18/4/314 Google Scholar CrossRef Search ADS   Niu Q., Revil A., 2016. Connecting complex conductivity spectra to mercury porosimetry of sedimentary rocks, Geophysics , 81, E17– E32. https://doi.org/10.1190/geo2015-0072.1 Google Scholar CrossRef Search ADS   Niu Q., Prasad M., Revil A., Saidian M., 2016. Textural control on the quadrature conductivity of porous media, Geophysics , 81, E297– E309. https://doi.org/10.1190/geo2015-0715.1 Google Scholar CrossRef Search ADS   Osterman G., Keating K., Binley A., Slater L., 2016. A laboratory study to estimate pore geometric parameters of sandstones using complex conductivity and nuclear magnetic resonance for permeability prediction, Water Resour. Res. , 52, 4321– 4337. https://doi.org/10.1002/2015WR018472 Google Scholar CrossRef Search ADS   Park S.K., Van G.P., 1991. Inversion of pole–pole data for 3-D resistivity structure beneath arrays of electrodes, Geophysics , 56, 951– 960. https://doi.org/10.1190/1.1443128 Google Scholar CrossRef Search ADS   Perrier E., Rieu M., Sposito G., de Marsily G., 1996. Models of the water retention curve for soils with a fractal pore size distribution, Water Resour. Res. , 32, 3025– 3031. https://doi.org/10.1029/96WR01779 Google Scholar CrossRef Search ADS   Pidlisecky A., Haber E., Knight R., 2007. RESINVM3D: a 3D resistivity inversion package, Geophysics , 72, H1– H10. https://doi.org/10.1190/1.2402499 Google Scholar CrossRef Search ADS   Revil A., 2012. Spectral induced polarization of shaly sands: influence of the electrical double layer, Water Resour. Res. , 48, W02517, doi:10.1029/2011WR011260. https://doi.org/10.1029/2011WR011260 Revil A., 2013. Effective conductivity and permittivity of unsaturated porous materials in the frequency range 1 mHz–1 GHz, Water Resour. Res. , 49, 306– 327. https://doi.org/10.1029/2012WR012700 Google Scholar CrossRef Search ADS PubMed  Revil A., Glover P.W.J., 1998. Nature of surface electrical conductivity in natural sands, sandstones, and clays, Geophys. Res. Lett. , 25, 691– 694. https://doi.org/10.1029/98GL00296 Google Scholar CrossRef Search ADS   Revil A., Koch K., Holliger K., 2012. Is it the grain size or the characteristic pore size that controls the induced polarization relaxation time of clean sands and sandstones?, Water Resour. Res. , 48, W05602, doi:10.1029/2011WR011561. https://doi.org/10.1029/2011WR011561 Revil A., Florsch N., Camerlynck C., 2014. Spectral induced polarization porosimetry, Geophys. J. Int. , 198, 1016– 1033. https://doi.org/10.1093/gji/ggu180 Google Scholar CrossRef Search ADS   Revil A., Binley A., Mejus L., Kessouri P., 2015. Predicting permeability from the characteristic relaxation time and intrinsic formation factor of complex conductivity spectra, Water Resour. Res. , 51, 6672– 6700. https://doi.org/10.1002/2015WR017074 Google Scholar CrossRef Search ADS   Revil A. et al.  , 2017. Complex conductivity of soils, Water Resour. Res. , 53, 7121– 7147. https://doi.org/10.1002/2017WR020655 Google Scholar CrossRef Search ADS   Schwarz G., 1962. A theory of the low-frequency dielectric dispersion of colloidal particles in electrolyte solution, J. Phys. Chem. , 66, 2636– 2642. https://doi.org/10.1021/j100818a067 Google Scholar CrossRef Search ADS   Scott J.B.T., Barker R.D., 2003. Determining pore-throat size in Permo-Triassic sandstones from low-frequency electrical spectroscopy, Geophys. Res. Lett. , 30(9), 1450, doi:10.1029/2003GL016951. https://doi.org/10.1029/2003GL016951 Shi J.-Q., Xue Z., Durucan S., 2011. Supercritical CO2 core flooding and imbibition in Berea sandstone—CT imaging and numerical simulation, Energy Procedia , 4, 5001– 5008. https://doi.org/10.1016/j.egypro.2011.02.471 Google Scholar CrossRef Search ADS   Tarantola A., Valette B., 1982. Generalized nonlinear inverse problems solved using the least squares criterion, Rev. Geophys. , 20, 219– 232. https://doi.org/10.1029/RG020i002p00219 Google Scholar CrossRef Search ADS   Timur A., Hempkins W.B., Weinbrandt R.M., 1971. Scanning electron microscope study of pore systems in rocks, J. geophys. Res. , 76, 4932– 4948. https://doi.org/10.1029/JB076i020p04932 Google Scholar CrossRef Search ADS   Titov K., Tarasov A., Ilyin Y., Seleznev N., Boyd A., 2010. Relationships between induced polarization relaxation time and hydraulic properties of sandstone, Geophys. J. Int. , 180, 1095– 1106. https://doi.org/10.1111/j.1365-246X.2009.04465.x Google Scholar CrossRef Search ADS   Titterington D., 1985. General structure of regularization procedures in image reconstruction, Astron. Astrophys. , 144, 381– 387. Tong M., Li L., Wang W., Jiang Y., 2006. Determining capillary-pressure curve, pore-size distribution, and permeability from induced polarization of shaley sand, Geophysics , 71, N33– N40. https://doi.org/10.1190/1.2195989 Google Scholar CrossRef Search ADS   Tuli A., Kosugi K., Hopmans J.W., 2001. Simultaneous scaling of soil water retention and unsaturated hydraulic conductivity functions assuming lognormal pore-size distribution, Adv. Water Resour. , 24, 677– 688. https://doi.org/10.1016/S0309-1708(00)00070-1 Google Scholar CrossRef Search ADS   Vervoort R., Cattle S., 2003. Linking hydraulic conductivity and tortuosity parameters to pore space geometry and pore-size distribution, J. Hydrol. , 272, 36– 49. https://doi.org/10.1016/S0022-1694(02)00253-6 Google Scholar CrossRef Search ADS   Volkmann J., Klitzsch N., 2015. Wideband impedance spectroscopy from 1mHz to 10 MHz by combination of four- and two-electrode methods, J. Appl. Geophys. , 114, 191– 201. https://doi.org/10.1016/j.jappgeo.2015.01.012 Google Scholar CrossRef Search ADS   Volkmann J., Klitzsch N., 2016. Evaluation of low frequency polarization models using well characterized sintered porous glass samples, J. Appl. Geophys. , 124, 39– 53. https://doi.org/10.1016/j.jappgeo.2015.11.011 Google Scholar CrossRef Search ADS   Weigand M., Kemna A., 2016. Debye decomposition of time-lapse spectral induced polarisation data, Comput. Geosci. , 86, 34– 45. https://doi.org/10.1016/j.cageo.2015.09.021 Google Scholar CrossRef Search ADS   Wildenschild D., Sheppard A.P., 2013. X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems, Adv. Water Resour. , 51, 217– 246. https://doi.org/10.1016/j.advwatres.2012.07.018 Google Scholar CrossRef Search ADS   Zhuravlev L.T., 2000. The surface chemistry of amorphous silica. Zhuravlev model, Colloids Surf. A , 173, 1– 38. https://doi.org/10.1016/S0927-7757(00)00556-2 Google Scholar CrossRef Search ADS   APPENDIX A: RELAXATION TIME OF EDL POLARIZATION Consider a spherical surface Γ separating two domains Ωi and Ωe (Fig. A1). The charge density of the surface is QS, and the radius of the curvature of the surface is rc. Apply an oscillating electric field E with its magnitude in the form of E = E0eiωt where ω is the angular frequency. Define spherical coordinates r and θ with r = 0 at the centre of domain Ωi and θ = 0 in the direction of E. Figure A1. View largeDownload slide A spherical surface Γ separating two domains Ωi and Ωe. The radius of the curvature of the surface is rc. Figure A1. View largeDownload slide A spherical surface Γ separating two domains Ωi and Ωe. The radius of the curvature of the surface is rc. The potential distribution in the domains φi and φe follows Laplace equation,   $${\nabla ^2}\ {\varphi ^i} = \ 0 \left( {r < {r_c}} \right)$$ (A1)and   $${\nabla ^2}\ {\varphi ^e} = \ 0 \left( {r > {r_c}} \right).$$ (A2) Define the potential at the surface Γ as φS. At least, eqs (A1) and (A2) satisfy the following boundary conditions: φe → −Ercos θ as r → ∞, φi is finite as r → 0, and $${\varphi _S} = \mathop {\lim }\limits_{r = {r_c}} {\varphi ^e}{\rm{\ }} = \mathop {\lim }\limits_{r = {r_c}} {\varphi ^i}{\rm{\ }}$$. The solutions to eqs (A1) and (A2) are well known and can be expressed as   $${\varphi ^i} = \sum \left( {A_n^i{r^n} + B_n^i\frac{1}{{{r^{n + 1}}}}} \right)\ {P_n}\left( {\cos \theta } \right)$$ (A3)and   $${\varphi ^e} = \sum \left( {A_n^e{r^n} + B_n^e\frac{1}{{{r^{n + 1}}}}} \right)\ {P_n}\left( {\cos \theta } \right)$$ (A4)where $$A_n^i$$, $$B_n^i$$, $$A_n^e$$, $$B_n^e$$ are unknowns in the coefficients for the nth term of Legendre polynomials Pn(cos θ). Applying boundary condition (1), eq. (A4) can be simplified as   $${\varphi ^e} = \left( {A_1^er + B_1^e\frac{1}{{{r^2}}}} \right)\ \cos \theta .$$ (A5)Then, the potential at the surface can be obtained as   $${\varphi _S} = \left( {A_1^e{r_c} + B_1^e\frac{1}{{{r_c}^2}}} \right)\ \cos \theta .$$ (A6) If we expand φS using Legendre polynomials, φS has the following form:   $${\varphi _S} = \sum {\beta _n}{P_n}\left( {\cos \theta } \right),$$ (A7)where βn is the associated coefficient. Comparing eqs (A6) and (A7), we find that only the first term (n = 1) in eq. (A7) is non-zero. This result will be used later to determine the relaxation time of the change in QS in the surface Γ. Assume the charge on Γ (a thin EDL) only moves tangentially. The movement of charge on Γ is due to two driving forces: (1) electric field and (2) ion concentration gradient. Considering the continuity, the resulting time derivative of the surface charge density can be written as   $$\frac{{\partial {Q_S}}}{{\partial t}} = {e_0}{\rm{\ }}\nabla \cdot \left( {{j_E} + {j_D}} \right),$$ (A8)where e0 is the charge carried by an ion in the EDL, and jE and jD are surface fluxes of ions induced by electric field and ion diffusion, respectively. Consider the change in QS is ΔQS, which is much less than QS. Also, consider the expressions of jE and jD. Then, eq. (A8) can be expressed as   \begin{eqnarray}i\omega \ \Delta {Q_S} \!=\! \frac{{\mu kT}}{{{r_c}^2}}\ \frac{1}{{\sin \theta }}\frac{{\rm{d}}}{{{\rm{d}}\theta }}\left( {\sin \theta \frac{{{\rm{d}}\Delta {Q_S}}}{{{\rm{d}}\theta }} \!+\! \frac{{{e_0}{Q_S}}}{{kT}}\sin \theta \frac{{{\rm{d}}{\varphi _S}}}{{{\rm{d}}\theta }}} \right),\end{eqnarray} (A9)where μ is the mobility of ions in EDL, k is Boltzmann constant and T is the absolute temperature. We also expand the surface charge variation ΔQS using Legendre polynomials as   $$\Delta {Q_S} = \sum {\alpha _n}{P_n}\left( {\cos \theta } \right).$$ (A10) Insert eqs (A7) and (A10) into eq. (A9), and after some derivative calculations, eq. (A9) becomes   \begin{eqnarray}i\omega \ \sum {\alpha _n}{P_n} &=& \frac{{\mu kT}}{{{r_c}^2}}\ \sum {\alpha _n}\left( {\frac{{\cos \theta }}{{\sin \theta }}P_n^{\prime} + P_n^{\prime\prime}} \right) \nonumber\\ &&+ \frac{{\mu {e_0}{Q_S}}}{{{r_c}^2}}\sum {\beta _n}\left( {\frac{{\cos \theta }}{{\sin \theta }}P_n^{\prime} + P_n^{\prime\prime}} \right).\end{eqnarray} (A11) Considering the following equality,   $$P_n^{\prime\prime} + \frac{{\cos \theta }}{{\sin \theta }}P_n^\prime + n\left( {n + 1} \right)\ {P_n} = \ 0,$$ (A12)eq. (A11) can be rewritten as   \begin{eqnarray}\sum {\alpha _n}\left( {i\omega \!+\! \frac{{\mu kT}}{{{r_c}^2}}n\left( {n + 1} \right)} \right){P_n} \!=\! \ \sum {\beta _n}\frac{{\mu {e_0}{Q_S}}}{{{r_c}^2}}n\left( {n + 1} \right){P_n}.\nonumber\\ \end{eqnarray} (A13) Because of the orthogonality of Legendre polynomials, the following relation is valid for all the terms,   $${\alpha _n} = \ - \frac{1}{{1 + i\omega \frac{{{r_c}^2}}{{n\left( {n + 1} \right)\mu kT}}}}\frac{{{e_0}{Q_S}}}{{kT}}{\beta _n}.$$ (A14) As discussed, the coefficient βnfor the surface potential φS is non-zero only when n = 1, we have the following relation between ΔQS and φS:   $$\Delta {Q_S} = \ - \frac{1}{{1 + i\omega \tau }}\frac{{{e_0}{Q_S}}}{{kT}}{\varphi _S}$$ (A15) with   $$\tau \ = \frac{{{r_c}^2}}{{2\mu kT}}\ = \frac{{{r_c}^2}}{{2{D^ + }}}\$$ (A16)where D+ is the diffusion coefficient of the ions in the EDL. Eq. (A16) states that the surface charge density follows a simple Debye relaxation with relaxation time τ, which is controlled by the radius of the curvature of the surface Γ. In the derivation, the electrical properties of domains Ωi and Ωe were not specified, and therefore, eq. (A15) holds for both following cases: (1) domain Ωi is solid particle and domain Ωe is electrolyte; and (2) domain Ωi is electrolyte and domain Ωe is solid particle. In the former case, the radius of the curvature of the surface Γ is equal to the particle radius; for the latter, rc is equal to the pore size. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations