Real-Time Prediction of Anisotropies kh and kv in Dual Axial Probe Formation Testing

Wilson C. Chin ^*

1. Massachusetts Institute of Technology, M.I.T., Cambridge, MA, USA

* Correspondence: Wilson C. Chin

Academic Editor: Grigorios L. Kyriakopoulos

Received: May 25, 2026 | Accepted: August 28, 2026 | Published: September 09, 2026

Journal of Energy and Power Technology 2026, Volume 8, Issue 3, doi:10.21926/jept.2603016

Recommended citation: Chin WC. Real-Time Prediction of Anisotropies kh and kv in Dual Axial Probe Formation Testing. Journal of Energy and Power Technology 2026; 8(3): 016; doi:10.21926/jept.2603016.

© 2026 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction and Background

Drilling and exploration in oil and gas are important to societal progress, economic and infrastructure development. In technical terms, an understanding of the reservoir and how petroleum fluids flow as a function of viscosity and rock structure is critical. Anisotropy, the most relevant flow property, was recognized early on in Muskat [1], the industry’s classic book on Darcy analysis. Until the mid-1990s, estimates were obtained from pressure drop and flow rate measurements associated with steady flow field tests. However, as newly discovered resources offered heightened resistance to flow, such tests required increasing time resources and well logging expense. In 1997, this author introduced a real-time method in Halliburton’s U.S. Patent 5,703,286 (see Proett, Chin and Chen [2]), supporting a commercially successful service known as GeoTap™. The patent provided isotropic permeability estimates rapidly, conveniently and accurately, allowing operators to log oil wells with higher point densities and greater efficiency. This breakthrough reduced dependence on resistivity, acoustic and nuclear instruments, which provided only indirect clues about the formation. It spearheaded two decades of progress supporting drilling efficiency, energy exploration and oilfield development. However, the 1997 method applied strictly to isotropic media, and extensions to general anisotropies have not been forthcoming. Governing equations for transversely isotropic media with flowline storage and surface skin effects were presented in Proett, Chin and Mandal [3] but a strategy for their solution was not available at the time.

Our objective in predicting both individual permeabilities kh and kv from real-time Formation Testing While Drilling (FTWD) pressure data is years in the making. In 2024, the author contacted Mark Proett, co-inventor of Halliburton’s GeoTap™ method, regarding industry status on inverse methods. The 1997 method, using a formal partial differential equation formulation for transient Darcy liquid flow, improved upon the company’s earlier U.S. Patent 5,602,334 due to Proett and Waid [4]. This work employed simple exponential formulas based on heuristic engineering arguments. The model was similar to Baker Hughes’ “Formation Rate Analysis (FRA)” method, based on U.S. Patent 5,708,204 due to Kasap [5]-this is further discussed in Kasap et al. [6], which also used similar elementary physical arguments. All three methods applied to isotropic formations only, with k = kh = kv. Validations for the k_eff inverse procedure were not available until 2008, when Halliburton’s sponsorship of a Doctoral Thesis at The University of Texas at Austin was completed, e.g., see Lee [7]. Lee’s laboratory based effort, requiring years and consuming significant expense, demonstrated that predicted k_eff values were consistent with experimental observation. Still, predictions for both kh and kv would not be available, despite the company’s dual probe capabilities, and only now, as reported in this paper, is a real-time “kh and kv” available that is supported by validations drawing on the synthetic pressure forward-inverse strategy developed in Section 2.

Low Reynolds number Darcy flows through porous media are not new to the petroleum industry. Aside from Muskat’s seminal contributions early on, other well-known works are available, e.g., Dake [8], Scheidegger [9], Streeter [10] and Yih [11]. Within this broad flow category, “well testing,” that is, formation characterization by analyzing pressure transients in cylindrical wells, is a long established discipline. Key references include Horne [12], Peaceman [13], Raghavan [14], Sabet [15], Stanislav and Kabir [16], and Streltsova [17]. While qualitative similarities between well testing and formation testing are seen, the two cannot be more different. For one, the former flows are cylindrical while the latter are spherical for isotropic media and ellipsoidal for anisotropic formations. This difference is not superficial-for instance, under constant rate pumping, cylindrical flows are forever transient, while spherical and ellipsoidal flow reach steady equilibrium in finite times. The formation testing literature is also well developed, e.g., Goode and Thambynayagam [18], Joseph and Koedetitz [19], Moran and Finklea [20], Schlumberger Staff [21], SPWLA Staff [22] and SPWLA Staff [23].

While various authors have tackled the forward problem in different simplified limits, none had addressed the general problem embodying flowline storage, skin effects and arbitrary flow rates. It was not until Proett, Chin and Chen [2] that the forward problem was solved in “exact,” closed, analytical form in terms of “complex complementary error functions.” At the time, “exact” referred to a spherically symmetric source flows in isotropic media. When this solution was asymptotically evaluated in the “intermediate time” limit, the simplified exponential solutions in Proett and Waid [4] and Kasap [5], both based on heuristic arguments, were recovered, thus lending confidence to the new differential equation method and the exponential solutions already in commercial use. But heuristic arguments, unlike differential equation methods, are very limited. They do not provide guidance, for example, on extensions to anisotropic problems, flows with skin effect, and so on. The formal partial differential equation foundation established in Proett, Chin and Chen [2] did because one could mathematically incorporate additional general physical effects into the high-level formulation and in principle solve entire classes of important forward and inverse problems.

In the 2003-2005 time frame, the United States Department of Energy, through its Small Business Innovation Research (SBIR) organization, would fund Stratamagnetic Software, LLC, founded by this author, in extending the state-of-the-art in formation testing solutions to broader classes of problems. Its two initial contracts were subsequently extended by additional private equity funding. Extensive results were made public through six books published by John Wiley & Sons and Elsevier Science, namely, Chin [24,25], Chin [26], Chin et al. [27], Lu et al. [28] and Lu et al. [29]. These studies were mathematically and numerically oriented. The earliest 2014 reference provides numerous calculated examples and details related to solution strategy and purpose. The basic approaches laid the foundation for Halliburton’s GeoTap™ and China Oilfield Services Limited’s (COSL) EFDT™ faster “rational polynomial” methods developed by this author. For completeness, we give the references used in developing asymptotic methods, Laplace transform inversions and partial differential equation solutions. These are Abramowitz and Stegun [30], Bateman [31], Carnahan, Luther and Wilkes [32], Churchill [33], IMSL Library Reference Manual [34], Pipes [35] and Press et al. [36].

In Section 2, we summarize our solution for dual kh and kv prediction, providing the steps used in formulating and validating “forward” and “inverse” models. Before turning to details, we indicate that, in the 1990s, as documented in Proett, Chin and Mandal [3], commercial finite element software packages were used to formulate forward and inverse solutions, as shown in Figure 1 to include storage and skin effects. While solutions provided qualitatively useful insights, grid dependencies were problematic: change the grid, the answer changes. Moreover, numerical and round-off errors produce “artificial viscosity” effects, so-called by John von Neumann, the 1940s computing pioneer, that are difficult to quantify in terms of magnitude and rheology type. Figure 2 shows typical field operations personnel deploying formation testing instruments, while Figure 3 displays well logs with different pressure drawdown and buildup curves as well as permeability logs useful to infill drilling, hydraulic fracturing, wellbore stability, reservoir simulation, and of course, cash flow analysis and overall economic assessment.

Figure 1 Industry formation testers. Credit: Proett, Chin and Mandal [3].

Figure 2 Single and dual probe testers lowered into the well. Courtesy: COSL.

Figure 3 Pressure versus time records and well logs. Courtesy: Courtesy: Halliburton Energy Services [3].

2. Problem Formulation, Solution and Software Design

We organize this section into three readable parts, isotropic theory, anisotropic extensions, and complex complementary error function analysis. Only “zero skin” and “zero supercharge” models are described in this paper. More general models are given in Chin et al. [27] and later publications. Readers interested in greater mathematical detail than offered here should request an unabridged report from the author.

2.1 Isotropic Theory

Let P(r,t) represent transient Darcy fluid pressure, r and t being radial and time coordinates. P₀ denotes the constant initial and farfield pressure of the quiescent reservoir; further, φ, k, μ and c refer to porosity, permeability, liquid viscosity and compressibility. R_w is the effective spherical well radius, that is, actual radius corrected by a multiplicative geometric factor. V is the flowline volume internal to the formation tester while C represents fluid compressibility within V. Note C and c may differ by ten-fold due to phase segregation effects. Otherwise, C and c are often equal. Q(t) denotes the total volume flowrate produced by the source probe due to sandface production plus flowline volume storage effects arising from compressibility. The Initial-BVP without skin and supercharge effects is completely specified by

\[ \partial^2P(r,t)/\partial r^2+2/r\partial P/\partial r=(\phi\mu c/k)\partial P/\partial t \tag{1} \]

\[ P(r,t=0)=P_0 \tag{2} \]

\[ P(r=\infty,t)=P_0 \tag{3} \]

\[ (4\pi R_W^2k/\mu)\partial P(R_w,t)/\partial r-VC\ \partial P/\partial t=Q(t) \tag{4} \]

The “(4πR_w²k/μ)∂P(R_w,t)/∂r” represents the contribution to Q(t) from flow through the sandface. It is the product of the Darcy fluid velocity k/μ ∂P/∂r and the spherical surface area 4πR_w², while “VC ∂P/∂t” is the amount due to fluid expansion and decompression in the flowline. We also assume that initial and far field pressures are identical constants (the author’s post-2014 work allows borehole “over-balance” and “underbalance,” effects not considered here).

2.1.1 Dimensionless Solution Strategy

Equations 1-4 contain many physical parameters, however, dimensionless variables can be introduced which simplify mathematical analysis. Details are offered in Proett, Chin and Chen [2] and Chin et al. [27]. For our purposes, we note that a simplified, isotropic italicized BVP results, namely

\[ \partial^2p(r,t)/\partial r^2+2/r\ \partial p/\partial r=\partial p/\partial t \tag{5} \]

\[ p(r,0)=0 \tag{6} \]

\[ p(\infty,t)=0 \tag{7} \]

\[ \partial p(r_w,t)/\partial r-\partial p/\partial t=F(t) \tag{8} \]

This is free of explicit parameters except for the single dimensionless radius r_w = 4πR_w³φc/(VC) appearing only in the argument of Equation 8. To solve this system, the equations are first Laplace transformed in time, leading to simpler ordinary differential equations that can be inverted exactly in closed analytical form as highlighted in Figure 2. To further simplify analysis, constant fluid withdrawal (drawdown) or injection (buildup) rates are assumed. In this limit, the transform at the source probe r_w is p(r_w,s) = -1/{s(s + s^1/2+ r_w^-1)}. The transform p_exact(r,s) at any passive observation probe location r > r_w is also available. These must be “inverted” to recover actual dimensional pressures with time dependencies. Numerical inversions such as Stehfast’s method are inaccurate and unable to provide the exact “synthetic data” needed to validate the overall “forward and inverse” approach discussed later. Fortunately, it is possible to derive analytical solutions that are exact, and these are duplicated from Chin et al. [27], as follows,

\[ \begin{equation} \begin{aligned} p_{exact}(r_{w},t) &=\{1/(\beta_1-\beta_2)\}\{{\beta_1}^{-1}-{\beta_1}^{-1}\exp({\beta_1}^2t)\operatorname{erfc}(\beta_1\sqrt{t}) \\&-{\beta_2}^{-1}+{\beta_2}^{-1}\exp({\beta_2}^2t)\operatorname{erfc}(\beta_2\sqrt{t})\} \end{aligned} \end{equation} \tag{9} \]

\[ \begin{equation} \begin{aligned} & p_{exact}(r,t)=\{r_w/(r(\beta_1-\beta_2))\}\\& \times\{\{-\exp\left({\beta_1}^2t-\beta_1(r-r_w)\right)\mathrm{erfc}\left(\beta_1\sqrt{t}+(r-r_w)/(2\sqrt{t})\right)+\mathrm{erfc}\left((r-r_w)/2\sqrt{t})\right)\}/\beta_1 \\&-\{-\exp\left({\beta_2}^2t-\beta_2(r-r_w)\right)\mathrm{erfc}\left(\beta_2\sqrt{t}+(r-r_w)/(2\sqrt{t})\right)+\mathrm{erfc}\left((r-r_w)/2\sqrt{t})\right)\}/\beta_2\} \end{aligned} \end{equation} \tag{10} \]

Here $\beta_{1}=+1/2-1/2\sqrt{1-4r_{w}^{-1}}$ and $\beta_{2}=+1/2+1/2\sqrt{1-4r_{w}^{-1}}$ are complex constants. The “erfc” refers to “complex complementary error functions” with complex arguments. While the pressures in these transient formulas are formally exact, significant numerical evaluation errors can arise that are addressed later. While we do use exact numerical forward solutions in Section 3, it is also useful to derive “early, intermediate, and late-time” formulas for permeability, compressibility and porosity inverse prediction from pressure data. These are illustrated below, but for now, we cautiously note that asymptotic expansions of dimensionless results like Equations 9 and 10 can introduce ambiguities relating to pressure transient behavior in dimensional time. Asymptotic time regimes are defined dimensionlessly, “How early is early” or “how late is late” depend on fluid and formation properties which are, of course, unknown. Section 3 definitively addresses these questions using parameters common to field practice. The care taken to derive exact, analytical solutions, to carefully evaluate “erfc” functions, to introduce our use of accurate “synthetic” pressure data to validate inverse predictive methods, will lend significant credibility to our approach to providing both kh and kv using early time data.

2.1.2 Early Time Series Solution

Our use of complex complementary error functions in pressure transient analysis may be unfamiliar. To show that the above reduces to known conventional results at small times, we introduce power series approximations for the exponential and complementary error functions, that is, e^x = 1 + x + 1/2x² + 1/6x³ and $\mathrm{erfc}(x)=1-2x/\sqrt{\pi}+2x^3/(3\sqrt{\pi})+\cdots$ in the exact solution. For example, this yields

\[ \begin{aligned}p(r_{w},t)_{early-time\ from\ exact}&=-t+4t^{3/2}/3\sqrt{\pi}\\&+(1-r_w)t^2/(2r_w)+8(r_w-2)t^{5/2}/(15r_w\sqrt{\pi})+\cdots\end{aligned} \tag{11} \]

Equation 11 provides a formal power series solution in t valid for small times. If we retain the first “p(r_w,t) = -t” term only and return to dimensional variables, we obtain

\[ P(R_w,t)_{early-time\ from\ exact}=P_{0}-Q_0t/(VC) \tag{12} \]

This “very-early-time” solution describes a pressure linear in time; the slope depends on Q and storage VC only and no transport properties. In practice, pressure responses at very early time are used to predict flowline compressibility C. Volume V is known from hardware specifications, varying from tool to tool, and manufacturer to manufacturer. Dynamical effects due to viscosity and permeability are captured by retaining additional terms, but these transport properties are only significant at later times. Large values of VC can mask permeability effects at early times-thus, an interpretation model allowing accurate formation evaluation is especially invaluable.

2.1.3 Large Time Asymptotic Solution

For large times, we can approximate the complementary error function using the asymptotic result $\mathrm{erfc}(x)=\exp(-x^2)\{1-1/(2x^2)+\cdots\}/\{x\sqrt{\pi}\}$. Then, our exact solution reduces to

\[ \begin{aligned}p(r_w,t)_{early-time\ from\ exact}=&-r_w+{r_w}^2/\sqrt{\pi}t+(2-r_w){r_w}^3/(2t^{3/2}\sqrt{\pi})\\&+3{r_w}^4(r_w-1)(r_w-3)/4t^{5/2}\sqrt{\pi}+\cdots\end{aligned} \tag{13} \]

In the late time limit, we retain the leading “$p(r_{w},t)=-r_{w}+{r_w}^{2}/\sqrt{(\pi t)}$” terms only and return to dimensional variables to obtain

\[ P(R_w,t)_{early-time\ from\ exact}=P_0-Q_0\mu/4\pi R_wk+\{Q_0\mu/(4\pi k)\}\sqrt{\{\phi\mu c/(\pi kt)\}}+\cdotp\cdotp\cdotp \tag{14} \]

Late time solutions are independent of flowline storage, depending only on transport properties like viscosity and permeability. They also predict an algebraic “inverse-square-root” time wise decline in pressure. Equation 14 can also be solved for the porosity φ, once other parameters in the formula are available. We have used algebraic manipulation software to calculate 100 additional terms in Equations 11-14 to understand their algebraic structure. Note that arbitrary flow rates Q(t) can also be modeled, as explained in Chin et al. [27], using convolution integral superpositions for Equations 9 and 10. Large dimensionless time solutions form the subject of this paper. Section 3 demonstrates how, for parameter ranges common to oil exploration, large dimensionless times based on asymptotics can be identical to early dimensional times based on subjective experience. This observation is crucial to dual permeability kh and kv prediction, and hence, to drilling and well logging costs as well as reservoir production strategies. So far, we have only addressed isotropic media. As noted, the FTWD literature is also restricted in this regard. However, our differential equation approach lends itself to powerful extensions, particularly to general transversely isotropic media considered next, and also to overbalanced and underbalanced drilling, multi-rate pumping, and liquid and gas applications considered in our post-2014 work.

2.2 Anisotropic Ellipsoidal Flow with Storage

The above work considers isotropic spherical flows that are relatively simple to solve. It is possible to derive exact solutions for transversely isotropic ellipsoidal flows with flowline storage in homogeneous media. This involves some complexity but, interestingly, final results can be expressed in terms of isotropic solutions. In general, an “effective permeability” is desired where P(x,y,z,t) satisfies k_xP_xx + k_yP_yy + k_zP_zz = φμcP_t for slightly compressible liquids. Often, scale changes like $x=\{(k_xk_yk_z)^{1/6}/\sqrt{k_x}\}x$, $y=\{(k_xk_yk_z)^{1/6}/\sqrt{k_y}\}y$ and $z=\{(k_xk_yk_z)^{1/6}/\sqrt{k_z}\}z$, are used. The assumption “P(x,y,z,t) = P(x,y,z,t)” leads to (k_xk_yk_z)^1/3(P_xx + P_yy + P_zz = φμcP_t), showing that (k_xk_yk_z)^1/3 is the effective permeability.

In “transversely isotropic” media, this reduces to (k_vk_h²)^1/3 since k_v = k_z and k_h = k_x = k_y. Hence, the claim is often made that isotropic predictions apply to anisotropic application with “k” replaced by k_h^2/3k_v^1/3. This is true, but only fortuitously, since scale arguments applied to partial differential equations alone provide only a partial proof. It is necessary to additionally rescale flow rate conditions, and then, demonstrate how a similar dimensionless formulation exists that reduces to Equations 5-8. This proof, offered in Chin et al. [27], is now summarized. Earlier we studied isotropic flow using the spherical equation ∂²P/∂r² + 2/r ∂P/∂r = (φμc/k)∂P/∂t. This was solved together with (4πR_w²k/μ)∂P(R_w,t)/∂r - VC ∂P/∂t = Q(t). In both of these equations, “k” represented an isotropic permeability constant in all directions. The volume flux (4πR_w²k/μ) ∂P/∂r, by virtue of spherical symmetry, was simply the product between the Darcy velocity (k/μ) ∂P/∂r and the spherical surface area 4∂R_w²associated with an effective radius R_w.

For anisotropic flow, complications related to differential equation and flow rate condition arise. For the former, consider k_v∂²P/∂z2 + k_h(∂²P/∂x2 + ∂²P/∂y²) = φμc ∂P/∂t. We first re-scale x, y and z so that a term proportional to “∂²P/∂x2 + ∂²P/∂y² + ∂²P/∂z2” appears. This is converted to “∂²P/∂r² + 2/r ∂P/∂r” using spherical curvilinear coordinate transforms. This shows how spherical volumes become ellipsoidal, with two k_h axes and one k_y. The second problem, relating to “(4πR_w²k/μ) ∂P/∂r,” requires more than multiplying the velocity k/μ) ∂P/∂r by the area 4πR_w². The “scalar, dot product” between (i) the Darcy velocity q (expressed in terms of spherical ∂P components), and (ii) the unit normal n to the ellipsoidal surface “z = f(x,y)” must be taken at all ellipsoidal surface points and integrated over the enclosed volume. The algebra is tedious, but a dimensionless p can be found that satisfies Equations 5-8 although with different normalizations. Hence algorithms solving isotropic flows will solve anisotropic ones with minor modification. The analytical, closed form, exact nature of our solutions offers significant advantages. That is, accuracy can be improved incrementally if limitations identified at each stage of any required evaluations can be removed. This process is explained in our presentation. A sketch of a transversely isotropic flow domain appears in Figure 4.

Figure 4 Transversely isotropic ellipsoidal domain (from Chin et al. [27]).

2.3 Permeability Strategy and Algorithm Requirements

Our objective is real-time dual kh and kv permeability prediction in FTWD applications. We will demonstrate our inverse methods are correct, for ranges of permeability encountered in oil exploration, using derived numerical methods that are highly accurate. How is this achieved? The strategy requires a “Forward Simulator,” which creates exact pressure versus time results, or “synthetic data,” when k_h, k_v, φ, μ, c, V, C, Q, R_w and P₀ inputs are specified. Two “Inverse Simulators,” based on dimensional re-expressions of Equations 9 and 10, must reproduce k_h, k_v and P₀ to sufficient accuracy without knowledge of φ, c, V and C. This cannot be achieved with conventional numerical methods, e.g., finite difference, finite element, or other methods, because mesh-dependency leads to truncation errors masquerading as “artificial viscosities.” This is true of commercial CFD models and second-order schemes developed by the author.

In fact, each component strategy requires rigorously designed procedures. For example, we solved Equations 1-4 for time-varying pressures at source and observation points in terms of exact “complex complementary error functions” erect(z). We derived short and late time asymptotic pressure versus time solutions for compressibility, permeability and porosity (this paper focuses on permeability analysis, with details offered in Section 3). When any three (p,t) points are reasonably chosen, our objective is accurate k_h^2/3k_v^1/3 prediction for single, source probe only tools. If dual probe pressures are available at a specific “late” time, we seek accurate values for both k_h and k_v. We know whether or not our predictions are correct because these were chosen to create the synthetic data used to create exact transient pressure histories. Moreover, our method applies to all dip angles. In our Section 3 examples, we assumed a 45 deg dip so that our forward and inverse algorithms would execute all internal software logic loops-successful predictions, in this sense, imply reliable, “bug-free” program code. During development, we found that constant rate withdrawals did not always imply physically required monotonic time pressure drawdowns at observation probes. Slight rises were observed on occasion while oscillations were never found. This behavior was traced to errors in scientific libraries accessed by language compilers used to develop code. Errors arose from unpredictable and indeterminate products of “very small numbers, M” and “very large numbers, N”. To remove this problem, commonly used erfc routines were replaced by improved algorithms focusing on finite “M × N” entities.

2.3.1 Complex Complementary Error Function

All of our models were formally developed in terms of erfc(z), that is, the complex complementary error function with complex arguments. For example, the dimensionless source radius r_w = 4πR_w³ φc/(VC) arises in isotropic flow. This appears in the complex $\beta_1=+1/2-1/2\sqrt{(1-4r_w^{-1})}$ and $\beta_2=+1/2+1/2\sqrt{(1-4r_w^{-1})}$ constants. These parameters are found in dimensionless source and observation probe pressures. Typically, the β’s occur in terms like $\mathrm{erfc}(\beta\sqrt{t})$ and $\mathrm{erfc}(\beta\sqrt t+(r-r_w)/(2\sqrt t))$. Although our formal solutions are exact, for many parameters of practical formation testing interest, common computer evaluations using standard software libraries of erfc(z) can be time-consuming or divergent, particularly for observation probe calls. This is the case for certain flowline volumes and permeabilities encountered in field application. This possibility is not acceptable in oil field interpretation, where thousands of accurate evaluations may be required at a single depth in real-time, and economic consequences can be significant. Detailed analysis appears in Chin et al. [27]. We simply quote results, referring to Figure 5 and Figure 6.

Figure 5 Standard evaluation of erfc(z) for x > 0, y > 0 (x and y are not space variables).

Figure 6 Improved numerical approach.

Figure 5 displays Quadrant 1 ranges of z where x > 0, y > 0, not to be confused with spatial coordinates, for which erfc(x + iy) can and cannot be computed by standard algorithms. The lower left (blue) zone shows where solutions cannot be computed due to overflow, while elevated right (green) highlights where solutions can be found. For solutions to be useful, robust numerics are required. Our computations focus on the product “exp(-z²) erfc(-iz)” using efficient series and continued fractions. The desired erfc(z) is evaluated by back-calculation. This strategy allows a much broader range of arguments for successful calculation, and as shown in Figure 6, enables highly efficient regression analysis for permeability and anisotropy prediction. The far left (blue) zone still represents solutions that cannot be computed; fortunately, the corresponding arguments do not represent parameters often used in formation testing. The far right (green) domain indicates solutions possible using conventional approaches, while the new and sizable middle (red) zone depicts an increased range of newly available computations for erfc(z). Accuracy and speed are emphasized in the method, e.g., hundreds of forward evaluations per second are possible on typical personal computers.

2.3.2 Numerical Methods, Analytical Solutions and Errors

Our formation testing objectives focused on kh and kv prediction in anisotropic, homogeneous media using single and dual axial probe data. Applications include kh >> kv important to common sedimentary layers, kh << kv in vertically fractured media, at general dip angles, and isotropic kh = kv problems. In particular, two methods under continual refinement over the past three decades are Inverse Simulator #1 and Inverse Simulator #2 explained and validated in Section 3. Oilfield publications often cite modeling successes by reporting concurrent petroleum discoveries. The logic is flawed True validations require agreement with exact datasets whose accuracies are above reproach, with analytics explained and source code changes available for examination. In the late 1990s, Proett and Waid [4] and Kasap [5] provided reasonable heuristic arguments leading to acceptable results. Lee [7] gave lab results showing that our 1997 methods were accurate. However, rigorous validations were impossible as quality control standards (e.g., synthetic datasets) were difficult to define. Without a differential equation approach, extensions to anisotropy, supercharge, underbalanced drilling and multi-rate superpositions would not be possible. The author’s development philosophy was also shaped by his early experiences.

For example, in developing signal processing methods for reflection removal in Measurement While Drilling (MWD) telemetry, good datasets offering high accuracies were required against which filtering schemes were evaluated. Numerically created standards were inaccurate-an input speed of 5,000 ft/sec might create 4,500-5,500 ft/sec signals depending on space and time step sizes selected. In MWD siren design, commercial simulators offered credible solutions through impressive color graphics. However, solutions changed as grid sizes and aspect ratios changed, a consequence of what John von Neumann, the 1940s computing pioneer, termed “artificial viscosity.” In short, CFD methods introduce round-off and truncation errors that contaminate solutions, not only with variable viscosities, but unpredictable nonlinear rheologies. In turbine design, the same commercial simulators would predict inaccurate no-load torques and power levels not useful for design. Conversations with algorithm designers revealed an unawareness of “Kutta’s condition,” a principle well known in aerodynamics.

With respect to the exact, closed form analytical solution in Proett, Chin and Chen [2], there is no question that, if solutions were evaluated accurately, the work would have been useful in job planning and timely inverse applications. However, the published isotropic solution was not studied until the author’s formation testing research was funded by the United States Department of Energy’s Small Business Innovation Research (SBIR) unit in 2003-2005. This author had assumed company programmers would deploy the new results in its software. That this was not so implied additional delays that led to problems uncovered and solved in Figure 5 and Figure 6. In the meantime, the author had shown that the transversely isotropic flows considered in this paper, satisfy an identical dimensionless boundary value problem for a different set of normalizations, an observation that provided greater impetus to validating error function properties and inverse solutions once and for all.

2.3.3 Solutions and Error Analysis

There are no errors in the exact, closed form analytical solutions, but limitations exist in available scientific program libraries used to evaluate complex error functions. Most limitations have been removed by us although, as Figure 6 shows, some remain. Aside from these exceptions, our inverse validation strategy was straightforward. Enter k_h, k_v and other parameters into the forward simulator to create synthetic pressure versus time data, and then, enter arbitrarily selected computed pressure and time pairs into our inverse simulators to predict effective k_eff, and k_h and k_v. If agreement is found, the method works. If not, it doesn’t. This test is more demanding than it appears. For example, the forward model requires c, C, V, P₀ and φ inputs, which our inverse simulators should not and must not at intermediate times. If agreement is found, and this is often, this provides excellent evidence that our method not only solves the right equations, but that these embody the correct physics.

Finally, several useful reminders. The forward simulator also assists in oilfield well logging work. Operationally, it provides means for “job planning” activities. Geologists and drillers often know, based on prior experience, approximate fluid and rock properties. Forward simulations predict pressures as they also depend on flow rate, dip angle and hardware parameters. If insufficient, simply calculate the increased flow rate or decreased nozzle size needed. If inverse methods will be used for permeability prediction, strong pressures are required at both source and observation probes. How is this requirement guaranteed? The answer is simple. Again, vary rate or nozzle size. A second reminder relates flowline storage effects to pressure waveform. In field work and design, VC is deliberately small in value, since “large values distort the pressure.” However, it can be shown that despite significant changes to shape, permeability predictions are unaffected. An example calculation is studied in Section 3, in which a 300 cc flowline is increased to 30,000 cc and permeabilities remain unchanged. This is obviously important to truly evaluating permeability. We have used this property recently in hardware design. In one application, pressure equilibrium was achieved so rapidly that the small time measurement window was insufficient for accurate data collection. The problem was solved by increasing the VC product, guided by the software in this paper. Finally, users should ensure that exact analytical formulations such as ours are used, noting that the guidance offered in Figure 5 and Figure 6 is crucial to quality control and not simply academic exercises.

3. Discussion-Simulations and Validations for Sedimentary Layers and Vertical Fractures Over Wide Flowline Volume Ranges for Data Quality Control

Our “transversely isotropic” model for two horizontal principal axes having identical kh permeabilities, and a third with kv, does not impose restrictions on assumed flow symmetries. We support tools oriented at any dip angle, whether vertical, deviated or horizontal. We will describe the analysis and design methods for forward and two inverse simulators before proceeding with suites of rigorous validations. We cover simulations for sedimentary layers with kh >> kv, vertical fracture applications with kh << kv, at a dip angle of 45 degrees so that all source code logic loops are evaluated, and finally, simple isotropic problems with kh = kv. Henceforth, subscript “h” and “v” notations will be omitted for clarity in calculated results. In addition, we explore dynamical properties associated with the flowline storage product VC. Low to high VC values change the shape of the pressure waveform within the measurement window. However, these differences do not affect the value of the predicted permeabilities k_eff, kh or kv. This can be proven mathematically, but is not well known operationally. If steady state pressures are achieved too rapidly for accurate measurement, the problem is solved by extending formation tester internal tubing length or storage volumes. If time scales are excessively long that wait times at the drilling rig are lengthy and costly, the flowline should be shortened. VC flow properties are easily evaluated with our software with the advantages implied by Figure 6. We will now explain our three simulators using Figures 7-10. For convenience, the mathematical symbols used on our derivations and examples are summarized in Table 1 below.

Figure 7 Exact forward pressure versus time simulator for job planning and synthetic data creation for inverse algorithm evaluation.

Figure 8 Dual axial probe formation tester tool at δ dip angle, conventions for Forward Simulator, and Inverse Simulators #1 and #2.

Figure 9 Inverse Simulator #1 for real-time kh and kv prediction for dual axial probe operation.

Figure 10 Inverse Simulator #2 for real-time k_eff prediction for single source probe operation.

Table 1 Formation Testing Input Parameters.

Forward Simulator-Creates synthetic source and observation probe data, when formation kh and kv are given, for use in evaluating inverse kh and kv predictions results based on asymptotic methods. While this software uses the exact, closed form analytical pressure transient algorithm derived in above Section 2, Figure 6 shows that the standard complex complementary error function evaluation in many scientific libraries is not entirely accurate. Physically, observation probe pressures due to constant rate fluid withdrawal should decrease monotonically with time. In the great majority of cases, our inverse methods have proven successful. In some instances, minor anomalies are apparent from plotted curves. Not too often, an expected drawdown curve may increase somewhat before continuing its decline with time. These are noted in calculated examples below. For such cases, improvements to complex erfc accuracy are needed and under active investigation. At this time, Figure 6 illustrates improvements achieved over Figure 5. The forward simulator menu is shown in Figure 7 and menu conventions are given in Figure 8.

Inverse Simulator #1-Predict individual kh and kv values from very late dimensionless time transient pressure drops ∆P_S and ∆P_d at source and observation probes. The exact analytical solutions for forward simulator above are used, that is, Equations 9 and 10 for source probe and observation probe positions, respectively. The two dimensionless equations are rewritten in dimensional physical variables. They are not interpreted as pressure solutions, but now viewed as implicit relations for kh and kv in terms of pressure drop variables. Within this framework, analysis shows that kh, kv and the anisotropy ratio kh/kv satisfy cubic equations of the form ( ) kh³ + ( ) kh + ( ) = 0, ( ) kv^3/2 + ( ) kv + ( ) = 0 and ( ) (kh/kv) + ( ) (kh/kv)^1/3 + ( ) = 0 where ( ) contain coefficients dependent on Q, μ, δ, R_w, ∆P_s, ∆P_d and L where a section of the formation tester appears in Figure 8. It is important to note, from the menu in Figure 9, that c, C, V and φ inputs are not needed in intermediate time expansions although these are entered the complete forward solution. This is consistent with overall physical assumptions since the flowline storage product VC is important only at very early time, while φc, arising at large times asymptotically, disappears at steady-state. Note that root (kh, kv) pairs may be entirely real or may contain complex conjugates. Only positive real numbers are physical significant, while complex roots may be acceptable if imaginary parts are minor, arising from inaccurate tool calibration or measured pressures.

Inverse Simulator #2-Predict “spherical” or “effective k_eff = kh^2/3kv^1/3” from intermediate time source probe data. As in Inverse Simulator #1 above, c, C, V and φ are not required inputs, the only required parameters being three sets of (p,t) source probe pressure drops, Q and R_w. The menu in Figure 10 allows “supercharge pressures” found in overbalanced drilling, for the authors post-2014 extended models. As this is not the case here, we set P_over = 0. Predictions are given below the yellow diagram, for pore pressure, effective mobility = k_eff/μ and compressibility. The latter is provided for reference only, as the time data used in the examples of this section are not compatible with very early-time physics. As derived earlier in Equation 12, the formula P(R_w,t)_{early-time from exact} = P₀ - Q₀ t/(VC) provides more accurate C results. Note that the “pore,” or hydrostatic pressures for forward analysis and predictions below are not dynamically significant. In our validations, they are chosen large enough so that drawdown pressures are always positive to allow clearer line plots. In real logging situations, the forward simulator and its synthetic pressures would not be used, and calculated pore pressures would be geologically important. The inverse solution strategy differs from Inverse Simulator #1. The intermediate time solution, again derived from the exact erfc solution, now takes on a different math structure, with coefficients A and B defining a nonlinear transcendental equation which is solved iteratively. In this case, A = μQ₀/(4πR_wk_h^2/3k_v^1/3) has dimensions of pressure while B = 4πR_wk_h^2/3k_v^1/3/(μVC) has units of inverse time. Detailed analysis shows how the kh^2/3kv^1/3, known as “effective permeability,” arises naturally, while for prior Inverse Simulator #1, kh and kv are individually solved. While some industry inverse methods require as many as ten pairs of (p,t) data plus additional smoothing, our method is quite robust as three (p,t) pairs generally suffice.

3.1 Validation Problems (Arranged, from High to Low Permeabilities)

This listing summarizes the sequence of validation problems considered from high to low permeabilities. Again, kh >> kv, kh < kv and kh = kv are considered, mostly with V = 300 cc and once a highly exaggerated V = 30,000 cc to demonstrate an important physical property. A 45 deg dip angle ensures that all forward and inverse software loops are executed and perform correctly. All of these cases are required to cover the parameter space relevant to our methods.

Example 1. High kh = 1,000 md, kv = 100 md (sedimentary layer).

Case 1. Formation too permeable, equilibration too rapid, insufficient time to take measurements. This field problem is common.

Case 2. Flowline volume increased 100× for demo purposes. This does not affect permeability prediction, an important physical fact not well known but can be proven mathematically.

Example 2. Two studies, kh >> kv and kh << kv.

Case 1. Assumed, kh = 100 md, kv = 10 md (sedimentary layer). Permeabilities can be predicted in first 10 sec of logging.

Case 2. Test example, kh = 100 md, kv = 1,000 md (vertical fractures).

Example 3. Assumptions, kh = 10 md, kv = 1 md (sedimentary layer).

Case 1. Permeabilities from t = 1,000 sec data.

Case 2. Permeabilities from t = 200 sec data.

Case 3. Permeabilities from t = 100 sec data.

Example 4. Assumptions, kh = 1 md, kv = 1 md (isotropic media).

Case 1. Permeability from t = 1,000 sec data.

Case 2. Permeability from t = 100 sec data.

Example 5. Assumptions, kh = 1 md, kv = 10 md (fractured media).

Case 1. Permeabilities from t = 1,000 sec data.

Case 2. Permeabilities from t = 200 sec data.

Example 6. Assumptions, kh = 0.1 md and kv = 0.01 md (sedimentary layer, very low permeability).

Example 7. Field and experimental applications.

Case 1. Review of effective permeability k_eff = kh^2/3kv^1/3 from early time source probe data.

Case 2. Anisotropy kh and kv from Inverse Simulator #2 (see above) plus field or lab measurements.

Case 3. Anisotropy (kh,kv) from Inverse Simulator #1 model for “dimensionless late time” conditions.

Example 8. Additional very low permeability isotropic validation.

3.2 Example 1: High kh = 1,000 md, kv = 100 md (Sedimentary Layer)

3.2.1 Case 1

Formation too permeable, equilibration too rapid, insufficient time to take measurements. This field problem is very common. Refer to Figure 11 for the menu inputs used.

Forward Solution.

Figure 11 Assumptions, source and observation probe drawdown pressures versus time.

Equilibrium is attained so rapidly that there is insufficient time to perform accurate, practical measurements. Also, pressure drops are too small. Inverse calculations cannot be performed with the assumed data. Let us next increase the flowline volume of 300 cc by 100 times, to 30,000 cc, which is highly exaggerated to demonstrate an important physical and mathematical property discussed earlier.

3.2.2 Case 2

Flowline volume increased 100× for demo purposes. This does not affect permeability prediction, an important physical fact not well known but can be proven mathematically and demonstrated below. Refer to Figure 12 for the menu inputs used.

Forward Solution.

Figure 12 Assumptions, dual probe drawdown pressures versus time (the slight upward trend in observation probe pressure versus time follows from erfc errors noted earlier).

These “forward solution” results for 30,000 cc are distorted relative to those in Case 1 but turn out to be usable for inverse calculations. However, the assumed kh = 1,000 md and kv = 100 md is very accurately predicted. In petroleum engineering, engineers avoid using long flowlines in tools because these “will distort the ‘true’ pressures and therefore permeabilities.” However, it can be shown that predicted permeabilities are independent of the flowline storage product VC. In highly mobile applications where pressures equilibrate too rapidly for accurate measurements, an excellent, proven strategy simply lengthens the flowline as necessary.

Using the final row of synthetic data above, we only need the two underlined pressure drops at the source and passive observation probes for the “inverse model” below in Figure 13. The inverse model will predict both kh and kv individually, without any knowledge of the compressibilities c and C, the flowline volume V, and the formation porosity φ. This is the main advantage of our inverse model, that is, the ability to distill information about anisotropic permeability without needing additional unimportant parameters.

Inverse Model.

Figure 13 Inverse model assumptions.

It is important to note that the correct anisotropic permeabilities are obtained even though flowline distortion, given the exaggerated V = 30,000 cc assumed, is significant. Also, multiple mathematical solutions may exist. Obviously, negative permeabilities and complex permeabilities with imaginary parts are disallowed, although minor discrepancies are permissible that may arise from experimental error. At other times, both kh >> kv and kh << kv are found. Ruling out either requires additional well logging data or trained geologist opinion.

3.3 Example 2: Two Studies, kh >> kv and kh << kv

3.3.1 Case 1

Assumed, kh = 100 md, kv = 10 md (sedimentary layer). Permeabilities can be predicted in first 10 sec of logging. Refer to Figure 14 and Figure 15 below for the menu inputs used.

Forward Solution.

Figure 14 Assumptions, dual probe drawdown pressures versus time.

Inverse Model.

Figure 15 Inverse model assumptions.

3.3.2 Case 2

Test example, kh = 100 md, kv = 1,000 md (vertical fractures). Refer to Figure 16 and Figure 17 for the menu inputs used.

Forward Solution.

Figure 16 Assumptions, dual probe drawdown pressures versus time (note unlikely observation probe pressure increases at right).

Inverse Model.

Figure 17 Inverse model assumptions.

3.4 Example 3: Assumptions, kh = 10 md, kv = 1 md (Sedimentary Layer)

3.4.1 Case 1

Permeabilities from t = 1,000 sec data. Refer to Figure 18 and Figure 19 for the menu inputs used.

Forward Solution.

Figure 18 Assumptions, dual probe drawdown pressures versus time.

Inverse Model (This Model Uses Underlined Data).

Figure 19 Inverse model assumptions.

3.4.2 Case 2

Permeabilities from t = 200 sec data. Refer to Figure 20 and Figure 21 for the menu inputs used.

Forward Solution.

Figure 20 Assumptions, dual probe drawdown pressures versus time.

Inverse Model (This Model Uses Underlined Data).

Figure 21 Inverse model assumptions.

3.4.3 Case 3

Permeabilities from t = 100 sec data. Refer to Figure 22 and Figure 23 for the menu inputs used.

Forward Solution.

Figure 22 Assumptions, dual probe drawdown pressures versus time.

Inverse Model (This kh and kv Model Uses Underlined Data).

Figure 23 Inverse model assumptions.

The above calculation for individual kh and kv uses pressure drops at both source and observation probes, and importantly, does not rely on compressibilities c and C, flowline volume V or formation porosity φ. When only highly transient data from a single source probe tool is available, this paper explains that the “spherical” or “effective permeability” k_eff = kh^2/3kv^1/3 can be predicted. This is untrue of conventional methods which, as noted, are restricted to isotropic problems. In this Example 3, the forward solution predicts synthetic pressures used to evaluate inverse permeability methods. For the kh = 10 md and kv = 1 md, we have exactly k_eff = 10^2/31^1/3 md or 4.642 md. The above inverse method gives an approximate k_eff = 10.7^2/30.877^1/3 = 4.648 md for a 0.1% error. The second “drawdown or buildup only” Inverse Simulator #2 (Figure 24) will not give kh and kv individually, but only k_eff = kh^2/3kv^1/3 = 4.713 md (bottom right, beneath yellow diagram) as shown. The compressibility given is not accurate as compressibility is dominant only at very early times – we had noted this and encouraged use of the more accurate Equation 12. The error in this case is 1.5%. Data used in Figure 24 are obtained at 0, 20 and 98 sec and shown in bolded lettering in the forward solution tabulation following Figure 22. While it is satisfying that all three effective permeabilities agree closely, this average is physically meaningful only when kh and kv are comparable in magnitude. In cases where kh and kv differ substantially, it is preferable for engineers in hydraulic fracturing, infill drilling, wellbore stability and reservoir engineering to work with individual kh and kv values from Inverse Simulator #1.

Inverse FTWD Model (Model Uses 0, 20 and 98 Sec) Underlined Data).

Figure 24 Inverse Simulator #2 model assumptions.

3.5 Example 4: Assumptions, kh = 1, kv = 1 (Isotropic Media)

3.5.1 Case 1

Permeability from t = 1,000 sec data. Refer to Figure 25 and Figure 26 for the menu inputs used.

Forward Solution.

Figure 25 Assumptions, dual probe drawdown pressures versus time.

Inverse Model (This Uses Underlined Pressure Data).

Figure 26 Inverse model assumptions.

3.5.2 Case 2

Permeability from t = 100 sec data. Refer to Figure 27 and Figure 28 for the menu inputs used.

Forward Solution.

Figure 27 Assumptions, dual probe drawdown pressures versus time.

Inverse Model (This Uses Underlined Pressure Data).

Figure 28 Inverse model assumptions.

3.6 Example 5: Assumptions, kh = 1 md, kv = 10 md (Fractured Media)

3.6.1 Case 1

Permeabilities from t = 1,000 sec data. Refer to Figure 29 and Figure 30 for the menu inputs used.

Forward Solution.

Figure 29 Assumptions, dual probe drawdown pressures versus time.

Inverse Model (This Uses Underlined Pressure Data).

Figure 30 Inverse model assumptions.

3.6.2 Case 2

Permeabilities from t = 200 sec data. Refer to Figure 31 and Figure 32 for the menu inputs used.

Forward Solution.

Figure 31 Assumptions, dual probe drawdown pressures versus time.

Inverse Model (This Uses Underlined Pressure Data).

Figure 32 Inverse model assumptions.

3.7 Example 6: Assumptions, kh = 0.1 md and kv = 0.01 md (Sedimentary Layer, Very Low Permeability)

Please refer to Figure 33 and Figure 34 for the menu inputs used in this example.

3.7.1 Forward Solution

Figure 33 Assumptions, dual probe drawdown pressures versus time.

3.7.2 Inverse Model (This Used Underlined Pressure Data)

Please refer to Figure 34 for the menu inputs used in this example.

Figure 34 Inverse model assumptions.

3.8 Example 7: Field and Experimental Applications

3.8.1 Case 1: Review of Effective Permeability k_eff = kh^2/3kv^1/3 from Early Time Source Probe Data

Again, to evaluate kh-kv inverse prediction accuracy, we require exact transient synthetic pressure data. This is obtained from the Forward Simulator, which provides accurate pressure versus time data from an exact, closed form, analytical solution. This solution is expressed in terms of erfc functions. Figure 35 summarizes all input parameters, emphasizing the anisotropic assumption kh = 10 md, kv = 1 md and a flow rate of 1 cc/s. Exact pressure responses for the first 10 seconds at both probes are plotted, with results for 110 sec later in our discussions.

Figure 35 Forward model assumptions, with source (top) and observation probe data (bottom), assuming zero (vertical well) dip angle.

At 0, 5 and 9.8 sec data, dual probe pressures without hydrostatic effects are shown above. This is pressure drawdown or fluid withdrawal, associated with negative pressure. Three (p,t)’s are used for early time inverse predictions based on the first 10 sec of data. These are (0.0, 0.0), (5.0, -248.70) and (9.8, -249.47). imes are widely separated, focused on obtaining permeability. Compressibility depends on very early time physics; a more accurate prediction uses Equation 12. Now we use these three drawdown points in Inverse Simulator #2 to predict k_eff in Figure 36. We set Pover = 0 as there is no supercharge pressure in this problem. The exact value from Forward Simulator inputs is k_eff = 10^2/3 1^1/3 = 4.6416.

Figure 36 Inverse Simulator #2 predicted effective mobility, middle block at bottom right, noting dip angle is not required for source probe only problems.

The predicted k_eff = 4.7353 implies 4.7353/4.6416 = 1.0202 gives a small 2% error. The pore pressure is correct at 0, since it is the actual “20,000-20,000 (hydrostatic)” = 0. Again, k_eff follows from intermediate time data and only requires source point pressures from the first few seconds. To obtain kh and kv individually, we can turn to Inverse Simulator #1 as illustrated in Examples 1-6.

3.8.2 Case 2: Anisotropy kh and kv from Inverse Simulator #2 (Case 1) Plus Field or Lab Measurements

From theory, early time source probe only data gives k_eff = kh^2/3kv^1/3, an average meaningful physically only if kh and kv are close. In many applications, e.g., hydraulic fracturing, infill drilling, reservoir simulation and wellbore stability, kh and kv are needed individually. This was demonstrated earlier using Inverse Simulator #1 which uses final pressure drop values. Another physically based method for kh and kv is possible at relatively early times. We now view “kh^2/3kv^1/3 = k#” as a constraint with a known value relating kh to kv, which we can re-write as kh = (k#)^3/2/(kv)^1/2. For the above calculation, this simplifies to k# = 4.7353. This leads to the relationship kh = (k#)^3/2/(kv)^1/2 = (4.7353)^3/2/kv^1/2 = 10.30/kv^1/2, allowing us to calculate a wider range of allowed (kh, kv) pairs, e.g.,

We can now construct a more detailed (kh, kv) table to guide the plotting exercise that follows,

We next create pressure versus time plots at the observation probe for each (kh,kv) pair using the Forward Simulator. These curves are very different, as shown in Figure 37. Only six illustrative pressure drawdown examples are shown due to space limitations (similar ideas apply to pressure buildups with opposite flow rate signs).

Figure 37 Possible “pressure versus time” curves at observation probe.

Which curves correspond to an actual (kh, kv) solution? In this example, they are determined from measured observation probe data as illustrated in Figure 38. When experimental results like the “circle” and “square” points coincide with calculated curves, the sought (kh, kv) pair is circled and listed at the top of the figure. A third experiment shows “triangle” points. This scatter is more likely in practice, and the “best fit” line then defines the correct (kh, kv). This procedure requires history matching over a range of times, rather than at just one point, preferable since empirical data is rarely accurate enough that a single data point suffices. We do not recommend automated least squares and curve-fitting algorithms because blind fits contain no physical correlation. The engineer’s insights in outlier removal are preferable to perceived advantages in time saved. Whereas Examples 1-6 validate kh and kv predictions through agreement with input permeabilities to the Forward Simulator, the present application illustrates actual field situations where our foregoing synthetic data is now replaced by real data. This history matching may be performed at any time, early, intermediate or late.

Figure 38 Possible kh-kv solutions are determined by well logging measurements.

3.8.3 Case 3: Anisotropy (kh, kv) from Inverse Simulator #1 Model for “Dimensionless Late Time” Conditions

In Case 2, we used Inverse Simulator #2 with source probe data to define the constraint kh^2/3kv^1/3 = k#. This relation was further constrained using experimental field data to provide kh and kv individually. In fact, we used the entire measured pressure versus time curve at the observation probe at a time value that is arbitrary. Here, we closely examine Inverse Simulator #1. Its underlying derivation uses late time asymptotic expansions for the dimensionless solution. However, dimensional and dimensionless parameters are not identical – they often differ significantly as shown in Section 2. What is late time dimensionlessly may, in fact, be early time dimensionally. The converse is often true in mathematics. From calculus, for example, the early time Taylor expansion sinx ≈ x - x³/3! + x⁵/5! - x⁷/7! + … is valid over -∞ < x < +∞. To study this, we re-run our exact Forward Simulator for a longer duration, extending the earlier 10 sec to 110 sec in Figure 39. The following synthetic pressure results are obtained.

Figure 39 Synthetic pressures for source (left) and observation (right) probes.

Pressures always vary rapidly at source probes and equilibrate, but much more slowly at distant observation probes since diffusion is large. A long time is needed to reach steady-state. Equilibrium is not achieved even for t > 100 sec. Let’s use 100 sec (about 2 min) data for practical purposes to evaluate Inverse Simulator #1.

Run 1: Use Forward Simulator t = 108 sec Data in Inverse Simulator #1.

Please refer to Figure 40 for the menu inputs used in this run.

Figure 40 Inverse Simulator #1 data.

This solution compares favorably with exact Forward Simulator inputs, that is, kh = 10.9 vs 10 md for a 9% error, and kv = 0.842 vs 1 md for 16% error.

Run 2: Use Forward Simulator t = 33 sec Data in Inverse Simulator #1.

Predictions for kh and kv individually are given next. The input data are shown in Figure 41 and results appear immediately thereafter.

Figure 41 Inverse Simulator #1 data.

This solution compares favorably with exact Forward Simulator inputs, that is, kh = 11.8 vs 10 md for an 18% error, and kv = 0.725 vs 1 md for 28% error.

Run 3: Use Late-Time Forward Simulator “Early-Time” t = 11 sec Data in Inverse Simulator #1.

Please refer to Figure 42 for the menu inputs used in this run.

Figure 42 Inverse Simulator #1 data.

This solution compares favorably with exact Forward Simulator inputs, that is, kh = 13.6 vs 10 md for a 36% error, and kv = 0.554 vs 1 md for a 45% error. To log analysts with field experience, permeability predictions having much higher errors are not uncommon.

In conclusion, Inverse Simulator #1 calculations for t = 108, 33 and 11 sec using the exact, late dimensionless time, asymptotic solution for (kh, kv) produce acceptable permeability predictions. This is true even when “early” dimensional t = 11 sec data is used. A logical inconsistency? No, not really. “Early time, 11 secs” relative to human perception and prejudices, affected by drilling rig costs or environmental conditions, may well be consistent with late dimensionless times. After all, “early, intermediate and late” in dimensionless asymptotic expansions are not clearly connected with human or work-perceived physical times. This distinction, while clear in retrospect, is never clearly emphasized in mathematics books.

3.9 Example 8: Additional Very Low Permeability Isotropic Validation

In math modeling, asymptotic expansions operating on dimensionless solutions are used for simplicity. However, these do not define “how early is early” or “how late is late” in terms of physical experience or actual times. We consider a very low permeability isotropic formation to further explore this question. Forward synthetic data was created for an isotropic kh = kv = 1 md run, with inputs and calculated source and observation probe pressures shown in Figure 43.

Figure 43 Forward Simulator inputs, source (top) and observation probe (bottom) results.

From the above exact probe results, steady conditions are rapidly attained at 20 sec, certainly “early” by any practical (non-mathematical) engineering measure. We now use t = 20 sec data in our asymptotic late dimensionless time Inverse Simulator #1 model in Figure 44.

Figure 44 Inverse Simulator #1 data.

This solution compares favorably with exact Forward Simulator inputs, that is, kh = 1.04 vs 1 md for a 4% error, and kv = 0.953 vs 1 md for a 5% error. The accuracy is welcome given that permeabilities of 1 md still represent highly diffusive phenomena in many geological applications.

4. Conclusion and Closing Remarks

We have developed a method to predict both individual kh and kv in uniform transversely isotropic media that is rapid, stable, which also offers pore pressure. We emphasize that the “exact” model assumes a centered formation tester nozzle responsible for spherical flow (in isotropic limits) and ellipsoidal flow (in the anisotropic case). Wellbore effects, e.g., radius, invasion, overbalance and supercharge are excluded. Nozzle details are introduced using a multiplicative “geometric factor” correction as is commonly done. The method presented offers convenience, providing results in seconds, and because computing resource demands are minimal, is suitable for downhole real-time implementation. Our prior works have addressed effects not addressed here, but have not offered kh and kv individually. For example, borehole effects as noted are addressed in Chin [26] to include pressure and contamination coupling in the miscible multiphase limit. Overbalanced and underbalanced borehole pressures are studied in Chin [25] although not specifically with respect to packer applications.

Our highly accurate forward simulator calculates pressures versus time at source and observation probes, plus all field locations, useful for logging job planning and for use in validating inverse k_eff, kh and kv predictions when arbitrarily selected pressure versus time data are used. In general, we select first time point, last and an arbitrary intermediate time point for best results over the intermediate time interval. Theoretical results were provided in complete mathematical detail, and inverse strategies were used in developing capabilities for effective permeability k_eff and individual kh and kv predictions, using data from the first seconds to minutes of FTWD well logging. A comprehensive set of validation examples was given. The methodology, applicable to all dip angles in transversely isotropic media, is now available to the petroleum industry for reservoir characterization applications (use of a 45 deg dip example ensures that all forward and inverse software loops are executed properly). Our modeling approach adds to broad industry capabilities in formation evaluation and should make workflows run more smoothly, cost-effectively and productively.

Formation testers used commercially take numerous forms, e.g., Figure 45 displays a subset of present commercial offerings. To account for differences in geometric detail, finite difference curvilinear grid methods of finite element methods such as those used in Figure 1 are needed. While they provide qualitatively useful results, CFD methods are highly grid dependent: change the grid, and the answer changes. Moreover, such methods introduce “artificial viscosity” errors due to round-off and numerical truncation effects, and thus affect permeability predictions. For these reasons, geometric idealizations such as those used here are necessary, but also because they provide rapid and stable predictions useful for trend analysis. In these approaches, a multiplicative “geometric factor” correction is used. For instance, a source probe nozzle with an inner radius of 1 cm may be corrected to 0.9 or 1.1 cm, a decision made based on experimental data, computational suggestions, or simply physical intuition. This modification to radius is simple and easily performed, but is not unique, nor is the approach rigorous mathematically.

Figure 45 Conventional formation testers with different nozzle geometries. Credit: Google Images.

Again, spherical (isotropic) and ellipsoidal (anisotropic) approximations are used for simplicity and the need to avoid completely computational solutions. Thus, cylindrical well effects like those in Figure 46 are represented by simpler ideal physical domains. This is not unreasonable since, in the distant farfield, the physical problem does not “see” the formation tester, but rather, a “black hole” into which fluid is injected or withdrawn. Thus, practical situations like those in Figure 46 can also be solved with such representations. A problem not addressed in this paper are problems associated with pressure overbalance and underbalance in the well, in the former case, associated with mudcake growth, fluid invasion and supercharging. Such problems require a change in initial conditions used in the boundary value problem formulation, and hence, a completely different solution strategy and process. This problem has also been addressed, but because this differs from the model at hand, will be presented in a separate but forthcoming publication. This concludes our presentation on kh and kv methods.

Figure 46 Formation testing in high overbalance pressure environment, associated with dynamic mudcake growth, fluid invasion and supercharging.

Author Contributions

Wilson C. Chin, Ph.D., MIT and M.Sc., Caltech, contributed solely to the authorship and methods reported in this paper. He is author of six formation testing books published by John Wiley and Sons, New York and Elsevier Science, Amsterdam, and developer of real-time permeability prediction software used by Halliburton Energy Services (GeoTap™) and China Oilfield Services Limited (EFDT™). In 2005, Wilson was awarded two national awards by the United States Department of Energy’s “Small Business Innovation Research” (SBIR) program to pursue formation testing methods development. His petroleum interests also include MWD, reservoir engineering, electromagnetic logging, and drilling and cementing rheology.

Funding

Funding for the work reported here was supported entirely by Stratamagnetic Software, LLC.

Competing Interests

The author declares that there are no competing interests or conflicts of interest.