Interpretation of deep directional resistivity measurements acquired in high-angle and horizontal wells using 3-D inversion

Interpretation of deep directional resistivity measurements acquired in high-angle and horizontal... SUMMARY The interpretation of resistivity measurements acquired in high-angle and horizontal wells is a critical technical problem in formation evaluation. We develop an efficient parallel 3-D inversion method to estimate the spatial distribution of electrical resistivity in the neighbourhood of a well from deep directional electromagnetic induction measurements. The methodology places no restriction on the spatial distribution of the electrical resistivity around arbitrary well trajectories. The fast forward modelling of triaxial induction measurements performed with multiple transmitter–receiver configurations employs a parallel direct solver. The inversion uses a pre-conditioned gradient-based method whose accuracy is improved using the Wolfe conditions to estimate optimal step lengths at each iteration. The large transmitter–receiver offsets, used in the latest generation of commercial directional resistivity tools, improve the depth of investigation to over 30 m from the wellbore. Several challenging synthetic examples confirm the feasibility of the full 3-D inversion-based interpretations for these distances, hence enabling the integration of resistivity measurements with seismic amplitude data to improve the forecast of the petrophysical and fluid properties. Employing parallel direct solvers for the triaxial induction problems allows for large reductions in computational effort, thereby opening the possibility to invert multiposition 3-D data in practical CPU times. Controlled source electromagnetics (CSEM), Electromagnetic theory, Inverse theory, Numerical approximations and analysis, Numerical modelling INTRODUCTION With hydrocarbon exploration and production continuously moving into areas of increasing geological and fluid complexity, the oil and gas industry is adopting new technologies to improve the efficiency of the well placement. One of the most important properties of geological formations is their electrical resistivity, which has high sensitivity to both porosity and hydrocarbon saturation. Given the limitations of current electromagnetic (EM) technologies such as logging-while-drilling (LWD), new monitoring and exploration methods were recently introduced. A new generation of deep directional resistivity (DDR) tools provides accurate estimations of various petrophysical properties of hydrocarbon-producing rock formations and has an improved depth of investigation of more than 30 m from the wellbore (Ezioba & Denichou 2014; Seydoux et al. 2014; Dupuis & Denichou 2015). This technology, often referred to as reservoir mapping-while-drilling, several times exceeds the spatial coverage of conventional LWD tools and offers an unprecedented opportunity to map the near-wellbore reservoir structure. In turn, the improved ability to navigate well trajectories with respect to geological and fluid boundaries using real-time information allows operators to reduce drilling risk, maximize reservoir exposure, and increase the overall production potential. With recent advances in drilling technology, long high-angle (HA) and horizontal (HZ) wells have become common in the oil and gas industry (Ijasan et al. 2014). The accuracy of landing HZ wells in a reservoir embedded in a complex geology often determines the overall financial viability. Many of the conventional LWD modelling methods were developed for vertical wells and are not suitable for interpretation of DDR measurements in HA/HZ wells. The real-time processing of DDR measurements in complex geologic media requires a full 3-D inversion method. An efficient inversion will open the possibility of estimating 3-D distributions of electrical conductivity. However, one of the main practical limitations of this approach is its high computational cost. Because the DDR tool is moving in the well, modelling needs to consider a large number of transmitter and receiver positions (in the range of thousands of samples when the position and the three independent emission directions are considered) for a long HA/HZ well. The inversion process, in turn, involves a high number of forward modelling operations. Thus, the total computational cost is mainly influenced by the cost of the forward modelling. Multiple 3-D modelling methods of multicomponent borehole EM induction logging measurements have been developed during the past decades (e.g. Newman & Alumbaugh 2002; Davydycheva et al. 2003; Hou et al. 2006; among others). However, the complexity of 3-D models dramatically increases the computational requirements for modelling and hence for inversion. The size of these problems has often been considered excessively large to be solved in real time for practical applications, thereby preventing the widespread use of full 3-D methods for the interpretation of borehole induction measurements. Many of the published borehole induction modelling methods take advantage of the axial symmetry of formation geometry and employ 2-D or 2.5D modelling and inversion algorithms (Pardo et al. 2006; Abubakar et al. 2008; Wang et al. 2009; Pardo et al. 2013; Wiese et al. 2015). Other approximations have been employed as well to limit the large number of model parameters and to mitigate the computational cost of inversion. For example, recently, Thiel et al. (2016) documented a fast 2-D inversion that utilized pointwise 1-D forward simulations for real-time interpretation of DDR measurements for MWD. Bensdorp et al. (2016) described the inversion for 3-D conductivity distributions based on single spherical scatterer approximations. Bakr et al. (2017) used 1.5D simulations for fast inversion of 3-D transversely isotropic formations approximated by planar layers. With the development of high-performance computing methods and the increase in the availability of computer resources, 3-D inversion of geophysical data sets no longer poses an insurmountable computational challenge. While still requiring significant computational resources for large-scale exploration problems, inversion of local data sets can be performed on a workstation or a small computer cluster in practical CPU times, that is, hours or even minutes. In particular, recent advances in parallel direct solvers allow us to improve the efficiency of simulations of multiposition and multicomponent triaxial induction measurements in complex reservoir formations with extreme resistivity contrasts (Puzyrev & Torres-Verdín 2016). The main limitation of these solution techniques is their large computational and memory demand; however, many state-of-the-art EM modelling problems involve up to several million unknowns and thus can be efficiently solved with direct solvers. Recent work by Shantsev et al. (2017) presents a large-scale 3-D controlled-source EM problem with more than twenty million unknowns solved with a block low-rank multifrontal direct solver. In this paper, we describe a full 3-D parallel inversion method to interpret multicomponent DDR measurements acquired in HA/HZ and apply it to several synthetic models. The remainder of this paper is organized as follows. In the next section, we summarize the fundamentals of the inverse problem. Then, we discuss the numerical approach used for forward modelling of multicomponent DDR measurements and specific details of the minimization algorithm such as pre-conditioning and optimal step length estimation. In the Examples section, we study the sensitivity of DDR tools to determine the oil-water contact during production (monitoring problem) and to locate favourable sandstone targets (exploration problem). Finally, the last section summarizes the conclusions and implications of our study. THEORY The inverse problem is formulated as a regularized nonlinear minimization problem with the following quadratic cost function:   \begin{equation}\phi {{\bf (m)}} = \frac{1}{2}\left\| {{{\bf D}}\left( {{{\bf F}}({{\bf m}}) - {{{\bf d}}_{{{\bf obs}}}}} \right)} \right\|_2^2\,\, + \,\,\frac{1}{2}\lambda \,{{\bf R}}({{\bf m}}),\end{equation} (1)where the first term in the right-hand side is the L2-norm based data misfit, dobs is the vector of measurements (observations), F(m) is the forward problem mapping, D is a diagonal data weighting matrix, R (m) is the stabilizing functional, and λ is the Lagrange multiplier. In this study, we use the pre-conditioned nonlinear conjugate gradient (NLCG) method to minimize the misfit between the measurements and their numerical simulations. This algorithm is efficient for parallel computers and has a small memory footprint compared to second-order minimization schemes (Rodi & Mackie 2001; Nocedal & Wright 2006). The stabilizing functional ensures the well-posedness of the nonlinear inverse problem in the presence of noise and/or inadequate measurements. Various regularization methods exist in the literature (Farquharson & Oldenburg 2004; Zhdanov 2009). In the examples described below, we use the Tikhonov minimum-norm regularization, given by   \begin{equation}{{\bf R}}\,({{\bf m}}) = \left\| {{{\bf m}} - {{{\bf m}}_{{{\bf ref}}}}} \right\|_2^2,\end{equation} (2)where vectors m and mref designate the parameters of the current and the reference models, respectively. The model parameters are commonly chosen as the logarithm of the conductivity σ, thus decreasing large gradient variations   \begin{equation}{{\bf m}} = \ln {{\boldsymbol \sigma }}, \quad {{{\bf m}}_{{{\bf ref}}}} = \ln {{{\boldsymbol \sigma }}_{{{\bf ref}}}}.\end{equation} (3) The reference model contains a priori information about subsurface properties, which is especially important for spatially complex geologies. Because one can often choose a good reference model mref and expect sharp resistivity contrasts in the conditions under analysis, regularization (2) has been found more suitable than, for example, maximum smoothness methods based on L2 or L1 norms. The inversion algorithm estimates formation resistivity by progressively minimizing the misfit function (1). The data gradient (i.e. the derivative of the data term in eq. (1) with respect to the model parameters, m) is given by (Newman & Boggs 2004)   \begin{equation}{{{\bf g}}_{{\bf d}}} = \nabla {\phi _{{\bf d}}}({{\bf m}}) = - {\mathop{\rm Re}\nolimits} \left[ {{{\left( {{{\bf DJ}}} \right)}^T}{{\left[ {{{\bf D}}\left( {{{\bf F}}({{\bf m}}) - {{{\bf d}}_{{{\bf obs}}}}} \right)} \right]}^*}} \right],\end{equation} (4)where * stands for complex conjugation. Deriving a computationally efficient form of eq. (4) is critical for a robust inverse solution (McGillivray & Oldenburg 1990). The adjoint method to calculate the data gradient has been used previously in EM (Newman & Alumbaugh 1997) and seismic full waveform inversion (Plessix 2006). This procedure allows us to avoid the explicit calculation and storage of the Jacobian matrix, J. When a direct solver is employed, the adjoint field can be calculated at a relatively low cost since the EM inverse problem is self-adjoint. At each iteration step, after the frequency-dependent forward problem F(m) is solved, we used the system matrix factorization to calculate the misfit gradient using the adjoint method. Thus, the expensive factorization is performed once for each frequency at the beginning of each iteration and then it is applied in the forward and the adjoint solves throughout this iteration. The data misfits are added over all frequencies and pre-conditioned to calculate the conjugate gradient and the search direction uk + 1 given by   \begin{equation}{{{\bf u}}_{{{\bf k + 1}}}} = - {{{\bf g}}_{{{\bf k + 1}}}} + {\beta _k}{{{\bf u}}_{{\bf k}}},\quad {{{\bf u}}_{{\bf 0}}} = - {{{\bf g}}_{{\bf 0}}}.\end{equation} (5)Here, the scalar βk is the conjugate gradient update parameter. Different conjugate gradient methods correspond to different choices for βk. We refer the reader to Hager & Zhang (2006) for a detailed chronological description of the update parameters. For large-scale problems, we prefer a βk that does not require the evaluation of the Hessian matrix. We use the standard Polak–Ribiere formula in our geophysical inversion algorithm (Polak & Ribiere 1969). Finally, the optimal step length αk is estimated as described below and the model is updated as   \begin{equation}{{{\bf m}}_{{{\bf k + 1}}}} = {{{\bf m}}_{{\bf k}}} + {\alpha _k}{{{\bf u}}_{{\bf k}}}.\end{equation} (6) This iterative process starts from the starting model m0 which can be (but not necessarily is) chosen the same as the reference model. FORWARD MODELLING To obtain the simulated measurements d = F(m), we solve the forward problem F(m) that is described by the time-harmonic diffusive Maxwell equations   \begin{equation}\nabla \times {{\bf E}} = - \,i\omega {\mu _0}{{\bf H}} - {{{\bf M}}_{{\bf s}}},\quad \nabla \times {{\bf H}} = \bar{\sigma }{{\bf E}} + {{{\bf J}}_{{\bf s}}},\end{equation} (7)where ω is angular frequency, $$\bar{\sigma }({{\bf r}})$$ is electric conductivity tensor, μ0 is magnetic permeability of free space, Js and Ms are the electric and magnetic sources, respectively (Ward & Hohmann 1988). Taking the curl of the first equation in eq. (7) and substituting it in the second one gives the curl-curl equation for the electric field   \begin{equation}\nabla \times \nabla \times {{\bf E}} + i\omega {\mu _0}\bar{\sigma }{{\bf E}} = - i\omega {\mu _0}{{{\bf J}}_{{\bf s}}} - \nabla \times {{{\bf M}}_{{\bf s}}}.\end{equation} (8) Reliable inversion and interpretation of multicomponent induction data in complex geological environments requires accurate and efficient forward modelling. In our implementation, we approximate eq. (8) by finite differences based on a staggered-grid method that is commonly applied in geophysical inversion (e.g. Newman & Alumbaugh 2002; Hou et al. 2006; Abubakar et al. 2008). The scattered-field formulation of eq. (8) is used to avoid excessive meshing near the sources. Our method handles large-scale problems in parallel and was successfully tested on complex isotropic and transversely isotropic reservoir models with extreme resistivity contrasts. The results were in good agreement with those obtained from a semi-analytic 1-D approach, thus confirming the accuracy and reliability of the simulation method. Arbitrary anisotropic media could be accurately simulated by using Lebedev grids (Davydycheva et al. 2003; Jaysaval et al. 2016). Modern DDR tools provide multicomponent measurements at various frequencies and transmitter–receiver offsets to estimate the complexity of the surrounding formations with better spatial resolution. Triaxial induction tools introduced at the beginning of the century (Kriegshauser et al. 2000) are now routinely used to estimate the resistivity anisotropy and formation dip in rock formations. These tools comprise three mutually orthogonal transmitter coils located at the same position and a similar configuration of receiver coils separated by an offset distance. Using this configuration, all nine couplings between each transmitter and receiver are measured, giving the most complete data set possible. One set of triaxial induction tool measurements includes nine orthogonal magnetic field components, these are the three orthogonal magnetic fields Hx, Hy, Hz separately excited by each of three mutually orthogonal magnetic sources Mx, My, Mz; the Hα component excited by Mβ source is denoted by Hαβ. The sensitivity of triaxial induction instruments to both vertical and horizontal resistivities of the formations has led to the wide use of these tools in transversely isotropic and general anisotropic media (Wang et al. 2009). A high number of distinct measurements allows for full 3-D inversion of a complex formation structure at an increased depth of investigation. Typical transmitter–receiver spacing for the current generation of DDR technology (e.g. Schlumberger's GeoSphere tool) ranges from 5 m to 35 m and matches the scales of many reservoir environments and overlaps with the surface-seismic resolution (Dupuis & Denichou 2015). A rule of thumb is that the effective depth of investigation is equal to the maximum transmitter–receiver spacing in high resistivity environments and is halved in conductive formations. Multiple frequencies that vary by about an order of magnitude provide different sensitivities and depths of investigation thus improving the information content and making the inversion more robust. The simulation of a triaxial induction tool moving in a long HA/HZ well is a computationally demanding task since a large number of transmitter positions with three emission directions at each position. For a typical production well, the total number of forward modelling tasks can be in the thousands per frequency. In order to make the 3-D DDR inversion computationally feasible, one requires a forward solver that can accurately and efficiently simulate the forward problem for multiple transmitter positions. If a single mesh can include all source–receiver positions in the wellbore, the discretization of eq. (8) results in a sparse linear system with multiple right-hand sides. Thus, we use a parallel direct solver to simulate all transmitter positions at a comparable cost to one matrix factorization. The matrix factors of the forward model allow us to compute the adjoint solution at a small extra cost at the end of the iteration. PRE-CONDITIONING A pre-conditioning for the NLCG method optimizes its convergence and provides a better detection of deeper targets by preventing their shielding by shallower objects. Due to the physics of the problem, the amplitude of the gradient is strongly influenced not only by geometric decay of the EM field, but also by galvanic attenuation in conductive media. Using a wide range of frequencies in inversion can partially alleviate this problem since both resolution and penetration depth depend on the frequency content of the data. When the number of different frequencies at which data are collected or frequency bandwidth are relatively small, quality of inversion can be greatly improved by a proper pre-conditioning scheme. Without such scheme, the model updates may focus in the near-wellbore zone and conductivity changes at large distances from the tool may not be determined correctly during the inversion process (even when the forward modelling could detect them). At each iteration step, the gradient for frequency f and source s is rescaled according to the following rule   \begin{equation}\widetilde {{\bf g}}(f,s) = {{\bf g}}(f,s) \cdot \exp \left( {\frac{{k(f) \cdot z(s)}}{{\delta (f)}}} \right).\end{equation} (9)This compensates the depth-attenuation of the diffusive EM signal (Zach et al. 2008). Similar depth weighting schemes have been used in inversion of magnetic and gravity data (Oldenburg et al. 1998). Here, z(s) is the vertical distance to the current transmitter position in the wellbore, δ(f) is the skin depth for a given frequency, and k(f) is a parameter based on the near-wellbore resistivity that adjusts the scaling factor. The effect of such exponential rescaling of the gradient is shown in Fig. 1. The influence of deeper zones of the model is enhanced which is essential for inversion of the DDR measurements. Figure 1. View largeDownload slide Magnitude of the exponential scaling function (top) for the wellbore trajectory (black line); data gradient for two equal targets located at different depths from the wellbore (middle) and the same gradient as in the middle panel but rescaled (bottom). Figure 1. View largeDownload slide Magnitude of the exponential scaling function (top) for the wellbore trajectory (black line); data gradient for two equal targets located at different depths from the wellbore (middle) and the same gradient as in the middle panel but rescaled (bottom). WOLFE CONDITIONS Once the search direction uk is calculated at iteration k, the step length αk is chosen to update the model (6). The step lengths are typically chosen large for the first iterations (resulting in large model updates), then we decrease their size as the minimum is localized. Determining an optimal value for αk further decreases the misfit at the current iteration and thus reduces the total number of iterations, and avoids extra forward simulations at each iteration. At least one extra modelling is required to evaluate the sensitivity of the misfit function to the current step length and search direction. To determine whether a step length is sufficient, we use the Wolfe conditions (Wolfe 1969). When these conditions are satisfied for a given step length αk, it is considered suitable for the next model update. Otherwise, the step length is modified and Wolfe conditions are re-checked. Our method uses a parabolic interpolation technique to determine the optimal step length through a reduced number of evaluations. The first and the second Wolfe conditions are given by (Wolfe 1969)   \begin{equation}\phi ({{{\bf m}}_k} + {\alpha _k}{{{\bf u}}_k}) \le \phi ({{{\bf m}}_k}) + {c_1}{\alpha _k}{{\bf u}}_k^T\nabla \phi ({{{\bf m}}_k}),\end{equation} (10)  \begin{equation}{{\bf u}}_k^T\nabla \phi ({{{\bf m}}_k} + {\alpha _k}{{{\bf u}}_k}) \ge {c_2}{{\bf u}}_k^T\nabla \phi ({{{\bf m}}_k}),\end{equation} (11)where c1 and c2 are constants between 0 and 1. c1 is usually chosen to be quite small (around 10−3–10−4), while c2 should be around 0.1 for the NLCG method (Nocedal & Wright 2006). The first Wolfe condition ensures that the step length αk decreases the misfit function sufficiently, and the second one is the curvature condition that verifies whether the slope has been reduced sufficiently. In practice, only the first Wolfe condition can be used (dos Santos & Pestana 2015) because checking both of them (or using stricter conditions, such as the strong Wolfe condition) requires the evaluation of the gradient of the misfit function for the next model ∇ϕ(mk + αkuk), thereby demanding an additional adjoint solve. When a direct solver is employed and the matrix is already factorized during the forward solve, this additional adjoint problem can be solved relatively cheaply. Other acceptance conditions can be used to evaluate the step length guesses as well (see e.g. the discussion in Hager & Zhang 2006). However, simpler conditions, such as decreasing the misfit function ϕ(mk + αkuk) < ϕ(mk) can result in accepting suboptimal step lengths, while more complicated rules introduce additional parameters to tune. The Wolfe conditions ensure a sufficient misfit reduction at every iteration. Once the step length is accepted, the forward modelling used to evaluate the first Wolfe condition is re-used in the next iteration requiring no extra computations. To improve the computational efficiency of the forward and inverse modelling, we employ a hybrid MPI/OpenMP parallelization scheme, which is well-suited for modern multicore processor architectures. The frequencies are processed independently on different computational nodes, so the parallel scalability is mainly limited by the size of the forward problems. In the numerical examples given below, we use one node per frequency when distributing the modelling tasks among the computational nodes. A shared-memory solver PARDISO (Schenk & Gartner 2004) is employed in this case. As the number of frequencies in inversion typically does not exceed 7–10, distributed-memory solvers can be used for inversion on a larger cluster. NUMERICAL RESULTS We illustrate the performance of the inversion scheme we describe above by simulating two synthetic examples that illustrate the main application conditions of DDR tools. The first example models reservoir monitoring during production seeking to define the dynamic location of the oil–water contact (OWC) and thus to estimate the thickness of an oil column and the well's production potential. In this case, the extended depths of investigation of DDR tools can detect early water breakthroughs and increase the field recovery factor by maximizing the production performance of HA/HZ wells. The second scenario, which we refer to as the exploration condition, aims to detect and delineate the hydrocarbon-bearing sandstones from low-resistivity formations and shales in the vicinity of the wellbore. Such detection is necessary for finding good-quality reservoir sands and drilling on target for real time geosteering. The latest generation of DDR tools can reveal subsurface bedding and fluid-contact details at distances of several tens of meters away from the wellbore. In the following examples, we examine this sensitivity in the inversion process. The optimal survey configurations that exhibit the highest sensitivity in both conditions were determined from prior modelling studies (Puzyrev & Torres-Verdín 2016). Different depths of investigation are better illuminated with different frequencies. This requires simultaneous inversion at several frequencies from the induction range (5–50 kHz) used by many commercial instruments. More frequencies improve the inversion quality while they increase the computational cost since we use a frequency domain forward solver. We choose the frequency intervals used in the inversion from the modelling studies, which mainly depend on the resistivity in the near-wellbore zone. In the monitoring condition, the transmitter is located in resistive zones and the goal is to determine the OWC (i.e. the boundary with the conductive media). In the exploration condition, resistivity of the near-wellbore zone is lower than that of the sandstone targets being located. Thus, the frequencies employed in this case are lower than those used in the monitoring condition. We estimate the resistivity of a fluid saturated rock using Archie's relations that relate it to the porosity, ϕ, the brine resistivity, Rw, and the water saturation, Sw, as follows:   \begin{equation}{R_{\,t}} = a \cdot {R_w} \cdot {\phi ^{ - m}} \cdot S_w^{ - n}.\end{equation} (12) We choose the following rock/fluid properties and Archie's parameters in the monitoring condition: Rw = 0.13 Ω m, ϕ = 0.25, and Sw = 0.15. We set the porosity exponent m and the saturation exponent n to 2, and Winsauer's factor a to 1. Thus, the resistivity contrasts in the examples below do not exceed two orders of magnitude. In the applications considered herein, the wells typically do not have steel casings which can seriously affect the EM measurements. Steel-cased wells can be modelled with finite differences using finer grids (Puzyrev et al. 2017). Possibly, a more efficient and robust way to simulate wells of arbitrary trajectories with the state-of-the-art simulation methods lies in the development of 3-D inversion based on unstructured meshes. MONITORING CONDITION Fig. 2 shows the 3-D synthetic reservoir model we use as a target model for inversion of the monitoring scenario. We choose the resistivities of the fluid saturated rocks according to Archie's equation (12) which vary in this example from approximately 1.2 to 90 Ω m. Fig. 3 shows 2-D vertical slices through the centre of the target and initial models. The target model has two zones where an early water breakthrough should be expected. The main goal of this monitoring exercise is to detect these zones and thus diminish the chance of early water breakthrough and increase the well's producing lifetime. Figure 2. View largeDownload slide Target resistivity model of the reservoir monitoring condition. A hydrocarbon reservoir is penetrated by an HA/HZ well. Figure 2. View largeDownload slide Target resistivity model of the reservoir monitoring condition. A hydrocarbon reservoir is penetrated by an HA/HZ well. Figure 3. View largeDownload slide XZ slice at the centre of the true monitoring model (top) and the starting model (middle); YZ slices at Zones 1 and 2 of the starting model (bottom). The well path is indicated with the black line. Figure 3. View largeDownload slide XZ slice at the centre of the true monitoring model (top) and the starting model (middle); YZ slices at Zones 1 and 2 of the starting model (bottom). The well path is indicated with the black line. As a starting/reference model, we use an almost flat OWC representing the reservoir before the onset of production (Fig. 3). In this model, the OWC (i.e. the resistive/conductive boundary under the wellbore path) is located 12–18 m true vertical depth (TVD) below the well, while during production the OWC raises to 7 and 5.5 m in Zones 1 and 2, respectively. The top of the reservoir is located 6–8 m TVD above the well path. Fig. 4 shows the inversion results obtained using two multicomponent receivers located at offsets of 20 and 40 m from the transmitter. We choose this receiver spacing since alternatives lower that 20 m would be less sensitive in these monitoring conditions. This is due to the initial position of the OWC relative to the well which was 12–18 m. We use seven frequencies logarithmically spaced from 5 to 25 kHz in this case. At lower frequencies, the skin depth is too large compared to the dimensions of the problem, while at higher ones the signal falls below the noise floor at the longest offsets. We add a 2 per cent Gaussian noise to the data at each frequency. The OWC in both zones is successfully detected with a slight overestimation of its TVD and underestimation of the water conductivity in Zone 2. The strongest OWC response is observed at the 40 m offset for the source positions directly above the water-rise zones. The inversion process is not able to detect the resistivity variations in the middle part of the model, most likely because of the shadowing of this region by the water-rise zones. Figure 4. View largeDownload slide Resistivity variations in the monitoring model estimated with the inversion using long offsets (20 and 40 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) water breakthrough zones. The white box shows the inversion area. Figure 4. View largeDownload slide Resistivity variations in the monitoring model estimated with the inversion using long offsets (20 and 40 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) water breakthrough zones. The white box shows the inversion area. Fig. 5 shows the monitoring model inverted using the same configuration but with larger source–receiver offsets (30 and 60 m). These results are slightly better than the ones shown in Fig. 4. The position and TVD of the OWC in the water breakthrough zones, as well as their resistivity, are more precise. Both inversions show that this configuration provides accurate estimations of the formation resistivity in the monitoring scenario. The triaxial induction measurements produce significantly more useful information than coaxial measurements, thus leading to better inversion results. Figure 5. View largeDownload slide Resistivity variations in the monitoring model estimated with the inversion using very long offsets (30 and 60 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) water breakthrough zones. The white box shows the inversion area. Figure 5. View largeDownload slide Resistivity variations in the monitoring model estimated with the inversion using very long offsets (30 and 60 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) water breakthrough zones. The white box shows the inversion area. The finest grid used in the modelling has 220 × 120 × 80 cells resulting in about six million unknowns in the forward problem. The number of model parameters to be determined by inversion is approximately 1.1 million. We analyse 180 positions of the triaxial induction source with a spacing of 0.6 m in the well, resulting into 540 modelling tasks. While these numbers would pose a serious challenge for iterative schemes, the main limitation of using a direct solver in 3-D problems is not the number of sources, but rather the size of the resulting linear system. The computational effort for forward and backward substitutions of 540 right-hand-side vectors is still small compared to the cost of the factorization of a matrix of this size (Puzyrev et al. 2016). The inversion depicted in Fig. 5 took about 30 hr of CPU time using one computational node equipped with 16 Intel Xeon cores per frequency. Fig. 6 shows the convergence plot of the NLCG method with and without the pre-conditioning scheme. The misfit values are normalized to the initial misfit. The convergence is fast at the beginning of the inversion, when we determine the main structures of the target model. The misfit curve flattens after 15–20 iterations as only small changes in the model are performed. The computational cost of one iteration with or without the pre-conditioner is roughly the same since the gradient rescaling (9) is a very cheap operation compared to the forward solution. Figure 6. View largeDownload slide Convergence of the NLCG method for the monitoring problem of Fig. 5. Pre-conditioned and non-pre-conditioned normalized misfits are shown as functions of iteration number. Figure 6. View largeDownload slide Convergence of the NLCG method for the monitoring problem of Fig. 5. Pre-conditioned and non-pre-conditioned normalized misfits are shown as functions of iteration number. EXPLORATION CONDITION Herein, we test the accuracy and detectability of our inversion algorithm in exploration conditions. Fig. 7 shows the exploration model that includes several small resistive bodies (100 Ω m) located in the shales of varying resistivity (6–8 Ω m) below the wellbore. The dimensions of the targets range from 12 to 20 m in the horizontal direction and from 3 to 5.5 m in depth. In the previous studies, we confirmed that DDR tools are sensitive to small targets located at these distances from the well trajectory. Hence, we expect the DDR inversion to correctly map and characterize these resistive targets (possible hydrocarbon reservoirs). Fig. 8 shows the true and reference explorations models discretized on a quite coarse grid with about 0.4 million of cells. As a reference model in this example, we use a model with homogeneous 6 Ω m resistivity below the wellbore. Figure 7. View largeDownload slide The target resistivity model of the exploration condition. Several resistive targets (100 Ω m) are located in the shales of varying resistivity (6–8 Ω m) below the wellbore. Figure 7. View largeDownload slide The target resistivity model of the exploration condition. Several resistive targets (100 Ω m) are located in the shales of varying resistivity (6–8 Ω m) below the wellbore. Figure 8. View largeDownload slide XZ slice at the centre of the true exploration model (top) and YZ slices at Targets 1 and 2 (middle). As a starting model, we use the model with homogeneous 6 Ω m resistivity below the wellbore (bottom). Figure 8. View largeDownload slide XZ slice at the centre of the true exploration model (top) and YZ slices at Targets 1 and 2 (middle). As a starting model, we use the model with homogeneous 6 Ω m resistivity below the wellbore (bottom). Figs 9 and 10 show the results from two inversions using short (10 and 20 m) and medium (15 and 30 m) transmitter–receiver offsets, respectively. Both figures show that Target 1 (located several metres below the wellbore) produces a much stronger response than the deeper resistor (Target 2) buried at 12–16 m. The presence of the latter is not visible at the short offsets (under 20 m) in this exploration condition. Inversion using longer offsets detects a resistive object and recovers accurately its size and depth, but not the horizontal position. The edge of Target 3 is visible in both configurations. These results confirm that the inversion method yields accurate estimates of electrical resistivity for the exploration condition. Three frequencies from the 5–12 kHz range and 100 positions of the source (spacing of 0.5 m) are used in these simulations. Such dense transmitter coverage is not strictly necessary in this case and one can employ only a subset of the transmitter–receiver positions for the inversion. The size of the forward modelling problems ranges from approximately 1.2 to 2.5 million of unknowns and the number of model parameters is approximately 0.4 million. The CPU time required is less than 3 hr using one computational node per frequency. Figure 9. View largeDownload slide Inversion results for the exploration model using short offsets (10 and 20 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) targets. The white box shows the inversion area. Figure 9. View largeDownload slide Inversion results for the exploration model using short offsets (10 and 20 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) targets. The white box shows the inversion area. Figure 10. View largeDownload slide Inversion results for the exploration model using medium offsets (15 and 30 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) targets. The white box shows the inversion area. Figure 10. View largeDownload slide Inversion results for the exploration model using medium offsets (15 and 30 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) targets. The white box shows the inversion area. DISCUSSION DDR tools provide multi-spacing and multi-frequency data sets and thus are no longer restricted to the detection of a single boundary with simplified inversion methods. These logging tools enable reliable estimation of spatial distributions of electrical resistivity several tens of metres away from the well. The interpretation of DDR data in complex multilayered structures should rely on full 3-D modelling and inversion. However, the high computational cost of 3-D numerical methods often constrains the development and usage of inversion-based interpretation. Our numerical results show that by making use of recent advances in linear algebra and parallel computing, the full 3-D inversion of DDR data can be performed in practical CPU times. The inversion of the fine-scale monitoring condition took about 30 hr of wall clock time using one computational node per frequency. The exploration test case was inverted by using smaller grids in 3 hr of CPU time. With appropriate computing resources, 3-D multisource inversion similar to our reservoir monitoring example can be performed within a few hours. A computationally efficient method based on unstructured grids would allow for further improvements in performance. Optimal transmitter–receiver configurations vary for different application conditions and frequencies. With the acquisition geometry being limited when compared to surface EM surveys, deeper zones are often shadowed by nearby resistivity inhomogeneities. Large source–receiver offsets (double of the target depth for the second receiver) provide advantageous inversion results compared to short offsets. The combination of resistivity measurements acquired with the DDR tools and seismic amplitude measurements will make possible the optimal steering of wells in geologically complex reservoirs. CONCLUSIONS We introduce and benchmark a new inversion method for 3-D DDR measurements in complex heterogeneous reservoirs penetrated by high-angle and horizontal wells. The key ingredients of our numerical scheme are the following. The minimization scheme employs the NLCG algorithm whose accuracy and performance is improved by invoking a pre-conditioning scheme. We use the Wolfe conditions to evaluate whether the step length is sufficient to provide near optimal convergence. Forward modelling of multicomponent triaxial induction data for multiple transmitter–receiver configurations uses a parallel direct solver. All components of the scheme, namely forward modelling, inversion algorithm, regularization, and pre-conditioner, are implemented in parallel and are well suited for execution on modern computer architectures. The workflow is validated on synthetic data sets that include monitoring and exploration borehole induction logging conditions. Numerical results verify that large transmitter–receiver offsets increase the depths of investigation to several tens of metres. Inversion of the monitoring example shows that the developed inversion-based method can determine the position of the oil-water contact in early water breakthrough zones, as well as the resistivity of the nearest formations. At the same time, the resistivity variations located at larger distances from the wellbore are not detected. These variations are most likely shadowed by water-rise zones. In the exploration conditions, the inversion algorithm provides accurate estimations of formation resistivity as well. The detection of deep (12–16 m away from the wellbore) targets requires offsets of 30 m to be reliably captured. These zones are often hidden by the resistivity variations in the vicinity of the well. These limitations of the acquisition geometry can be mitigated, for example, by using additional structural information available from other data sets, such as seismic amplitudes, to constrain the target formations. The synthetic data sets used in these simulations include 2 per cent Gaussian noise. Additional simulations with higher levels of noise verify the reliability of the inversion process. In particular, when some of the resistivity measurements are significantly noisier than the others (or missing at all). The forward modelling uses the finite difference method with a direct solver. For a multisource triaxial induction problem, this method enables a large reduction in the computational cost when compared to schemes based on iterative solvers. Future work will study the reliability and efficiency of the finite element method for this class of problems. Extending the inversion-based method to unstructured grids to reduce the number of unknowns in the resulting linear systems and thus to use solvers in a more efficient way is the subject of a forthcoming study. ACKNOWLEDGEMENTS This publication was made possible in part by the CSIRO Professorial Chair in Computational Geoscience at Curtin University and the Deep Earth Imaging Enterprise Future Science Platforms of the CSIRO. Additional support was provided by the National Priorities Research Program grant 7-1482-1-278 from the Qatar National Research Fund (a member of The Qatar Foundation) and by the European Union's Horizon 2020 Research and Innovation Program of the Marie Skłodowska-Curie grant agreement No. 644202. The J. Tinsley Oden Faculty Fellowship Research Program at the Institute for Computational Engineering and Sciences (ICES) of the University of Texas at Austin has partially supported the visits of VMC to ICES. Computational resources were provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia. The authors also acknowledge the University of Texas at Austin's Research Consortium on Formation Evaluation, jointly sponsored by Anadarko, Aramco, Baker-Hughes, BHP Billiton, BP, China Oilfield Services LTD., ConocoPhillips, Det Norske, ENI, ExxonMobil, Hess, Paradigm, Petrobras, Repsol, DEA, Schlumberger, Shell, Statoil, TOTAL, Wintershall and Woodside Petroleum Limited. The authors would like to thank Ute Weckmann, Colin Farquharson and one anonymous reviewer for their constructive suggestions, which helped to improve the manuscript. REFERENCES Abubakar A., Habashy T.M., Druskin V.L., Knizhnerman L., Alumbaugh D., 2008. 2.5 D forward and inverse modeling for interpreting low-frequency electromagnetic measurements, Geophysics , 73( 4), F165– F177. Google Scholar CrossRef Search ADS   Bakr S.A., Pardo D., Torres-Verdín C., 2017. Fast inversion of logging-while-drilling resistivity measurements acquired in multiple wells, Geophysics , 82( 3), E111– E120. Google Scholar CrossRef Search ADS   Bensdorp S., Petersen S.A., Olsen P.A., Antonsen F., van den Berg P.M., Fokkema J.T., 2016. An approximate 3D inversion method for inversion of single-well induction-logging responses, Geophysics , 81( 1), E43– E56. Google Scholar CrossRef Search ADS   Davydycheva S., Druskin V., Habashy T., 2003. An efficient finite-difference scheme for electromagnetic logging in 3D anisotropic inhomogeneous media, Geophysics , 68( 5), 1525– 1536. Google Scholar CrossRef Search ADS   dos Santos A.W., Pestana R.C., 2015. Time-domain multiscale full-waveform inversion using the rapid expansion method and efficient step-length estimation, Geophysics , 80( 4), R203– R216. Google Scholar CrossRef Search ADS   Dupuis C., Denichou J.M., 2015. Automatic inversion of deep-directional-resistivity measurements for well placement and reservoir description, Leading Edge , 34( 5), 504– 512. Google Scholar CrossRef Search ADS   Ezioba U., Denichou J.M., 2014. Mapping-while-drilling system improves well placement and field development, J. Pet. Technol. , 66( 08), 32– 35. Google Scholar CrossRef Search ADS   Farquharson C.G., Oldenburg D.W., 2004. A comparison of automatic techniques for estimating the regularization parameter in non-linear inverse problems, Geophys. J. Int. , 156( 3), 411– 425. Google Scholar CrossRef Search ADS   Hager W.W., Zhang H., 2006. A survey of nonlinear conjugate gradient methods, Pac. J. Optim. , 2( 1), 35– 58. Hou J., Mallan R.K., Torres-Verdín C., 2006. Finite-difference simulation of borehole EM measurements in 3D anisotropic media using coupled scalar-vector potentials, Geophysics , 71( 5), G225– G233. Google Scholar CrossRef Search ADS   Ijasan O., Torres-Verdín C., Preeg W.E., Rasmus J., Stockhausen E.J., 2014. Inversion-based interpretation of logging-while-drilling resistivity and nuclear measurements: field examples of application in high-angle and horizontal wells, Petrophysics , 55( 05), 374– 391. Jaysaval P., Shantsev D.V., de Ryhove S.D.L.K., Bratteland T., 2016. Fully anisotropic 3-D EM modelling on a Lebedev grid with a multigrid pre-conditioner, Geophys. J. Int. , 207( 3), 1554– 1572. Google Scholar CrossRef Search ADS   Kriegshauser B., Fanini O., Forgang S., Itskovich G., Rabinovich M., Tabarovsky L., Yu L., Epov M., 2000. A new multi-component induction logging tool to resolve anisotropic formations, in 40th Annual Logging Symposium , Society of Petrophysicists and Well Log Analysts. 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( 5), 499– 524. Google Scholar CrossRef Search ADS   Newman G.A., Alumbaugh D.L., 1997. Three-dimensional massively parallel electromagnetic inversion—I. Theory, Geophys. J. Int. , 128( 2), 345– 354. Google Scholar CrossRef Search ADS   Newman G.A., Alumbaugh D.L., 2002. Three-dimensional induction logging problems: Part 2—a finite-difference solution, Geophysics , 67( 2), 484– 491. Google Scholar CrossRef Search ADS   Newman G.A., Boggs P.T., 2004. Solution accelerators for large-scale three-dimensional electromagnetic inverse problems, Inverse Problems , 20( 6), S151– S170. Nocedal J., Wright S., 2006. Numerical Optimization , Springer Science & Business Media. Oldenburg D.W., Li Y., Farquharson C.G., Kowalczyk P., Aravanis T., King A., Zhang P., Watts A., 1998. Applications of geophysical inversions in mineral exploration, Leading Edge , 17( 4), 461– 465. Google Scholar CrossRef Search ADS   Pardo D., Demkowicz L., Torres-Verdín C., Paszynski M., 2006. Two-dimensional high-accuracy simulation of resistivity logging-while-drilling (LWD) measurements using a self-adaptive goal-oriented hp finite element method, SIAM J. Appl. Math. , 66( 6), 2085– 2106. Google Scholar CrossRef Search ADS   Pardo D., Matuszyk P.J., Torres‐Verdín C., Mora A., Muga I., Calo V.M., 2013. Influence of borehole-eccentred tools on wireline and logging‐while‐drilling sonic logging measurements, Geophys. Prospect. , 61( s1), 268– 283. Google Scholar CrossRef Search ADS   Plessix R.E., 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications, Geophys. J. Int. , 167( 2), 495– 503. Google Scholar CrossRef Search ADS   Polak E., Ribiere G., 1969. Note sur la convergence de méthodes de directions conjuguées, Rev. Fr. Inform. Rech. Oper., série rouge , 3( 1), 35– 43. Puzyrev V., Torres-Verdín C., 2016. Inversion-based interpretation of 3D borehole directional resistivity measurements acquired in high-angle and horizontal wells, in 2016 SEG Annual Meeting, SEG, Expanded Abstracts , pp. 724– 728. Puzyrev V., Koric S., Wilkin S., 2016. Evaluation of parallel direct sparse linear solvers in electromagnetic geophysical problems, Comput. Geosci. , 89, 79– 87. Google Scholar CrossRef Search ADS   Puzyrev V., Vilamajo E., Queralt P., Ledo J., Marcuello A., 2017. Three-dimensional modeling of the casing effect in onshore controlled-source electromagnetic surveys, Surv. Geophys. , 38( 2), 527– 545. Google Scholar CrossRef Search ADS   Rodi W., Mackie R.L., 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion, Geophysics , 66( 1), 174– 187. Google Scholar CrossRef Search ADS   Schenk O., Gartner K., 2004. Solving unsymmetric sparse systems of linear equations with PARDISO, Future Gener. Comput. Syst. , 20, 475– 487. Google Scholar CrossRef Search ADS   Seydoux J., Legendre E., Mirto E., Dupuis C., Denichou J.M., Bennett N., Kutiev G., Kuchenbecker M., Morriss C., Yang L., 2014. Full 3D deep directional resistivity measurements optimize well placement and provide reservoir-scale imaging while drilling, 55th Annual Logging Symposium , Society of Petrophysicists and Well-Log Analysts. Shantsev D.V., Jaysaval P., de Ryhove S.D.L.K., Amestoy P.R., Buttari A., L’Excellent J.Y., Mary T., 2017. Large-scale 3-D EM modelling with a Block Low-Rank multifrontal direct solver, Geophys. J. Int. , 209( 3), 1558– 1571. Google Scholar CrossRef Search ADS   Thiel M., Dupuis C., Omeragic D., 2016. Fast 2D inversion for enhanced real-time reservoir mapping while drilling, in SPE Annual Technical Conference and Exhibition , Society of Petroleum Engineers. Wang G.L., Torres-Verdín C., Gianzero S., 2009. Fast simulation of triaxial borehole induction measurements acquired in axially symmetrical and transversely isotropic media, Geophysics , 74( 6), E233– E249. Google Scholar CrossRef Search ADS   Ward S.H., Hohmann G.W., 1988. Electromagnetic theory for geophysical applications, in Electromagnetic Methods in Applied Geophysics , Vol. 1, pp. 131– 311, ed. Nabighian M.N., Society of Exploration Geophysicists. Wiese T., Greenhalgh S., Zhou B., Greenhalgh M., Marescot L., 2015. Resistivity inversion in 2-D anisotropic media: numerical experiments, Geophys. J. Int. , 201( 1), 247– 266. Google Scholar CrossRef Search ADS   Wolfe P., 1969. Convergence conditions for ascent methods, SIAM Rev. , 11( 2), 226– 235. Google Scholar CrossRef Search ADS   Zach J.J., Bjørke A.K., Støren T., Maaø F., 2008. 3D inversion of marine CSEM data using a fast finite-difference time-domain forward code and approximate Hessian-based optimization, in 2008 SEG Annual Meeting, SEG, Expanded Abstracts , pp. 614– 618. Zhdanov M.S., 2009. New advances in regularized inversion of gravity and electromagnetic data, Geophys. Prospect. , 57( 4), 463– 478. Google Scholar CrossRef Search ADS   © The Author(s) 2018. 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

Interpretation of deep directional resistivity measurements acquired in high-angle and horizontal wells using 3-D inversion

Loading next page...
 
/lp/ou_press/interpretation-of-deep-directional-resistivity-measurements-acquired-rhn6f1wwdw
Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggy047
Publisher site
See Article on Publisher Site

Abstract

SUMMARY The interpretation of resistivity measurements acquired in high-angle and horizontal wells is a critical technical problem in formation evaluation. We develop an efficient parallel 3-D inversion method to estimate the spatial distribution of electrical resistivity in the neighbourhood of a well from deep directional electromagnetic induction measurements. The methodology places no restriction on the spatial distribution of the electrical resistivity around arbitrary well trajectories. The fast forward modelling of triaxial induction measurements performed with multiple transmitter–receiver configurations employs a parallel direct solver. The inversion uses a pre-conditioned gradient-based method whose accuracy is improved using the Wolfe conditions to estimate optimal step lengths at each iteration. The large transmitter–receiver offsets, used in the latest generation of commercial directional resistivity tools, improve the depth of investigation to over 30 m from the wellbore. Several challenging synthetic examples confirm the feasibility of the full 3-D inversion-based interpretations for these distances, hence enabling the integration of resistivity measurements with seismic amplitude data to improve the forecast of the petrophysical and fluid properties. Employing parallel direct solvers for the triaxial induction problems allows for large reductions in computational effort, thereby opening the possibility to invert multiposition 3-D data in practical CPU times. Controlled source electromagnetics (CSEM), Electromagnetic theory, Inverse theory, Numerical approximations and analysis, Numerical modelling INTRODUCTION With hydrocarbon exploration and production continuously moving into areas of increasing geological and fluid complexity, the oil and gas industry is adopting new technologies to improve the efficiency of the well placement. One of the most important properties of geological formations is their electrical resistivity, which has high sensitivity to both porosity and hydrocarbon saturation. Given the limitations of current electromagnetic (EM) technologies such as logging-while-drilling (LWD), new monitoring and exploration methods were recently introduced. A new generation of deep directional resistivity (DDR) tools provides accurate estimations of various petrophysical properties of hydrocarbon-producing rock formations and has an improved depth of investigation of more than 30 m from the wellbore (Ezioba & Denichou 2014; Seydoux et al. 2014; Dupuis & Denichou 2015). This technology, often referred to as reservoir mapping-while-drilling, several times exceeds the spatial coverage of conventional LWD tools and offers an unprecedented opportunity to map the near-wellbore reservoir structure. In turn, the improved ability to navigate well trajectories with respect to geological and fluid boundaries using real-time information allows operators to reduce drilling risk, maximize reservoir exposure, and increase the overall production potential. With recent advances in drilling technology, long high-angle (HA) and horizontal (HZ) wells have become common in the oil and gas industry (Ijasan et al. 2014). The accuracy of landing HZ wells in a reservoir embedded in a complex geology often determines the overall financial viability. Many of the conventional LWD modelling methods were developed for vertical wells and are not suitable for interpretation of DDR measurements in HA/HZ wells. The real-time processing of DDR measurements in complex geologic media requires a full 3-D inversion method. An efficient inversion will open the possibility of estimating 3-D distributions of electrical conductivity. However, one of the main practical limitations of this approach is its high computational cost. Because the DDR tool is moving in the well, modelling needs to consider a large number of transmitter and receiver positions (in the range of thousands of samples when the position and the three independent emission directions are considered) for a long HA/HZ well. The inversion process, in turn, involves a high number of forward modelling operations. Thus, the total computational cost is mainly influenced by the cost of the forward modelling. Multiple 3-D modelling methods of multicomponent borehole EM induction logging measurements have been developed during the past decades (e.g. Newman & Alumbaugh 2002; Davydycheva et al. 2003; Hou et al. 2006; among others). However, the complexity of 3-D models dramatically increases the computational requirements for modelling and hence for inversion. The size of these problems has often been considered excessively large to be solved in real time for practical applications, thereby preventing the widespread use of full 3-D methods for the interpretation of borehole induction measurements. Many of the published borehole induction modelling methods take advantage of the axial symmetry of formation geometry and employ 2-D or 2.5D modelling and inversion algorithms (Pardo et al. 2006; Abubakar et al. 2008; Wang et al. 2009; Pardo et al. 2013; Wiese et al. 2015). Other approximations have been employed as well to limit the large number of model parameters and to mitigate the computational cost of inversion. For example, recently, Thiel et al. (2016) documented a fast 2-D inversion that utilized pointwise 1-D forward simulations for real-time interpretation of DDR measurements for MWD. Bensdorp et al. (2016) described the inversion for 3-D conductivity distributions based on single spherical scatterer approximations. Bakr et al. (2017) used 1.5D simulations for fast inversion of 3-D transversely isotropic formations approximated by planar layers. With the development of high-performance computing methods and the increase in the availability of computer resources, 3-D inversion of geophysical data sets no longer poses an insurmountable computational challenge. While still requiring significant computational resources for large-scale exploration problems, inversion of local data sets can be performed on a workstation or a small computer cluster in practical CPU times, that is, hours or even minutes. In particular, recent advances in parallel direct solvers allow us to improve the efficiency of simulations of multiposition and multicomponent triaxial induction measurements in complex reservoir formations with extreme resistivity contrasts (Puzyrev & Torres-Verdín 2016). The main limitation of these solution techniques is their large computational and memory demand; however, many state-of-the-art EM modelling problems involve up to several million unknowns and thus can be efficiently solved with direct solvers. Recent work by Shantsev et al. (2017) presents a large-scale 3-D controlled-source EM problem with more than twenty million unknowns solved with a block low-rank multifrontal direct solver. In this paper, we describe a full 3-D parallel inversion method to interpret multicomponent DDR measurements acquired in HA/HZ and apply it to several synthetic models. The remainder of this paper is organized as follows. In the next section, we summarize the fundamentals of the inverse problem. Then, we discuss the numerical approach used for forward modelling of multicomponent DDR measurements and specific details of the minimization algorithm such as pre-conditioning and optimal step length estimation. In the Examples section, we study the sensitivity of DDR tools to determine the oil-water contact during production (monitoring problem) and to locate favourable sandstone targets (exploration problem). Finally, the last section summarizes the conclusions and implications of our study. THEORY The inverse problem is formulated as a regularized nonlinear minimization problem with the following quadratic cost function:   \begin{equation}\phi {{\bf (m)}} = \frac{1}{2}\left\| {{{\bf D}}\left( {{{\bf F}}({{\bf m}}) - {{{\bf d}}_{{{\bf obs}}}}} \right)} \right\|_2^2\,\, + \,\,\frac{1}{2}\lambda \,{{\bf R}}({{\bf m}}),\end{equation} (1)where the first term in the right-hand side is the L2-norm based data misfit, dobs is the vector of measurements (observations), F(m) is the forward problem mapping, D is a diagonal data weighting matrix, R (m) is the stabilizing functional, and λ is the Lagrange multiplier. In this study, we use the pre-conditioned nonlinear conjugate gradient (NLCG) method to minimize the misfit between the measurements and their numerical simulations. This algorithm is efficient for parallel computers and has a small memory footprint compared to second-order minimization schemes (Rodi & Mackie 2001; Nocedal & Wright 2006). The stabilizing functional ensures the well-posedness of the nonlinear inverse problem in the presence of noise and/or inadequate measurements. Various regularization methods exist in the literature (Farquharson & Oldenburg 2004; Zhdanov 2009). In the examples described below, we use the Tikhonov minimum-norm regularization, given by   \begin{equation}{{\bf R}}\,({{\bf m}}) = \left\| {{{\bf m}} - {{{\bf m}}_{{{\bf ref}}}}} \right\|_2^2,\end{equation} (2)where vectors m and mref designate the parameters of the current and the reference models, respectively. The model parameters are commonly chosen as the logarithm of the conductivity σ, thus decreasing large gradient variations   \begin{equation}{{\bf m}} = \ln {{\boldsymbol \sigma }}, \quad {{{\bf m}}_{{{\bf ref}}}} = \ln {{{\boldsymbol \sigma }}_{{{\bf ref}}}}.\end{equation} (3) The reference model contains a priori information about subsurface properties, which is especially important for spatially complex geologies. Because one can often choose a good reference model mref and expect sharp resistivity contrasts in the conditions under analysis, regularization (2) has been found more suitable than, for example, maximum smoothness methods based on L2 or L1 norms. The inversion algorithm estimates formation resistivity by progressively minimizing the misfit function (1). The data gradient (i.e. the derivative of the data term in eq. (1) with respect to the model parameters, m) is given by (Newman & Boggs 2004)   \begin{equation}{{{\bf g}}_{{\bf d}}} = \nabla {\phi _{{\bf d}}}({{\bf m}}) = - {\mathop{\rm Re}\nolimits} \left[ {{{\left( {{{\bf DJ}}} \right)}^T}{{\left[ {{{\bf D}}\left( {{{\bf F}}({{\bf m}}) - {{{\bf d}}_{{{\bf obs}}}}} \right)} \right]}^*}} \right],\end{equation} (4)where * stands for complex conjugation. Deriving a computationally efficient form of eq. (4) is critical for a robust inverse solution (McGillivray & Oldenburg 1990). The adjoint method to calculate the data gradient has been used previously in EM (Newman & Alumbaugh 1997) and seismic full waveform inversion (Plessix 2006). This procedure allows us to avoid the explicit calculation and storage of the Jacobian matrix, J. When a direct solver is employed, the adjoint field can be calculated at a relatively low cost since the EM inverse problem is self-adjoint. At each iteration step, after the frequency-dependent forward problem F(m) is solved, we used the system matrix factorization to calculate the misfit gradient using the adjoint method. Thus, the expensive factorization is performed once for each frequency at the beginning of each iteration and then it is applied in the forward and the adjoint solves throughout this iteration. The data misfits are added over all frequencies and pre-conditioned to calculate the conjugate gradient and the search direction uk + 1 given by   \begin{equation}{{{\bf u}}_{{{\bf k + 1}}}} = - {{{\bf g}}_{{{\bf k + 1}}}} + {\beta _k}{{{\bf u}}_{{\bf k}}},\quad {{{\bf u}}_{{\bf 0}}} = - {{{\bf g}}_{{\bf 0}}}.\end{equation} (5)Here, the scalar βk is the conjugate gradient update parameter. Different conjugate gradient methods correspond to different choices for βk. We refer the reader to Hager & Zhang (2006) for a detailed chronological description of the update parameters. For large-scale problems, we prefer a βk that does not require the evaluation of the Hessian matrix. We use the standard Polak–Ribiere formula in our geophysical inversion algorithm (Polak & Ribiere 1969). Finally, the optimal step length αk is estimated as described below and the model is updated as   \begin{equation}{{{\bf m}}_{{{\bf k + 1}}}} = {{{\bf m}}_{{\bf k}}} + {\alpha _k}{{{\bf u}}_{{\bf k}}}.\end{equation} (6) This iterative process starts from the starting model m0 which can be (but not necessarily is) chosen the same as the reference model. FORWARD MODELLING To obtain the simulated measurements d = F(m), we solve the forward problem F(m) that is described by the time-harmonic diffusive Maxwell equations   \begin{equation}\nabla \times {{\bf E}} = - \,i\omega {\mu _0}{{\bf H}} - {{{\bf M}}_{{\bf s}}},\quad \nabla \times {{\bf H}} = \bar{\sigma }{{\bf E}} + {{{\bf J}}_{{\bf s}}},\end{equation} (7)where ω is angular frequency, $$\bar{\sigma }({{\bf r}})$$ is electric conductivity tensor, μ0 is magnetic permeability of free space, Js and Ms are the electric and magnetic sources, respectively (Ward & Hohmann 1988). Taking the curl of the first equation in eq. (7) and substituting it in the second one gives the curl-curl equation for the electric field   \begin{equation}\nabla \times \nabla \times {{\bf E}} + i\omega {\mu _0}\bar{\sigma }{{\bf E}} = - i\omega {\mu _0}{{{\bf J}}_{{\bf s}}} - \nabla \times {{{\bf M}}_{{\bf s}}}.\end{equation} (8) Reliable inversion and interpretation of multicomponent induction data in complex geological environments requires accurate and efficient forward modelling. In our implementation, we approximate eq. (8) by finite differences based on a staggered-grid method that is commonly applied in geophysical inversion (e.g. Newman & Alumbaugh 2002; Hou et al. 2006; Abubakar et al. 2008). The scattered-field formulation of eq. (8) is used to avoid excessive meshing near the sources. Our method handles large-scale problems in parallel and was successfully tested on complex isotropic and transversely isotropic reservoir models with extreme resistivity contrasts. The results were in good agreement with those obtained from a semi-analytic 1-D approach, thus confirming the accuracy and reliability of the simulation method. Arbitrary anisotropic media could be accurately simulated by using Lebedev grids (Davydycheva et al. 2003; Jaysaval et al. 2016). Modern DDR tools provide multicomponent measurements at various frequencies and transmitter–receiver offsets to estimate the complexity of the surrounding formations with better spatial resolution. Triaxial induction tools introduced at the beginning of the century (Kriegshauser et al. 2000) are now routinely used to estimate the resistivity anisotropy and formation dip in rock formations. These tools comprise three mutually orthogonal transmitter coils located at the same position and a similar configuration of receiver coils separated by an offset distance. Using this configuration, all nine couplings between each transmitter and receiver are measured, giving the most complete data set possible. One set of triaxial induction tool measurements includes nine orthogonal magnetic field components, these are the three orthogonal magnetic fields Hx, Hy, Hz separately excited by each of three mutually orthogonal magnetic sources Mx, My, Mz; the Hα component excited by Mβ source is denoted by Hαβ. The sensitivity of triaxial induction instruments to both vertical and horizontal resistivities of the formations has led to the wide use of these tools in transversely isotropic and general anisotropic media (Wang et al. 2009). A high number of distinct measurements allows for full 3-D inversion of a complex formation structure at an increased depth of investigation. Typical transmitter–receiver spacing for the current generation of DDR technology (e.g. Schlumberger's GeoSphere tool) ranges from 5 m to 35 m and matches the scales of many reservoir environments and overlaps with the surface-seismic resolution (Dupuis & Denichou 2015). A rule of thumb is that the effective depth of investigation is equal to the maximum transmitter–receiver spacing in high resistivity environments and is halved in conductive formations. Multiple frequencies that vary by about an order of magnitude provide different sensitivities and depths of investigation thus improving the information content and making the inversion more robust. The simulation of a triaxial induction tool moving in a long HA/HZ well is a computationally demanding task since a large number of transmitter positions with three emission directions at each position. For a typical production well, the total number of forward modelling tasks can be in the thousands per frequency. In order to make the 3-D DDR inversion computationally feasible, one requires a forward solver that can accurately and efficiently simulate the forward problem for multiple transmitter positions. If a single mesh can include all source–receiver positions in the wellbore, the discretization of eq. (8) results in a sparse linear system with multiple right-hand sides. Thus, we use a parallel direct solver to simulate all transmitter positions at a comparable cost to one matrix factorization. The matrix factors of the forward model allow us to compute the adjoint solution at a small extra cost at the end of the iteration. PRE-CONDITIONING A pre-conditioning for the NLCG method optimizes its convergence and provides a better detection of deeper targets by preventing their shielding by shallower objects. Due to the physics of the problem, the amplitude of the gradient is strongly influenced not only by geometric decay of the EM field, but also by galvanic attenuation in conductive media. Using a wide range of frequencies in inversion can partially alleviate this problem since both resolution and penetration depth depend on the frequency content of the data. When the number of different frequencies at which data are collected or frequency bandwidth are relatively small, quality of inversion can be greatly improved by a proper pre-conditioning scheme. Without such scheme, the model updates may focus in the near-wellbore zone and conductivity changes at large distances from the tool may not be determined correctly during the inversion process (even when the forward modelling could detect them). At each iteration step, the gradient for frequency f and source s is rescaled according to the following rule   \begin{equation}\widetilde {{\bf g}}(f,s) = {{\bf g}}(f,s) \cdot \exp \left( {\frac{{k(f) \cdot z(s)}}{{\delta (f)}}} \right).\end{equation} (9)This compensates the depth-attenuation of the diffusive EM signal (Zach et al. 2008). Similar depth weighting schemes have been used in inversion of magnetic and gravity data (Oldenburg et al. 1998). Here, z(s) is the vertical distance to the current transmitter position in the wellbore, δ(f) is the skin depth for a given frequency, and k(f) is a parameter based on the near-wellbore resistivity that adjusts the scaling factor. The effect of such exponential rescaling of the gradient is shown in Fig. 1. The influence of deeper zones of the model is enhanced which is essential for inversion of the DDR measurements. Figure 1. View largeDownload slide Magnitude of the exponential scaling function (top) for the wellbore trajectory (black line); data gradient for two equal targets located at different depths from the wellbore (middle) and the same gradient as in the middle panel but rescaled (bottom). Figure 1. View largeDownload slide Magnitude of the exponential scaling function (top) for the wellbore trajectory (black line); data gradient for two equal targets located at different depths from the wellbore (middle) and the same gradient as in the middle panel but rescaled (bottom). WOLFE CONDITIONS Once the search direction uk is calculated at iteration k, the step length αk is chosen to update the model (6). The step lengths are typically chosen large for the first iterations (resulting in large model updates), then we decrease their size as the minimum is localized. Determining an optimal value for αk further decreases the misfit at the current iteration and thus reduces the total number of iterations, and avoids extra forward simulations at each iteration. At least one extra modelling is required to evaluate the sensitivity of the misfit function to the current step length and search direction. To determine whether a step length is sufficient, we use the Wolfe conditions (Wolfe 1969). When these conditions are satisfied for a given step length αk, it is considered suitable for the next model update. Otherwise, the step length is modified and Wolfe conditions are re-checked. Our method uses a parabolic interpolation technique to determine the optimal step length through a reduced number of evaluations. The first and the second Wolfe conditions are given by (Wolfe 1969)   \begin{equation}\phi ({{{\bf m}}_k} + {\alpha _k}{{{\bf u}}_k}) \le \phi ({{{\bf m}}_k}) + {c_1}{\alpha _k}{{\bf u}}_k^T\nabla \phi ({{{\bf m}}_k}),\end{equation} (10)  \begin{equation}{{\bf u}}_k^T\nabla \phi ({{{\bf m}}_k} + {\alpha _k}{{{\bf u}}_k}) \ge {c_2}{{\bf u}}_k^T\nabla \phi ({{{\bf m}}_k}),\end{equation} (11)where c1 and c2 are constants between 0 and 1. c1 is usually chosen to be quite small (around 10−3–10−4), while c2 should be around 0.1 for the NLCG method (Nocedal & Wright 2006). The first Wolfe condition ensures that the step length αk decreases the misfit function sufficiently, and the second one is the curvature condition that verifies whether the slope has been reduced sufficiently. In practice, only the first Wolfe condition can be used (dos Santos & Pestana 2015) because checking both of them (or using stricter conditions, such as the strong Wolfe condition) requires the evaluation of the gradient of the misfit function for the next model ∇ϕ(mk + αkuk), thereby demanding an additional adjoint solve. When a direct solver is employed and the matrix is already factorized during the forward solve, this additional adjoint problem can be solved relatively cheaply. Other acceptance conditions can be used to evaluate the step length guesses as well (see e.g. the discussion in Hager & Zhang 2006). However, simpler conditions, such as decreasing the misfit function ϕ(mk + αkuk) < ϕ(mk) can result in accepting suboptimal step lengths, while more complicated rules introduce additional parameters to tune. The Wolfe conditions ensure a sufficient misfit reduction at every iteration. Once the step length is accepted, the forward modelling used to evaluate the first Wolfe condition is re-used in the next iteration requiring no extra computations. To improve the computational efficiency of the forward and inverse modelling, we employ a hybrid MPI/OpenMP parallelization scheme, which is well-suited for modern multicore processor architectures. The frequencies are processed independently on different computational nodes, so the parallel scalability is mainly limited by the size of the forward problems. In the numerical examples given below, we use one node per frequency when distributing the modelling tasks among the computational nodes. A shared-memory solver PARDISO (Schenk & Gartner 2004) is employed in this case. As the number of frequencies in inversion typically does not exceed 7–10, distributed-memory solvers can be used for inversion on a larger cluster. NUMERICAL RESULTS We illustrate the performance of the inversion scheme we describe above by simulating two synthetic examples that illustrate the main application conditions of DDR tools. The first example models reservoir monitoring during production seeking to define the dynamic location of the oil–water contact (OWC) and thus to estimate the thickness of an oil column and the well's production potential. In this case, the extended depths of investigation of DDR tools can detect early water breakthroughs and increase the field recovery factor by maximizing the production performance of HA/HZ wells. The second scenario, which we refer to as the exploration condition, aims to detect and delineate the hydrocarbon-bearing sandstones from low-resistivity formations and shales in the vicinity of the wellbore. Such detection is necessary for finding good-quality reservoir sands and drilling on target for real time geosteering. The latest generation of DDR tools can reveal subsurface bedding and fluid-contact details at distances of several tens of meters away from the wellbore. In the following examples, we examine this sensitivity in the inversion process. The optimal survey configurations that exhibit the highest sensitivity in both conditions were determined from prior modelling studies (Puzyrev & Torres-Verdín 2016). Different depths of investigation are better illuminated with different frequencies. This requires simultaneous inversion at several frequencies from the induction range (5–50 kHz) used by many commercial instruments. More frequencies improve the inversion quality while they increase the computational cost since we use a frequency domain forward solver. We choose the frequency intervals used in the inversion from the modelling studies, which mainly depend on the resistivity in the near-wellbore zone. In the monitoring condition, the transmitter is located in resistive zones and the goal is to determine the OWC (i.e. the boundary with the conductive media). In the exploration condition, resistivity of the near-wellbore zone is lower than that of the sandstone targets being located. Thus, the frequencies employed in this case are lower than those used in the monitoring condition. We estimate the resistivity of a fluid saturated rock using Archie's relations that relate it to the porosity, ϕ, the brine resistivity, Rw, and the water saturation, Sw, as follows:   \begin{equation}{R_{\,t}} = a \cdot {R_w} \cdot {\phi ^{ - m}} \cdot S_w^{ - n}.\end{equation} (12) We choose the following rock/fluid properties and Archie's parameters in the monitoring condition: Rw = 0.13 Ω m, ϕ = 0.25, and Sw = 0.15. We set the porosity exponent m and the saturation exponent n to 2, and Winsauer's factor a to 1. Thus, the resistivity contrasts in the examples below do not exceed two orders of magnitude. In the applications considered herein, the wells typically do not have steel casings which can seriously affect the EM measurements. Steel-cased wells can be modelled with finite differences using finer grids (Puzyrev et al. 2017). Possibly, a more efficient and robust way to simulate wells of arbitrary trajectories with the state-of-the-art simulation methods lies in the development of 3-D inversion based on unstructured meshes. MONITORING CONDITION Fig. 2 shows the 3-D synthetic reservoir model we use as a target model for inversion of the monitoring scenario. We choose the resistivities of the fluid saturated rocks according to Archie's equation (12) which vary in this example from approximately 1.2 to 90 Ω m. Fig. 3 shows 2-D vertical slices through the centre of the target and initial models. The target model has two zones where an early water breakthrough should be expected. The main goal of this monitoring exercise is to detect these zones and thus diminish the chance of early water breakthrough and increase the well's producing lifetime. Figure 2. View largeDownload slide Target resistivity model of the reservoir monitoring condition. A hydrocarbon reservoir is penetrated by an HA/HZ well. Figure 2. View largeDownload slide Target resistivity model of the reservoir monitoring condition. A hydrocarbon reservoir is penetrated by an HA/HZ well. Figure 3. View largeDownload slide XZ slice at the centre of the true monitoring model (top) and the starting model (middle); YZ slices at Zones 1 and 2 of the starting model (bottom). The well path is indicated with the black line. Figure 3. View largeDownload slide XZ slice at the centre of the true monitoring model (top) and the starting model (middle); YZ slices at Zones 1 and 2 of the starting model (bottom). The well path is indicated with the black line. As a starting/reference model, we use an almost flat OWC representing the reservoir before the onset of production (Fig. 3). In this model, the OWC (i.e. the resistive/conductive boundary under the wellbore path) is located 12–18 m true vertical depth (TVD) below the well, while during production the OWC raises to 7 and 5.5 m in Zones 1 and 2, respectively. The top of the reservoir is located 6–8 m TVD above the well path. Fig. 4 shows the inversion results obtained using two multicomponent receivers located at offsets of 20 and 40 m from the transmitter. We choose this receiver spacing since alternatives lower that 20 m would be less sensitive in these monitoring conditions. This is due to the initial position of the OWC relative to the well which was 12–18 m. We use seven frequencies logarithmically spaced from 5 to 25 kHz in this case. At lower frequencies, the skin depth is too large compared to the dimensions of the problem, while at higher ones the signal falls below the noise floor at the longest offsets. We add a 2 per cent Gaussian noise to the data at each frequency. The OWC in both zones is successfully detected with a slight overestimation of its TVD and underestimation of the water conductivity in Zone 2. The strongest OWC response is observed at the 40 m offset for the source positions directly above the water-rise zones. The inversion process is not able to detect the resistivity variations in the middle part of the model, most likely because of the shadowing of this region by the water-rise zones. Figure 4. View largeDownload slide Resistivity variations in the monitoring model estimated with the inversion using long offsets (20 and 40 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) water breakthrough zones. The white box shows the inversion area. Figure 4. View largeDownload slide Resistivity variations in the monitoring model estimated with the inversion using long offsets (20 and 40 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) water breakthrough zones. The white box shows the inversion area. Fig. 5 shows the monitoring model inverted using the same configuration but with larger source–receiver offsets (30 and 60 m). These results are slightly better than the ones shown in Fig. 4. The position and TVD of the OWC in the water breakthrough zones, as well as their resistivity, are more precise. Both inversions show that this configuration provides accurate estimations of the formation resistivity in the monitoring scenario. The triaxial induction measurements produce significantly more useful information than coaxial measurements, thus leading to better inversion results. Figure 5. View largeDownload slide Resistivity variations in the monitoring model estimated with the inversion using very long offsets (30 and 60 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) water breakthrough zones. The white box shows the inversion area. Figure 5. View largeDownload slide Resistivity variations in the monitoring model estimated with the inversion using very long offsets (30 and 60 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) water breakthrough zones. The white box shows the inversion area. The finest grid used in the modelling has 220 × 120 × 80 cells resulting in about six million unknowns in the forward problem. The number of model parameters to be determined by inversion is approximately 1.1 million. We analyse 180 positions of the triaxial induction source with a spacing of 0.6 m in the well, resulting into 540 modelling tasks. While these numbers would pose a serious challenge for iterative schemes, the main limitation of using a direct solver in 3-D problems is not the number of sources, but rather the size of the resulting linear system. The computational effort for forward and backward substitutions of 540 right-hand-side vectors is still small compared to the cost of the factorization of a matrix of this size (Puzyrev et al. 2016). The inversion depicted in Fig. 5 took about 30 hr of CPU time using one computational node equipped with 16 Intel Xeon cores per frequency. Fig. 6 shows the convergence plot of the NLCG method with and without the pre-conditioning scheme. The misfit values are normalized to the initial misfit. The convergence is fast at the beginning of the inversion, when we determine the main structures of the target model. The misfit curve flattens after 15–20 iterations as only small changes in the model are performed. The computational cost of one iteration with or without the pre-conditioner is roughly the same since the gradient rescaling (9) is a very cheap operation compared to the forward solution. Figure 6. View largeDownload slide Convergence of the NLCG method for the monitoring problem of Fig. 5. Pre-conditioned and non-pre-conditioned normalized misfits are shown as functions of iteration number. Figure 6. View largeDownload slide Convergence of the NLCG method for the monitoring problem of Fig. 5. Pre-conditioned and non-pre-conditioned normalized misfits are shown as functions of iteration number. EXPLORATION CONDITION Herein, we test the accuracy and detectability of our inversion algorithm in exploration conditions. Fig. 7 shows the exploration model that includes several small resistive bodies (100 Ω m) located in the shales of varying resistivity (6–8 Ω m) below the wellbore. The dimensions of the targets range from 12 to 20 m in the horizontal direction and from 3 to 5.5 m in depth. In the previous studies, we confirmed that DDR tools are sensitive to small targets located at these distances from the well trajectory. Hence, we expect the DDR inversion to correctly map and characterize these resistive targets (possible hydrocarbon reservoirs). Fig. 8 shows the true and reference explorations models discretized on a quite coarse grid with about 0.4 million of cells. As a reference model in this example, we use a model with homogeneous 6 Ω m resistivity below the wellbore. Figure 7. View largeDownload slide The target resistivity model of the exploration condition. Several resistive targets (100 Ω m) are located in the shales of varying resistivity (6–8 Ω m) below the wellbore. Figure 7. View largeDownload slide The target resistivity model of the exploration condition. Several resistive targets (100 Ω m) are located in the shales of varying resistivity (6–8 Ω m) below the wellbore. Figure 8. View largeDownload slide XZ slice at the centre of the true exploration model (top) and YZ slices at Targets 1 and 2 (middle). As a starting model, we use the model with homogeneous 6 Ω m resistivity below the wellbore (bottom). Figure 8. View largeDownload slide XZ slice at the centre of the true exploration model (top) and YZ slices at Targets 1 and 2 (middle). As a starting model, we use the model with homogeneous 6 Ω m resistivity below the wellbore (bottom). Figs 9 and 10 show the results from two inversions using short (10 and 20 m) and medium (15 and 30 m) transmitter–receiver offsets, respectively. Both figures show that Target 1 (located several metres below the wellbore) produces a much stronger response than the deeper resistor (Target 2) buried at 12–16 m. The presence of the latter is not visible at the short offsets (under 20 m) in this exploration condition. Inversion using longer offsets detects a resistive object and recovers accurately its size and depth, but not the horizontal position. The edge of Target 3 is visible in both configurations. These results confirm that the inversion method yields accurate estimates of electrical resistivity for the exploration condition. Three frequencies from the 5–12 kHz range and 100 positions of the source (spacing of 0.5 m) are used in these simulations. Such dense transmitter coverage is not strictly necessary in this case and one can employ only a subset of the transmitter–receiver positions for the inversion. The size of the forward modelling problems ranges from approximately 1.2 to 2.5 million of unknowns and the number of model parameters is approximately 0.4 million. The CPU time required is less than 3 hr using one computational node per frequency. Figure 9. View largeDownload slide Inversion results for the exploration model using short offsets (10 and 20 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) targets. The white box shows the inversion area. Figure 9. View largeDownload slide Inversion results for the exploration model using short offsets (10 and 20 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) targets. The white box shows the inversion area. Figure 10. View largeDownload slide Inversion results for the exploration model using medium offsets (15 and 30 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) targets. The white box shows the inversion area. Figure 10. View largeDownload slide Inversion results for the exploration model using medium offsets (15 and 30 m). XZ slice at the centre of the model (top) and two YZ slices through the first (bottom left) and the second (bottom right) targets. The white box shows the inversion area. DISCUSSION DDR tools provide multi-spacing and multi-frequency data sets and thus are no longer restricted to the detection of a single boundary with simplified inversion methods. These logging tools enable reliable estimation of spatial distributions of electrical resistivity several tens of metres away from the well. The interpretation of DDR data in complex multilayered structures should rely on full 3-D modelling and inversion. However, the high computational cost of 3-D numerical methods often constrains the development and usage of inversion-based interpretation. Our numerical results show that by making use of recent advances in linear algebra and parallel computing, the full 3-D inversion of DDR data can be performed in practical CPU times. The inversion of the fine-scale monitoring condition took about 30 hr of wall clock time using one computational node per frequency. The exploration test case was inverted by using smaller grids in 3 hr of CPU time. With appropriate computing resources, 3-D multisource inversion similar to our reservoir monitoring example can be performed within a few hours. A computationally efficient method based on unstructured grids would allow for further improvements in performance. Optimal transmitter–receiver configurations vary for different application conditions and frequencies. With the acquisition geometry being limited when compared to surface EM surveys, deeper zones are often shadowed by nearby resistivity inhomogeneities. Large source–receiver offsets (double of the target depth for the second receiver) provide advantageous inversion results compared to short offsets. The combination of resistivity measurements acquired with the DDR tools and seismic amplitude measurements will make possible the optimal steering of wells in geologically complex reservoirs. CONCLUSIONS We introduce and benchmark a new inversion method for 3-D DDR measurements in complex heterogeneous reservoirs penetrated by high-angle and horizontal wells. The key ingredients of our numerical scheme are the following. The minimization scheme employs the NLCG algorithm whose accuracy and performance is improved by invoking a pre-conditioning scheme. We use the Wolfe conditions to evaluate whether the step length is sufficient to provide near optimal convergence. Forward modelling of multicomponent triaxial induction data for multiple transmitter–receiver configurations uses a parallel direct solver. All components of the scheme, namely forward modelling, inversion algorithm, regularization, and pre-conditioner, are implemented in parallel and are well suited for execution on modern computer architectures. The workflow is validated on synthetic data sets that include monitoring and exploration borehole induction logging conditions. Numerical results verify that large transmitter–receiver offsets increase the depths of investigation to several tens of metres. Inversion of the monitoring example shows that the developed inversion-based method can determine the position of the oil-water contact in early water breakthrough zones, as well as the resistivity of the nearest formations. At the same time, the resistivity variations located at larger distances from the wellbore are not detected. These variations are most likely shadowed by water-rise zones. In the exploration conditions, the inversion algorithm provides accurate estimations of formation resistivity as well. The detection of deep (12–16 m away from the wellbore) targets requires offsets of 30 m to be reliably captured. These zones are often hidden by the resistivity variations in the vicinity of the well. These limitations of the acquisition geometry can be mitigated, for example, by using additional structural information available from other data sets, such as seismic amplitudes, to constrain the target formations. The synthetic data sets used in these simulations include 2 per cent Gaussian noise. Additional simulations with higher levels of noise verify the reliability of the inversion process. In particular, when some of the resistivity measurements are significantly noisier than the others (or missing at all). The forward modelling uses the finite difference method with a direct solver. For a multisource triaxial induction problem, this method enables a large reduction in the computational cost when compared to schemes based on iterative solvers. Future work will study the reliability and efficiency of the finite element method for this class of problems. Extending the inversion-based method to unstructured grids to reduce the number of unknowns in the resulting linear systems and thus to use solvers in a more efficient way is the subject of a forthcoming study. ACKNOWLEDGEMENTS This publication was made possible in part by the CSIRO Professorial Chair in Computational Geoscience at Curtin University and the Deep Earth Imaging Enterprise Future Science Platforms of the CSIRO. Additional support was provided by the National Priorities Research Program grant 7-1482-1-278 from the Qatar National Research Fund (a member of The Qatar Foundation) and by the European Union's Horizon 2020 Research and Innovation Program of the Marie Skłodowska-Curie grant agreement No. 644202. The J. Tinsley Oden Faculty Fellowship Research Program at the Institute for Computational Engineering and Sciences (ICES) of the University of Texas at Austin has partially supported the visits of VMC to ICES. Computational resources were provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia. The authors also acknowledge the University of Texas at Austin's Research Consortium on Formation Evaluation, jointly sponsored by Anadarko, Aramco, Baker-Hughes, BHP Billiton, BP, China Oilfield Services LTD., ConocoPhillips, Det Norske, ENI, ExxonMobil, Hess, Paradigm, Petrobras, Repsol, DEA, Schlumberger, Shell, Statoil, TOTAL, Wintershall and Woodside Petroleum Limited. The authors would like to thank Ute Weckmann, Colin Farquharson and one anonymous reviewer for their constructive suggestions, which helped to improve the manuscript. REFERENCES Abubakar A., Habashy T.M., Druskin V.L., Knizhnerman L., Alumbaugh D., 2008. 2.5 D forward and inverse modeling for interpreting low-frequency electromagnetic measurements, Geophysics , 73( 4), F165– F177. Google Scholar CrossRef Search ADS   Bakr S.A., Pardo D., Torres-Verdín C., 2017. Fast inversion of logging-while-drilling resistivity measurements acquired in multiple wells, Geophysics , 82( 3), E111– E120. Google Scholar CrossRef Search ADS   Bensdorp S., Petersen S.A., Olsen P.A., Antonsen F., van den Berg P.M., Fokkema J.T., 2016. An approximate 3D inversion method for inversion of single-well induction-logging responses, Geophysics , 81( 1), E43– E56. Google Scholar CrossRef Search ADS   Davydycheva S., Druskin V., Habashy T., 2003. An efficient finite-difference scheme for electromagnetic logging in 3D anisotropic inhomogeneous media, Geophysics , 68( 5), 1525– 1536. Google Scholar CrossRef Search ADS   dos Santos A.W., Pestana R.C., 2015. Time-domain multiscale full-waveform inversion using the rapid expansion method and efficient step-length estimation, Geophysics , 80( 4), R203– R216. Google Scholar CrossRef Search ADS   Dupuis C., Denichou J.M., 2015. Automatic inversion of deep-directional-resistivity measurements for well placement and reservoir description, Leading Edge , 34( 5), 504– 512. Google Scholar CrossRef Search ADS   Ezioba U., Denichou J.M., 2014. Mapping-while-drilling system improves well placement and field development, J. Pet. Technol. , 66( 08), 32– 35. Google Scholar CrossRef Search ADS   Farquharson C.G., Oldenburg D.W., 2004. A comparison of automatic techniques for estimating the regularization parameter in non-linear inverse problems, Geophys. J. Int. , 156( 3), 411– 425. Google Scholar CrossRef Search ADS   Hager W.W., Zhang H., 2006. A survey of nonlinear conjugate gradient methods, Pac. J. Optim. , 2( 1), 35– 58. Hou J., Mallan R.K., Torres-Verdín C., 2006. Finite-difference simulation of borehole EM measurements in 3D anisotropic media using coupled scalar-vector potentials, Geophysics , 71( 5), G225– G233. Google Scholar CrossRef Search ADS   Ijasan O., Torres-Verdín C., Preeg W.E., Rasmus J., Stockhausen E.J., 2014. Inversion-based interpretation of logging-while-drilling resistivity and nuclear measurements: field examples of application in high-angle and horizontal wells, Petrophysics , 55( 05), 374– 391. Jaysaval P., Shantsev D.V., de Ryhove S.D.L.K., Bratteland T., 2016. Fully anisotropic 3-D EM modelling on a Lebedev grid with a multigrid pre-conditioner, Geophys. J. Int. , 207( 3), 1554– 1572. Google Scholar CrossRef Search ADS   Kriegshauser B., Fanini O., Forgang S., Itskovich G., Rabinovich M., Tabarovsky L., Yu L., Epov M., 2000. A new multi-component induction logging tool to resolve anisotropic formations, in 40th Annual Logging Symposium , Society of Petrophysicists and Well Log Analysts. 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( 5), 499– 524. Google Scholar CrossRef Search ADS   Newman G.A., Alumbaugh D.L., 1997. Three-dimensional massively parallel electromagnetic inversion—I. Theory, Geophys. J. Int. , 128( 2), 345– 354. Google Scholar CrossRef Search ADS   Newman G.A., Alumbaugh D.L., 2002. Three-dimensional induction logging problems: Part 2—a finite-difference solution, Geophysics , 67( 2), 484– 491. Google Scholar CrossRef Search ADS   Newman G.A., Boggs P.T., 2004. Solution accelerators for large-scale three-dimensional electromagnetic inverse problems, Inverse Problems , 20( 6), S151– S170. Nocedal J., Wright S., 2006. Numerical Optimization , Springer Science & Business Media. Oldenburg D.W., Li Y., Farquharson C.G., Kowalczyk P., Aravanis T., King A., Zhang P., Watts A., 1998. Applications of geophysical inversions in mineral exploration, Leading Edge , 17( 4), 461– 465. Google Scholar CrossRef Search ADS   Pardo D., Demkowicz L., Torres-Verdín C., Paszynski M., 2006. Two-dimensional high-accuracy simulation of resistivity logging-while-drilling (LWD) measurements using a self-adaptive goal-oriented hp finite element method, SIAM J. Appl. Math. , 66( 6), 2085– 2106. Google Scholar CrossRef Search ADS   Pardo D., Matuszyk P.J., Torres‐Verdín C., Mora A., Muga I., Calo V.M., 2013. Influence of borehole-eccentred tools on wireline and logging‐while‐drilling sonic logging measurements, Geophys. Prospect. , 61( s1), 268– 283. Google Scholar CrossRef Search ADS   Plessix R.E., 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications, Geophys. J. Int. , 167( 2), 495– 503. Google Scholar CrossRef Search ADS   Polak E., Ribiere G., 1969. Note sur la convergence de méthodes de directions conjuguées, Rev. Fr. Inform. Rech. Oper., série rouge , 3( 1), 35– 43. Puzyrev V., Torres-Verdín C., 2016. Inversion-based interpretation of 3D borehole directional resistivity measurements acquired in high-angle and horizontal wells, in 2016 SEG Annual Meeting, SEG, Expanded Abstracts , pp. 724– 728. Puzyrev V., Koric S., Wilkin S., 2016. Evaluation of parallel direct sparse linear solvers in electromagnetic geophysical problems, Comput. Geosci. , 89, 79– 87. Google Scholar CrossRef Search ADS   Puzyrev V., Vilamajo E., Queralt P., Ledo J., Marcuello A., 2017. Three-dimensional modeling of the casing effect in onshore controlled-source electromagnetic surveys, Surv. Geophys. , 38( 2), 527– 545. Google Scholar CrossRef Search ADS   Rodi W., Mackie R.L., 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion, Geophysics , 66( 1), 174– 187. Google Scholar CrossRef Search ADS   Schenk O., Gartner K., 2004. Solving unsymmetric sparse systems of linear equations with PARDISO, Future Gener. Comput. Syst. , 20, 475– 487. Google Scholar CrossRef Search ADS   Seydoux J., Legendre E., Mirto E., Dupuis C., Denichou J.M., Bennett N., Kutiev G., Kuchenbecker M., Morriss C., Yang L., 2014. Full 3D deep directional resistivity measurements optimize well placement and provide reservoir-scale imaging while drilling, 55th Annual Logging Symposium , Society of Petrophysicists and Well-Log Analysts. Shantsev D.V., Jaysaval P., de Ryhove S.D.L.K., Amestoy P.R., Buttari A., L’Excellent J.Y., Mary T., 2017. Large-scale 3-D EM modelling with a Block Low-Rank multifrontal direct solver, Geophys. J. Int. , 209( 3), 1558– 1571. Google Scholar CrossRef Search ADS   Thiel M., Dupuis C., Omeragic D., 2016. Fast 2D inversion for enhanced real-time reservoir mapping while drilling, in SPE Annual Technical Conference and Exhibition , Society of Petroleum Engineers. Wang G.L., Torres-Verdín C., Gianzero S., 2009. Fast simulation of triaxial borehole induction measurements acquired in axially symmetrical and transversely isotropic media, Geophysics , 74( 6), E233– E249. Google Scholar CrossRef Search ADS   Ward S.H., Hohmann G.W., 1988. Electromagnetic theory for geophysical applications, in Electromagnetic Methods in Applied Geophysics , Vol. 1, pp. 131– 311, ed. Nabighian M.N., Society of Exploration Geophysicists. Wiese T., Greenhalgh S., Zhou B., Greenhalgh M., Marescot L., 2015. Resistivity inversion in 2-D anisotropic media: numerical experiments, Geophys. J. Int. , 201( 1), 247– 266. Google Scholar CrossRef Search ADS   Wolfe P., 1969. Convergence conditions for ascent methods, SIAM Rev. , 11( 2), 226– 235. Google Scholar CrossRef Search ADS   Zach J.J., Bjørke A.K., Støren T., Maaø F., 2008. 3D inversion of marine CSEM data using a fast finite-difference time-domain forward code and approximate Hessian-based optimization, in 2008 SEG Annual Meeting, SEG, Expanded Abstracts , pp. 614– 618. Zhdanov M.S., 2009. New advances in regularized inversion of gravity and electromagnetic data, Geophys. Prospect. , 57( 4), 463– 478. Google Scholar CrossRef Search ADS   © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

Journal

Geophysical Journal InternationalOxford University Press

Published: May 1, 2018

There are no references for this article.

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
discover and read the research
that matters to you.

Enjoy affordable access to
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.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off