A Bayesian Optimization Approach to the Simultaneous Extraction of Intrinsic Physical Parameters from T₁ and T₂ Relaxation Responses

NMR relaxation and diffusion measurements as probe of the microstructure of porous media have a long history. Focusing on relaxation responses, this includes the estimation of pore size distributions (Prammer 1994; Kenyon 1997), (pseudo) capillary pressure curves (Marschall et al. 1995; Slijkerman et al. 2001), wettability index (Fleury and Deflandre 2003; Looyestijn and Hofman 2005; Chen et al. 2006; Al-Garadi et al. 2022), and permeability (Timur 1968; Straley et al. 1987; Banavar and Schwartz 1987; Sen et al. 1990; Borgia et al. 1998; Dunn et al. 1999). Length scale interpretations of NMR relaxation responses from saturated rocks are based on the assumption of the fast diffusion limit (Brownstein and Tarr 1979) (i.e., there is a direct correspondence between a specific relaxation time component and pore size). The assumption is valid if a range of conditions are fulfilled, including fast diffusion and diffusionally isolated pores, insignificant internal gradients (absence of paramagnetic impurities), and uniform mineralogy. Under those conditions, the relationship between relaxation time and pore size may be described by a mono-exponential decay with relaxation rate (Torrey 1956)

where T₂ and $T2b$ are transverse relaxation time and bulk transverse relaxation time, respectively; $ρ2$ is transverse surface relaxivity; $S$ and $V$ are surface area and volume of an individual pore, respectively; $De$ is effective diffusion coefficient; $γ$ is the gyromagnetic ratio; $G$ is the average or effective internal gradient; and $tE$ is echo time. Summing over all pores results in the magnetization decay of the sample. Assuming a pore shape (e.g., spherical, and constant surface relaxivity), the magnetization decay can then be converted to a pore size. Practically, the prediction of pore sizes requires a calibration of field data to NMR and mercury intrusion porosimetry laboratory measurements using surface relaxivity as the constant conversion parameter. In literature, such a parameter is often (and correctly) called an effective surface relaxivity (Chen et al. 2006) or apparent surface relaxivity (Müller-Petke et al. 2015; Costabel et al. 2018) to emphasize the averaged nature of the parameter (Marschall et al. 1995). Depending on the method utilized for estimation of the surface relaxivity (e.g., using surface area measurements by gas adsorption analysis, accounting for surface roughness or not, utilizing short time limit restricted diffusion, by Bayesian magnetic resonance derived pore size distribution or derived by utilizing nonground eigenstates following Brownstein-Tarr theory), the literature-reported surface relaxivity values vary by an order of magnitude. For the benchmark rock we use in this work, Bentheimer sandstone, the reported values of $ρ2$ spread from 8 μm/s (Holland et al. 2013) to 26 μm/s (Afrough et al. 2019).

The application of the concept of constant effective surface relaxivity to realistic heterogeneous rocks leads to errors of various magnitude. Toumelin et al. (2003b) demonstrated possible errors approaching 50% in estimating immobile fluid fraction due to unaccounted diffusional coupling, Arns et al. (2006) reported a 30% spread of permeability estimates due to relaxivity heterogeneity, while Benavides et al. (2017) showed a very significant narrowing of the predicted pore size distribution (twice as narrow compared to a known reference on a standard log scale) if a uniform relaxivity value is utilized. Approaching this problem by introducing computational methods utilizing intrinsic phase-specific surface properties may allow one to quantify those uncertainties and in general enhance the interpretation of NMR relaxation responses.

Additional information about fluid content or compartmentalization of the pore space can be gained from T₁-T₂ correlation experiments (Song et al. 2002; Tinni et al. 2014). Pseudo T₁-T₂ correlations may be derived by utilizing independent T₁ and T₂ measurements (Williamson et al. 2018). As mentioned above, it is well known that NMR relaxation of fluids in porous media is very sensitive to the wetting properties of solid phases. The close relationship between affinity and surface relaxivity was utilized to develop an NMR wettability index under the assumption of uniform wetting properties. The concept was extended for the more practical case of heterogeneous wettability of rock [e.g., we demonstrated in Shikhov et al. (2018) the approach to determine fractional wetting properties of grains surfaces using a combination of transverse relaxation techniques, though the reference data for the rock of interest at two end wetting states in either case is required, limiting applications to laboratory core analysis]. Approaches relying on both types of relaxation (T₁ and T₂) may not require such calibration, allowing estimates of arbitrary wet rocks at field conditions (e.g., well logs). Korb et al. (2018) utilized a combination of 2D T₁-T₂ and T₁ field-dependent relaxation (dispersion) measurements to probe the wetting properties of fluids confined in the pore space of rocks without prior knowledge of endpoint wettability state responses. Gizatullin et al. (2019) practically demonstrated a high sensitivity of certain types of aromatic compounds to wetting conditions of sandstones and carbonates for T₁ dispersion measurements. This suggests that future NMR methods of formation wettability characterization from well logs are likely to be based on a combination of T₂ and T₁ field-dependent measurements. The interpretation of such measurements would require an ability to separate the responses from internal gradient effects when present and efficient extraction of multiple longitudinal and transverse surface relaxivity values $ρ1$ and $ρ2$⁠. The current work presents a significant step toward this goal.

Numerical simulations based on digitized materials offer a way to test such methods and potentially overcome some of their limitations. Basic random walk methods on digitized media were introduced in several studies (Mendelson 1990; Wilkinson et al. 1991; Bergman et al. 1995). Straley and Schwartz (1996) compared experiments on glass beads to numerical simulations with reasonable accuracy. Hidajat et al. (2002) introduced a first-passage random walk method for NMR responses and correlated T₂ log mean relaxation time to lattice Boltzmann calculations of permeability for synthetic random structures. Øren et al. (2002) modeled the NMR T₂ response of Fontainebleau sandstone utilizing a random walk method with fixed timestep. Toumelin et al. (2003b) considered a periodic array of spheres to capture the response of vuggy carbonates at different water saturations conceptually, where the microporosity itself also formed a periodic array of spheres. More algorithmic detail was given in Toumelin et al. (2003a)—the approach captures dephasing in the presence of an external constant background gradient and utilizes the first-passage algorithm. Valckenborg et al. (2003) extended the modeling of NMR responses to nonuniform fields. Song et al. (2000) measured the magnetization decay due to diffusion in the internal field, a technique that probes higher eigenmodes of the NMR T₂ relaxation response. It was first modeled by Arns (2004) on a series of micro-computed tomography (CT) images. Interestingly, the associated length scale estimates indicated a larger pore size as well as a wider pore size distribution than would have resulted from the direct use of the ground states. Arns et al. (2005) compared a range of length scale measures for permeability estimation including NMR-derived ones, showing strong correlations between NMR-based estimates and direct numerical simulation for Fontainebleau sandstone, in agreement with Sen et al. (1990). Numerical techniques operating on 3D representations of microstructure extracted from micro-CT images have since been refined further and were applied to higher-dimensional NMR measurements (e.g., $D$-T₂) (Arns et al. 2011). Internal field effects were accounted for by application of a dipole approximation with random walkers accumulating an effective phase; as in Arns (2004), a part of the dephasing behavior is thus accounted for explicitly, rather than by an effective relaxation rate.

Given the significant advances in numerical techniques, a standing problem is the determination of the correct intrinsic parameters associated with relaxation at surfaces or smaller-scale microstructure where bulk relaxation rates of effective materials may be specified. However, determination for those parameters using forward simulations requires an extensive search in the multidimensional solution space, posing an ill-conditioned inverse problem, the solution of which may lead to new NMR relaxation response interpretation techniques directly utilizing intrinsic physical properties. This is particularly so if both T₁ and T₂ responses are considered. There are two major issues related to the relaxation time inverse problem: The micro-CT image-based NMR simulation is expensive to evaluate and the solution sets are nonunique. For the computational burden, Bayesian optimization has become a powerful strategy for identifying the extrema of the expensive objective functions with a limited budget (Mockus 1994; Jones et al. 1998; Snoek et al. 2012), and we have developed ST-ISW demonstrating the efficiency and reproducibility of Bayesian optimization for T₂ relaxation time inverse problem in which multiple solution sets are identified simultaneously (Li et al. 2021). The key advantage of Bayesian optimization is that it models the objective function using GPs or kriging known in geostatistics as the cheap proxy, with which the most promising candidate is identified and proposed for exact function evaluation so that the computational cost is significantly reduced. GPs as a surrogate or proxy model have been applied to solving history matching problems (Horowitz et al. 2013; Sayyafzadeh 2015; Rana et al. 2018; Dachanuwattana et al. 2019; Li et al. 2015) and to recovering physical, geometrical, or structural parameters where expensive measurements or simulations are involved (Vargas-Hernández et al. 2019; Deng et al. 2020). For the identification of nonunique solution sets under box constraints, the objective function is intrinsically multimodal. Li et al. (2021) recovered both the global optimal solution and (major) local optimal solutions using a multimodal search strategy, comprising a multistart Box-constrained limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm (L-BFGS-B) optimizer searching for local optimum solutions and a global optimizer social-learning particle swarm optimizer for global optimum solutions. With this technique, more than one physically valid solution set is identified by computationally reproducing the transverse relaxation time (T₂) measurement using the NMR simulation. From the mathematical point of view, the solutions are considered more physically valid if the agreement between the simulated longitudinal relaxation time (T₁) distribution arising from the identified physical quantities and the measured T₁ distribution can also be achieved. From an NMR forward solver point of view, it is easy to output T₁ and T₂ relaxation responses with little extra computational effort, making it desirable to develop an inverse solution workflow that can use this advantage. However, compared to ST-ISW for T₂ of Bentheimer sandstone where the clay regions effective transverse relaxation time, effective diffusion coefficient, and transverse surface relaxivity $ρ2$ were determined, here the corresponding two longitudinal relaxation parameters would also be required, leading to a five-parameter inverse problem.

Incorporating the matching of T₁ distributions in addition to T₂ distributions following ST-ISW presented in Li et al. (2021) is straightforward: Both T₁ and T₂ fitting residuals can be jointly minimized by taking the sum of the two L2 norms of the fitting residuals as the objective. The inverse problem can then be solved by ST-ISW, except that incorporation of the T₁ distribution leads to additional physical constraints, which are helpful for the identification of the most realistic solution sets. This leads to a strategy mainly focused on extending the expected improvement (EI) or entropy search acquisition functions to the constrained optimization scenarios (e.g., the expected constrained improvement; Schonlau et al. 1998), and the predictive entropy search with constraints (Hernández-Lobato et al. 2015). The family of EI acquisition functions is superior for multimodal optimization as it automatically balances exploration and exploitation by taking into account both value information and uncertainty information. In this study, we identify the nonunique solution sets by multimodal maximization of the expected constrained improvement acquisition function using the multistart sequential quadratic programming (SQP).

While conceptually simple and easy to implement, ST-ISW fails to account for correlations between the two objectives, which may lead to a waste in computational resources. Historically, the (spatial) correlation between related objectives are mostly studied in geostatistics [e.g., Goovaerts (1997) used the cokriging approach to account for the spatial correlations between the concentration of the heavy metals that are expensive to evaluate and the variables that are easy to evaluate, for example, pH levels]. Cokriging with a linear model of coregionalization (LCM) is widely applied in geostatistics (Babak and Deutsch 2009; Ochie and Rotimi 2018; Al-Mudhafar 2018) and is then extensively studied and applied in transfer learning (Bonilla et al. 2007; Álvarez et al. 2012; Swersky et al. 2013). In this study, we incorporate a transfer learning ability to ST-ISW using LCM which jointly models the T₁ and T₂ objectives by leveraging information from each other. We propose two approaches that differentiate themselves by whether the two objectives (matching T₁ and T₂ distributions) are considered of equal importance: DT-ISW-S and DT-ISW-J, indicating if the two objectives are optimized sequentially or jointly; here DT stands for the dual-task. In sequential optimization, unknown physical quantities related to the objective with higher importance are identified first, setting the bounds for another objective, whereas in joint optimization all unknowns are identified simultaneously. DT-ISW-S is formulated for the scenario where one objective is considered of higher importance and it does not matter if the other objective strictly matches, whereas DT-ISW-J is for the scenario where the two objectives are considered of equal importance and good matches of the two objectives are desirable. We demonstrate the performance of DT-ISW-S and DT-ISW-J, with ST-ISW as a baseline by the identification of the five parameters of a Bentheimer sandstone sample: one shared quantity (i.e., the effective diffusion coefficient in clay regions), together with four unique quantities (i.e., surface relaxivity of quartz and effective relaxation time in clay regions for T₁ and T₂, respectively). We also provide the separate optimization (SEP) as the reference for which the physical quantities related to T₁ and T₂ are separately identified. Finally, we provide the POS, which allows operators to balance the preference of T₁ and T₂ data fits for the slightly conflicting objectives.

This paper is organized as follows: An overview of DT-ISW is presented, followed by a detailed description of the materials, NMR measurements, and NMR forward simulation. Next, a detailed description of DT-ISW for efficient modeling of the two related objectives using transfer learning is provided. Finally, the performances of DT-ISW-S and DT-ISW-J are investigated and compared against ST-ISW and SEP for a Bentheimer sandstone sample, determining five T₁ and T₂ relaxation-related intrinsic physical properties simultaneously.

There are five components in the developed DT-ISW for the inverse problem of T₁, T₂ parameter estimation, namely (a) the observables, (b) the forward NMR solver, (c) the cost-function, (d) the biobjective Bayesian optimization, and (e) the solution analysis, as presented in Fig. 1.

We briefly summarize these components in the itemized list below, while detailed descriptions follow in the next section.

(a) Measure NMR T₁ and T₂ relaxation decays $ml,exp$ and $mt,exp$⁠.

(b) Obtain the digitized representation of the sample via CT imaging (first iteration only), numerically derive the decays $ml,sim$ and $mt,sim$ for the chosen sets of unknown intrinsic parameters $x(T1)$ and $x(T2)$ (e.g., surface relaxivities of the resolved phase: quartz, effective diffusion coefficient of clay, effective relaxation times of clay regions only resolved on scanning electron miscroscope images) while honoring known quantities like porosity $ϕ$⁠, volumetric magnetic susceptibility $χν$⁠, bulk physical properties of water, and matching the physical kernel (see Eqs. 2 and 3) and regularization parameter $λ$⁠.

(c) Acquire the probability density distribution for both T₁ and T₂ using both simulated and measured decays [i.e., $ζsim(T1)$⁠, $ζexp(T1)$⁠, $ζsim(T2)$⁠, and $ζexp(T2)$]. Evaluate the cost-function and calculate the sum of the L2 norm of the fitting residuals for both T₁ and T₂ distributions.

(d) Jointly model the T₁ and T₂ objective functions using the vector-valued GP with correlation captured using the LCM, followed by maximum likelihood estimation of the model hyperparameters; then, the EI acquisition function is maximized to identify the promising candidates $x(T1)$ and $x(T2)$⁠. Finally, solutions within sparsely explored regions are evaluated using exact simulation.

(e) When the optimization budget is depleted, the objective functions are modeled with the updated model hyperparameters followed by a division of the solution space into various unique partitions (UPs) associated with different modes of the objective function using the solution space partitioning introduced in Li et al. (2021). Finally, for each UP, the top three solutions and the POS are reported.

NMR measurements were performed on a Bentheimer sandstone core plug (which is hereafter labeled as BHG2) fully saturated with 3 wt% NaCl brine (30,000 ppm). It was cored out of a rock block sourced from Bad Bentheim quarry (Lower Saxony, Germany) exhibiting 23.9% porosity. This sandstone is a remarkably clean quartz arenite, as five analyzed samples exhibit a narrow range of 97.94 to 98.07 wt% of silicon dioxide determined by X-ray fluorescence spectroscopy. The corresponding mineralogy associations are 95.67% quartz, 1.35% kaolinite, 2.39% feldspar, and 0.03% iron-rich components. Owing to the very small iron and manganese content, the rock is diamagnetic. The bulk volumetric magnetic susceptibility of three finely powdered samples (about 1 g each) measured using magnetic susceptibility balance (Sherwood Scientific, Cambridge, UK) is −7.45 × 10^–6 (SI units).

A 5-mm diameter core of Bentheimer sandstone was imaged on the UNSW Tyree X-ray CT facility’s HeliScan Micro-CT set equipped with a GE Phoenix NanoFocus X-ray source and a 3,072 × 3,072 flatbed detector. Radiographs were acquired over 12 hours with 80 kV X-ray energy using a 0.75-mm aluminum filter in double helix mode. The resulting spatial resolution of the 3D image was $ε=2.16$μm. Following reconstruction, the gray-scale tomogram was segmented into five phases: (macro) pore space, quartz, clay region (a microporous effective phase), feldspar, and an iron-rich dense mineral phase. The segmentation first identified pore space, intermediate intensity clay region, and high density (quartz and denser) in a standard converging active contour approach (Sheppard et al. 2004)—see Li et al. (2021) for the actual sample—followed by a further partitioning of the solid phase into quartz, feldspar, and an iron-rich phase. Feldspar and quartz grains were separated using grain partitioning by a watershed algorithm (Sheppard et al. 2004) followed by a thresholding based on average grain intensity, while the high-intensity iron-rich phase was extracted directly based on intensity. For reference, a slice through the tomogram and resultant phase segmentation is given in Fig. 2.

Statistics of the Bentheimer sandstone segmentations as well as magnetic susceptibility values for the calculation of the internal magnetic field at lattice resolution ε that follows is listed in Table 1.

Since we focus on NMR relaxation at low external field strength of 0.047 T and given the small amount of feldspar, we assume that the two phases, feldspar and quartz, share the same surface relaxivity and volumetric magnetic susceptibility, so that the increased magnetic susceptibility is assigned to the iron content in the high-density phase, leading to a higher $χν,h$⁠, with the clay region considered as a mixture of kaolinite and NaCl brine, as shown in Table 1.

Here we measured T₁ relaxation time using the standard inversion recovery pulse sequence (Carr and Purcell 1954). It is formed by a 180° pulse which inverts the magnetization, followed by recovery at time delay $τ$⁠, and concluded by a 90° pulse that projects the attenuated magnetization on the transverse plane for the acquisition of the free induction decay. By varying time intervals $τ$ in subsequent scans, a set of free induction decays with different amplitudes is obtained. T₁-inversion recovery measurements were performed with the following set of parameters: 64 time increments $τ$ were spaced logarithmically between 0.2 and 10,000 ms. The signal was acquired as an average of four scans following a phase cycling sequence to suppress systematic noise due to hardware imbalances. Recording the initial value of each free induction decay as a function of the recovery time results in a mono-, or in a general case, in a multiexponential decay described by the following expression with the T₁-inversion recovery kernel $1-2e-t/T1$as (Callaghan 1994; Song et al. 2002; Wang et al. 2017)

where $mT1(t)$ is the longitudinal magnetization decay as a function of time $t$⁠, $ηn$ is noise assumed white and additive, and our objective is the estimation of the probability density function of T₁ [i.e., $ζ(T1)$].

The standard method to measure transverse relaxation time T₂ is the Carr, Purcell, Meiboom, and Gill (CPMG) pulse sequence (Carr and Purcell 1954; Meiboom and Gill 1958). It consists of a 90° pulse that projects the initial magnetization on the transverse plane, followed by a series of 180° refocusing pulses separated by a time interval of $tE=2τE$ (interecho time). Each 180° pulse produces an attenuated signal (echo) with the maximum occurring at time $τE$ after the pulse. By acquiring the top of each echo, the exponential (or multiexponential, in a more general case of porous media) signal decay is obtained with the decaying rate equal to the inverse of the transverse relaxation time T₂, thus giving an $e-t/T2$ kernel for the multiexponential problem as (Callaghan 1994; Song et al. 2002)

where $mT2(t)$ is the transverse magnetization decay as a function of time $t$ and $ηn$ is additive white noise, and the objective is the estimation of the probability density function of T₂ [i.e., $ζ(T2)$]. All NMR experiments were performed using the Magritek 2 MHz Rock Core Analyzer. Being equipped with the P54 RF probe (54-mm bore opening), it is capable of achieving signal-to-noise ratio in excess of 50 from four scans for the most saturated rock samples of about 20 cm³ in volume. The CPMG measurements were recorded with an echo spacing of 0.2 ms and 50,000 pulses per echo-train; the echo amplitude was determined by averaging 16 points on top of each echo.

The resultant magnetization decays are converted to relaxation time distributions by solving the respective Fredholm integrals of the first kind (i.e., Eqs. 2 and 3) for the probability density functions $ζ(T1)$ and $ζ(T2)$⁠. Both are nonnegative least squares problems, for which various numerical solvers are available (Bertsekas 1982; Hansen 1992; Lawson and Hanson 1995). For consistency, all experimental and simulated decays use exactly the same kernels with 256 logarithmically spaced bins for both T₁ and T₂ and are inverted using the algorithm of Lawson and Hanson (1995).

The NMR simulations are carried out via a random walk method. Some modifications to our previous implementations have been made over time. Arns (2004) utilized a fixed displacement and fixed timestep random walk following Mendelson (1990) and Bergman et al. (1995) by directly using the natural resolution of the underlying micro-CT image, with the addition of accounting for dephasing. Arns et al. (2011) evolved that approach to one with fixed displacement, but variable timestep while also implementing a higher time resolution through a logical subgridding approach: Given the tomogram voxel size ε, random walks are carried out on a lattice with resolution $εw=ε/ls$⁠, where $εw$ is the step length of the random walk and l_s is a grid refinement factor. Furthermore, the approach allows for the existence of effective properties (e.g., voxels which are microporous and exhibit an effective diffusion coefficient). Alhwety et al. (2016) utilized the approach for a dual-scale synthetic laminated sphere pack, tracking the magnetization exchange of different microporosity domains with distance maps embedded into the NMR simulations for spatial analysis. SayedAkram et al. (2016) utilized effective medium approaches for describing the motion and NMR relaxation in Indiana limestone exhibiting significant microporosity, comparing experiment and simulation under mixed-wet conditions. Shikhov et al. (2017) measured and modeled the NMR T₁-z spatially resolved NMR relaxation response for a set of sandstones under centrifuge desaturation conditions resulting in good agreement between experiment and simulation. More recently, Zheng et al. (2018) considered the case of flow propagators and exchange of magnetization between microporous and macroporous regions of a glass-bead sample. The effective medium approach for nonresolved porosity was also utilized in several studies (Cui et al. 2021; Li et al. 2021). In this work, modifications are added, allowing the simultaneous extraction of T₁ and T₂ responses, which adds significant additional analysis options to the approach while also enabling the efficient determination of intrinsic parameters enabling these applications. For this reason, we summarize our method in the following in more detail.

Consider a microstructure given by a 3D regular grid from a micro-CT image segmented into $N$ phases. Then the material properties at voxel scale relevant for NMR response simulation can be presented by a set of arrays $T1b∈ℜ1×N$⁠, $T2b∈ℜ1×N$⁠, $χν∈ℜ1×N$⁠, $ρ1∈ℜN×N$⁠, $ρ2∈ℜN×N. De∈ℜN×N$⁠. $De∈ℜN×N$⁠. Here, $T1b$ is longitudinal relaxation time, $T2b$ is transverse relaxation time, $χν$ is volumetric susceptibility, $ρ1$ is longitudinal surface relaxivity, $ρ2$ is transverse surface relaxivity, and $De$ is the effective diffusion coefficient. We note that the array entries represent either intrinsic properties for fully resolved voxels or effective properties. For relaxation-related parameters, we have the element-wise conditions $T2b≤T1b$ and $ρ2≥ρ1$⁠. For a single fluid and a fully saturated sample considered here, random walks start only in voxels where porosity is greater than zero and with a probability proportional to that of porosity. Consequently, the introduced arrays are sparse—not all entries are meaningful and if a full notation is utilized for convenience, those entries may be set to zero. To translate the material properties into timesteps describing the motion and movement and stepwise signal damping factors, we introduce a set of additional quantities, namely $Δt∈ℜN×N$⁠, $pmov∈ℜN×N$⁠, $ST1∈ℜN×N$⁠, and $ST2∈ℜN×N$ as timestep, movement probability, and T₁ and T₂ signal damping factors. A random walker at step count $s$ with position $rs$ in a voxel of phase $p$ and internal clock t_s attempts a step in a random direction (one of six nearest neighbors). If the attempted motion is successful, the random walker changes to position $rs+1$ of phase $q$ with timestep $Δtpq=εw26Dpq$ and signal damping factor following Bergman et al. (1995) as

advancing the clock to $ts+1=ts+Δtpq$⁠. Here, the terms in round brackets correspond to survival probabilities for T₁ and T₂ surface relaxation processes, respectively, while the exponential terms account for bulk relaxation. The surface relaxivities are defined in the usual way (e.g., independently for a T₁ or T₂ experiment/simulation). This implies $ρ2≥ρ1$⁠, as a T₁ surface relaxation will also destroy phase coherence. When jointly simulating T₁ and T₂ relaxation, the relaxation already accounted for by T₁ processes therefore needs to be corrected for in the T₂ process, leading to the prefactor $1/ST1,pq$ in Eq. 5. If the attempted motion is not successful due to $pmov,pq§lt;1$⁠, the damping factors $ST1,pq$ and $ST2,pq$ are still applied for the attempted move. However, the timestep is taken to be $Δtpp$ and the internal field (see below) is considered to be the one of the current position. The above formulation allows one to write the signal decay of an individual walker $w$ carrying out $Ns$ steps with initial signal strength M₀ as a product, while recording the T₁ and T₂ decays simultaneously

Here, $ϕD$ is the accumulated phase shift due to diffusion in an inhomogeneous (offset) internal magnetic field $Bi$⁠, with total magnetic field $B(r)=Bi(r)+B0$⁠. The latter is calculated in the dipole approximation by convolution of the susceptibility field $χν(r)$ with the dipole field $Bdip$ as

Details of the implementation are given in Arns et al. (2011). Following the random path of a walker over steps $j$ and weighting the phase accumulation due to diffusion in the internal field by the timestep for changing locations, the total phase shift after $Ns$ steps can be written as

$r0$and $rj$ denote the location of the start and later positions of the random walker. The CPMG sequence is implemented by starting random walks with phase $ΔϕD=0$ and switching the sign of $ΔϕD$ at times $(nE+12)tE$⁠, $nE=1,2,…NE$⁠, where $tE$ is the echo spacing and $NE$ is the total number of echoes recorded. Recording at times $t=nEtE$⁠, $nE=1,2,…NE$⁠, and averaging over random walks result in the magnetization decays $MT1(t)$ and $MT2(t)$⁠. We note that $t$ is an integer multiple of $tE$⁠. As timesteps are allowed to vary, individual timesteps are split into two where necessary by linear interpolation for accurate timestepping of the CPMG sequence. Finally, a comment on the choice of p_mov: We assume continuity of the pore space (maximal correlation) at interfaces. This leads to $pmov,pq=1$ for $ϕsim,p≤ϕsim,q$⁠, $pmov,pq=ϕsim,q/ϕsim,p$ for $ϕsim,p§gt;ϕsim,q$⁠, and $pmov,pq=0$ otherwise, where $ϕsim,p$ denotes the calculated porosity of phase $p$ from simulation. Furthermore, we determine the off-diagonal entries of $De$ by taking the harmonic mean of the respective nonzero diagonal entries. This is the simplest assumption (e.g., treating the interface as two parallel half-spaces).

The specific problem of minimization of the sum of the L2 norm of the fitting residuals for both T₁ and T₂ distributions can be formulated as:

where $xl$ and $xu$ are the lower and upper bounds of the prespecified search domain and $f(T1)$ and $f(T2)$ are L2 norm of the fitting residuals for T₁ and T₂ distributions, respectively, which may be expressed as

with weights for all data points of the T₁ and T₂ distribution being equally 1. Each lower-dimension component $x(Ti),i∈{1,2}$ comprises the unique components $x(Ti,unique)=[ρi,q,log10⁡Tie,c]$ and the shared component $x(shared)=log10⁡De,c$⁠, and therefore Eq. 11 specifies the complete solution. The constraints Eqs. 14 and 15 apply because in most practical situations for saturated porous media, T₂ relaxation rate is faster than T₁ due to near-surface interactions specific to transverse relaxation. An important nature of the problem is that the NMR simulator is adapted in a way that both $f(T1)$ and $f(T2)$ can be simulated simultaneously in a single function evaluation, enabling individual modeling of the two objectives to better reveal the structure of the true objective.

GP is a distribution over functions and is characterized by its mean function $m(x)$ and covariance function $k(x,x′)$ as (Rasmussen and Williams 2005; Li et al. 2021)

where $x=(x1,x2,…,xD)⊤$ is a candidate in $𝒳$ of dimension $D×1$⁠, $f(x)$ is the process evaluated at location $x$⁠. In GP, the prior distribution over $f$ is expressed as

where $𝒩$ denotes normal distribution, $m0∈ℜn×1$ denotes the mean vector for observations, $(X,𝐲)={xi,yi}i=1n$ and $(x*,y*)$ denotes observed and predicted candidates, respectively, and $Kf(X,X)$ denotes the $n×n$ covariance matrix for $f$⁠, with the $(p,q)$ entry expressed by $kSE(xp,xq)$⁠. Using the Gaussian likelihood function to accommodate for the noise in observations $𝐲$ yields $y|f=N(f,σn2I)$⁠. As we are interested in prediction of a single candidate, the joint distribution of the observations $𝐲$ and a single prediction $y*$ can be expressed as

where $Kf(X,x*)$ is the $n×1$ covariance between observed and predicted candidates, $Kf(x*,x*)$ is the variance of the predicted candidate, $m0*$ is the mean for $y*$⁠, and $I$ is the $n×n$ identity matrix. By conditioning the predictions on observations from Eq. 21, the posterior distribution for a single predicted candidate becomes

where $𝒩$ denotes normal distribution, $θs$ is a vector of model hyperparameters for the scalar-valued GP, and

The vector-valued GP indicates GP with multiple outputs. We use $Q$ to represent the number of correlated objective functions to be modeled, with $Q=2$ applied throughout this study. To model the correlations between the $Q$ outputs of $f$⁠, the GP kernel in Eq. 19 has to be replaced by a vector-valued GP kernel, so that the vector-valued function $f$ can be expressed as (Bonilla et al. 2007; Swersky et al. 2013)

where $m(x)∈ℜQ×1$ denotes the stacked mean function ${mq(x)}q=1Q$ for each output and $𝐊$ denotes a matrix-valued function whose entries $(𝐊(x,x′))q,q′$ express the correlation not only between the input pairs $x$ and $x′$⁠, but also between the outputs pairs $fq(x)$ and $fq′(x′)$⁠. Similar as in Li et al. (2021), the null hypothesis is assumed as there is usually no prior knowledge about our problem; as a result, we use a constant mean function for each objective throughout this study (i.e., $m(x)=m0$⁠), with a stationary covariance function: parameterized squared exponential for each objective. The squared exponential kernel $Kf,q′(X,X)$ is parameterized using automatic relevance determination, that is,

where $σf2$ is the signal amplitude, $Λ=diag(ℓ)-2$ accommodates different levels of sensitivity against each component of $x$⁠, with $ℓ=(l1,…,lD)⊤$ being a vector of positive values, referred to as length scales. Automatic relevance determination has been incorporated to infer the sensitivity against each dimension of the input, which has been used successfully for removing irrelevant input associated with a very large inferred length scale (Williams and Rasmussen 1995; Rasmussen and Williams 2005; Wipf and Nagarajan 2007). $σf2$ is set to 1 when applied to vector-valued GP.

Various coregionalization models have been proposed in literature (Goovaerts 1997; Gelfand et al. 2004; Majumdar and Gelfand 2007; Krupskii and Joe 2015), and here we use the standard LCM originated from geostatistics to model $f$⁠, with the complete derivation of LCM shown in  Appendix A. Assume that there are $n$ distinct observed candidates and observations for each objective [i.e.,$𝒟obs=(X,𝐲)={(xt,i,yt,i)|i=1,…,n,q=1,…,Q}$]. Using LCM, the vector-valued GP kernel for the latent variable $f$ can be expressed as (Álvarez et al. 2012)

where $⊗$ is the Kronecker product between matrices, $Q′$ is the number of separable kernels, $Bq′$ is a $Q×Q$ covariance matrix measuring the relationship between the $Q$ outputs, and $Kf,q′(X,X)$ is an $n×n$ covariance matrix measuring the relationship between inputs, with the $(i,j)$ entry specified by $kSE(xi,xj)$⁠. Typically $Q′≤Q$ and Eq. 26 degrades to intrinsic coregionalization model if $Q′=1$⁠. Both matrices are symmetric and positive semidefinite. We use a full-rank LCM (i.e., $Q′=Q$⁠) throughout this study.

Eq. 26 is called the sum of separable kernels as the covariance function $𝐊f(X,X)$ is expressed as the sum of $Q′$ product of two covariance functions. It is observed from Eq. 26 that for each latent function, the correlation within tasks and between tasks are decoupled. Similar to the scalar-valued GP case, we use a Gaussian likelihood to account for the noisy observations (i.e., $𝐲|f=𝒩(f,Σn⊗I)$⁠, where $Σn=diag⁡(σn2)$ is a $Q×Q$ diagonal matrix accommodating task-dependent noise levels, and $σn2=(σn,12,…,σn,Q2)⊤$ is a vector of noise variances). Similar to the case of scalar-valued GP, the posterior distribution of mean and covariance for a single predicted candidate $x*$ can be expressed as

where $𝐲*=[𝐲1*,…,𝐲Q*]⊤$ concatenates $Q$ outputs for $x*$⁠, $θv$ is a vector of model hyperparameters where $v$ stands for vector-valued GP, and

where $m0*$ denotes the stacked constant mean function for each objective, $𝐊f(X,x*)$ is the covariance matrix for the latent variable $f$ evaluated at all pairs of observed candidates and predicted candidate, $𝐊f(x*,x*)$ is the covariance matrix for the predicted candidate, $μx*=[μ1(x*),…,μQ(x*)]⊤$ concatenates $Q$ outputs, and $Σx*∈ℜQ×Q$ is the covariance matrix for $Q$ tasks with variances on the diagonal. Fig. 3 shows a comparison of the prediction and uncertainty quantification between scalar-valued GP and vector-valued GP.

It is clear that given unevenly distributed observations for both tasks, the scalar-valued GP cannot truly reflect the unexplored region, whereas the vector-valued GP predicts the trend in unexplored region using knowledge from the related objective, leading to a more informative EI guiding the search.

Proposed by Mockus (1975), the EI acquisition function has been widely applied in Bayesian optimization with the standard GP, balancing exploration and exploitation by taking information on value and local uncertainty (Huang et al. 2006; Forrester et al. 2006; Quan et al. 2013). It is defined as the amount of improvement that the evaluation of a candidate is expected to induce.

Extending the EI to the multiobjective case is straightforward. For vector-valued GP, the EI acquisition function associated with the $qth$ component can be expressed as

where $xq-$ is the current best observed candidate for the $qth$ objective and $uq=yq*|X,𝐲,x*,θv$ is the short-hand symbol for the Gaussian distributed variable given by Eq. 27. Eqs. 27 and 29 yield

where $Φ(x)$ and $φ(x)$ are the cumulative distribution function and the probability density function of the standard normal distribution respectively.

In each iteration, the choice of the promising candidate is based on maximization of the EI acquisition function. The order of maximization requires a choice: Sometimes the T₁ data is only used for verification of the nonunique solution sets obtained from the T₂ inverse problem, whereas sometimes T₁ and T₂ data are considered of equal importance and it might be desirable to have both T₁ and T₂ distributions fitted of equal quality. Formally, we propose two routines for maximization of the EI acquisition function corresponding to these two scenarios:

• $aEI,T2$ is first maximized, followed by maximization of $aEI,T1$ in a subspace constrained by Eqs. 14 and 15. For clarity we denote DT-ISW with sequential maximization of $aEI,Tq$ as DT-ISW-S.

• $aEI,T1$ and $aEI,T2$ are jointly maximized by maximization of the sum of $aEI,T1$ and $aEI,T2$⁠. For clarity, we denote DT-ISW with joint maximization of $aEI,Tq$ as DT-ISW-J.

Throughout this study, $f(T2)$ is assigned more importance as the T₂ decays acquired using CPMG measurement consist of significantly more data points than the T₁ decays acquired during T₁ inversion recovery measurement, making T₂ distributions and the associated objective function more reliable. Maximization of $aEI,Tq$ under linear constraints can be tackled using SQP (Boggs and Tolle 1995). The initial guesses are placed randomly in the feasible domain, ensuring that both the global optimal solution and major local optimal solutions are identified. The performance of DT-ISW-J and DT-ISW-S is compared in the next section.

After the optimization budget is depleted, the model hyperparameters are updated and the T₁, T₂ objectives are jointly modeled under the new set of hyperparameters. The solution space of Eq. 10 is then divided into various UPs corresponding to various modes of the objective function. The solution with the lowest fitness value in each UP is referred to as the local minimum (LM). For details, see Li et al. (2021). As the two objectives involve noisy measurements which are never noise-free, it is expected that they will be slightly conflicting and there is no solution that simultaneously optimizes each objective precisely. To account for this, we provide the POS of the T₁ and T₂ data fits for the slightly conflicting objectives, as shown in Fig. 4.

POS is defined as the set of solutions that is not dominated by any other solutions. Generally speaking, $x1$ dominates $x2$ if $x1$ beats or ties $x2$ for each of the two objectives; mathematically, it is expressed as $f(i)(x1)≤f(i)(x2)$⁠, $∀i∈{1,2}$⁠, and $f(i)(x1)<f(i)(x2)$⁠, $∃i∈{1,2}$⁠.

In this study, we demonstrate the performance of the three approaches: DT-ISW-J, DT-ISW-S, and ST-ISW on Bentheimer sandstone, identifying five parameters by computationally producing T₁ and T₂ distributions almost identical with experiments using the dual-output NMR simulator. The individual optimization for the T₁ and T₂ are presented as reference methods. As was done in Li et al. (2021), practically the base-10 logarithm of the effective relaxation time and of the effective diffusion coefficient are used since they usually vary across orders of magnitude, and the search for $x$ is limited to a domain bounded by the following constraints:

• 0 $≤$$ρ1,q$/(μm/s), $ρ2,q$/(μm/s) $≤$ 16

• −3 $≤$$log10⁡(T1e,c/s)$⁠, $log10⁡(T2e,c/s)$$≤$ 0

• −6.28 $≤$$log10⁡[De,c/(cm2/s)]$$≤$ −4

and the corresponding 5D solution space $𝒳$ is rather large so that a grid search or random search is computationally intractable. We compare the performance of the following three approaches discussed in the previous section for the T₁, T₂ inverse problem: DT-ISW-S, DT-ISW-J, and ST-ISW. As the reference method, the SEP minimizes T₁ and T₂ independently in the 3D solution space as was done in Li et al. (2021), ignoring that $log10⁡De,c$ is a shared value. We summarize the description of the three approaches and the reference in Table 2 and note that the determined values $ρ1,q$ and $ρ2,q$ are the lattice-based discrete surface relaxivities of quartz, whereas the actual continuum values are $A$ times higher, and $A$ is 3/2 following Bergman et al. (1995). The optimization budget is 180 function evaluations.

Fig. 5a displays the best fitness value so far for $f(obj)$⁠, together with its two components $f(T1)$ and $f(T2)$ in Figs. 5b and 5c for the three approaches, with SEP plotted as a reference. In terms of the objective $f(obj)$⁠, both DT-ISW-S and DT-ISW-J are more efficient in the first 60 function evaluations, with the speed of convergence almost twice as fast as ST-ISW. In the following 120 steps, DT-ISW-J maintains the same rate of convergence until Step 83, where it finds the lowest fitness value, whereas DT-ISW-S continues to minimize $f(obj)$ and find a lower $f(obj)$ at Step 152. ST-ISW maintains the same convergence speed until the optimization budget is depleted. Minimization of $f(T1)$ is very similar to the case of minimization of $f(obj)$⁠: In the first 32 steps, DT-ISW-J immediately reduces $f(T1)$⁠, whereas due to low priority assigned to T₁ fitting, DT-ISW-S cannot reduce $f(T1)$ since a good proposal of $x(T1)$ is subjected to the valid proposal of $x(T2)$⁠. DT-ISW-S becomes effective since Step 32, and the overall rate of convergence is comparable to DT-ISW-J. However, ST-ISW, with a lower rate of convergence, steadily minimizes $f(T1)$ to the same level as DT-ISW-S, which is slightly higher than 1×10⁻³. All approaches identified lower $f(T1)$ than the reference methods.

In terms of $f(T2)$⁠, DT-ISW-J, DT-ISW-S, and SEP-T₂ maintain a similar speed of convergence until the fitness values fall below 1$×10-3$ before Step 70. In the following steps, both SEP and DT-ISW-S achieve lower fitness values because they can freely evaluate in the solution space without constraints. Though at a slower rate of convergence, ST-ISW steadily reduces $f(T2)$ to the level of 1$×10-3$ until its last step. It can be inferred from Fig. 5 that the speed of convergence for both DT-ISW-S and DT-ISW-J are more than twice as fast as ST-ISW. DT-ISW-J achieves good T₁ and T₂ distribution agreement, whereas DT-ISW-S achieves better agreement for T₂ and slightly worse agreement for T₁.

As is described in the previous section, ST-ISW caters for the multimodal nature of the objective functions, with solution space partitioning applied in the solution analysis to cluster all solutions to different UPs. Here we apply the same solution analysis procedure to both DT-ISW-S and DT-ISW-J. Table 3 shows identified parameter values of the LMs for the three approaches and the reference. LMs with $f(obj)$ above 1$×10-2$ are not listed. All three approaches identified a consistent value of $x$⁠, and all are capable of achieving lower $f(T1)$ values than the reference. SEP-T₂ is affected by the problem of nonunique solutions—although achieving a lower $f(T2)$⁠, the LM 1 identified by SEP-T₂ is less physically valid than LM 2. As a comparison, all approaches identified at least 2 UPs, but from Table 3 we see that for each approach there is only one UP within which the solutions are both mathematically sound and physically valid.