Mathematical Problems in Engineering
Volume 2018 (2018), Article ID 4570539, 10 pages
https://doi.org/10.1155/2018/4570539

Research Article

A New Way to Calculate Flow Pressure for Low Permeability Oil Well with Partially Penetrating Fracture

Xiong Ping,1,2 Liu Hailong,3,4 Hu Haixia,² and Wang Guan⁵

¹College of Earth Sciences, Yangtze University, Wuhan 430100, China
²College of Engineering and Technology, Yangtze University, Jingzhou 434020, China
³China Petrochemical Exploration and Development Research Institute, Beijing 100083, China
⁴Key Laboratory of Marine Oil & Gas Reservoirs Production, Sinopec, Beijing 100083, China
⁵Department of Geosciences, The University of Tulsa, Tulsa, OK 74104, USA

Correspondence should be addressed to Liu Hailong; 478277608@qq.com

Received 25 May 2017; Revised 29 August 2017; Accepted 21 March 2018; Published 30 May 2018

Academic Editor: Luca Heltai

Copyright © 2018 Xiong Ping et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

In order to improve the validity of the previous models on calculating flow pressure for oil well with partially perforating fracture, a new physical model that obeys the actual heterogeneous reservoir characteristics was built. Different conditions, including reservoir with impermeable top and bottom borders or the reservoir top which has constant pressure, were considered. Through dimensionless transformation, Laplace transformation, Fourier cosine transformation, separation of variables, and other mathematical methods, the analytical solution of Laplace domain was obtained. By using Stephenson numerical methods, the numerical solution pressure in a real domain was obtained. The results of this method agree with the numerical simulations, suggesting that this new method is reliable. The following sensitivity analysis showed that the pressure dynamic linear flow curve can be divided into four flow streams of early linear flow, midradial flow, advanced spherical flow, and border controlling flow. Fracture length controls the early linear flow. Permeability anisotropy significantly affects the midradial flow. The degree of penetration and fracture orientation dominantly affect the late spherical flow. The boundary conditions and reservoir boundary width mainly affect the border controlling flow. The method can be used to determine the optimal degree of opening shot, vertical permeability, and other useful parameters, providing theoretical guidance for reservoir engineering analysis.

1. Introduction

Low permeability carbonate reservoirs generally have the characteristics of large thickness and natural fractures. It is well known that the perforation completion is used in well completion, and the hydraulic fracturing is used to increase the well production. However, the degree of perforation is relatively small (a few meters), which can affect the wellbore’s pressure transmission and production wells’ productivity [1–3]. Therefore, it is crucial to analyze the seepage pressure of the fractured straight well in low permeability carbonate reservoirs. Several studies have been done on the partial (completing) degree of perforation for homogeneous (heterogeneous) reservoirs. In 1963, Warren built an expression of the seepage pressure for homogeneous reservoirs [4]. In 1973, Ramey and Gringarten used Green's function method to obtain the analysis solution of the partially penetrating straight well fractures [5, 6]; In 1987, Heucman and Abbaszaden analyzed the seepage pressure dynamic characteristics of the infinite large reservoirs [6]. In 1991, Raghavan and Ozkan obtained the point source solution to the partial seepage pressure [7]. In 2000, Bui created the pressure point source solution for the double medium homogeneous reservoirs [8, 9]; Valko and Amini presented a method of distributed volume sources to investigate a horizontal well with multiple transverse fractures [10]; Lin and Zhu developed a slab source method to evaluate performance of horizontal wells with or without fractures [11]. Chinese scholars also contributed a lot in this field. Feng Wenguang and Ge Jiaqing conducted a detailed study on single medium and double medium of non-Darcy flow and draw a map of dynamic pressure characteristic curve [12]; Fanhua and Ciqun obtained the unsteady seepage pressure expression with considering start-up pressure gradient [13]; Liu Qiguon and Wang Jianping used the boundary element method to analyze the problem of unstable reservoir seepage in complex reservoirs [14, 15]; Zhao, Y. obtained an analytical solution of the seepage pressure in Lattice space domain with taking reservoir heterogeneity into consideration [16, 17].

However, hydraulic fractures are assumed to be fully penetrating the formation in the previous studies. Limited efforts have been made to investigate the effects of partially penetrating fracture height on the performance of wells. In practice, fully penetrating fractures may lead to an early or immediate water or gas breakthrough in a reservoir with bottom water or gas cap in contact, whereas partially penetrating fractures provide a better way to prevent the early breakthrough. What is more, most of the previous methods are based on Gringarten and Ramey’s point source solution and Green function, whereas the original physical model established by Gringarten and Ramey only considered the upper and lower bounds, limiting the scope of the application. Therefore, it is necessary to study flow pressure for low permeability oil well with partially penetrating fracture.

This paper first presents the physical model and the mathematical model of the unstable seepage flow in the three-dimensional anisotropic rectangular reservoir under certain conditions, and then the solution of the mathematical model is obtained, which provides a new way to calculate flow pressure for low permeability oil well with partially penetrating fracture.

2. Model Establishing

The reservoir formation is fractured by hydraulic fracturing to form a plurality of fractures as shown in Figure 1. After hydraulic fracture, the half-length and height of the fracture constitute the fracture surface. Combining with the distance between the fractures, the one-dimensional diffusion volume of the reservoir fluid is built. The diffusion volume affects the flow rate of the reservoir fluid in the longitudinal direction. The fracture width and the fracture surface form the range of the linear flow in the fractures, which affects the whole straight well production. The influence of the fracture parameters on the seepage pressure can be analyzed according to the coupling model of the reservoir-fracture-wellbore, and through sensitivity analyzing, artificial fracture parameter combination can be obtained, which forms the optimal fracture grid to achieve the purpose of yield optimization. Therefore, in order to study the seepage pressure of fractured straight well, the seepage pressure of single fracture has to be studied, as shown in Figure 2, and then by using the superposition principle, the multifracture seepage pressure can be obtained. In order to simplify the physical model, the assumptions are made as follows:

(1)The oil production is a constant, and the formation of the reservoir is bounded and nonhomogeneous, with equal thickness.(2)Before producing, the reservoir pressure is equal to the original formation pressure.(3)The reservoir fluid is microcompressible, occurring single-phase and unstable seepage.(4)The intercrack interference is ignored, and the fluid flow in the fracture obeys the Darcy law.(5)The gravity of the fluid and capillary force is ignored, and the porosity and fluid viscosity is constant.(6)The fracture is partially penetrating the formation, and the reservoir fluid flows to the wellbore in a limited range.(7)Crossflow between the fracture and the matrix is ignored, and the fracture has infinite diversion capability.

Figure 1: The schematic of multistage fracturing vertical wells.

Figure 2: Single fracture schematic of section penetration.

The center of one single fracture is located at (); based on the above assumptions, the mathematical model of the fracture is established as follows:

where , , and are in the permeability in the , , and direction. is total compressibility, is porosity, is reservoir pressure, is the time, is formation volume factor, and is the flow rate of per unit area flowing through the fracture.

is the position of the source (sink) phase, and its expression is written as follows:

where is Dirac Delta function.

The initial conditions are written as follows:

The inner boundary condition is written as follows:

The outer boundary condition is written as follows:

3. Model Solving

3.1. Solution Research

According to the literature survey, there are mainly three methods to solve equation (1). First, by ignoring the pressure propagation time and using Green’s function method, the analytical solution of the Lagrangian space to (1) can be obtained [10, 11, 18]. Secondly, on the condition that the boundary extension is regarded as the function of time, by using steady state successive replacement method and differential discretization method, the numerical solution to (1) can be obtained [13, 16, 19]. Thirdly, by using numerical approximation method and series expansion method, the correlation between the leading edge of the pressure propagation and the time can be solved. A new way to solve (1) is present in this paper; that is, by means of mathematical methods such as dimensionless transformation, Laplace transform, Fourier cosine transformation, and separation variable method, (1) can be solved.

3.1.1. Dimensionless Transformation

The dimensionless transformation is a method of converting the seepage equation into a conventional mathematical equation. By dimensionless transformation, the number of comparisons can be greatly reduced, which makes the mathematical physics equation simple, neat, and easy to analyze and solve and help to check and verify the seepage equation [20]. The following dimensionless transformation is introduced in this paper:

where is the length of the firing reservoir [m].

So (1) can be written as

where is

Initial conditions are as follows:

Inner boundary conditions are as follows:

Outer boundary conditions are as follows:

3.1.2. Equation Solving

Laplace transform can eliminate the partial derivative of time from the unstable seepage equation and has been widely used to solve the problem of unstable seepage [20]. By using Laplace transformation, (7) can be written as

By using Fourier cosine transformation of , , and , (12) can be written as

where , , and are the solution of the following equations:

The Laplace space solution is obtained by using Fourier cosine inverse transformation. The Fourier cosine inverse transformation can be written as

where is the norm and its expression is

By using Fourier cosine inverse transformation of , , and , (13) can be written as

By introducing two equations,

(17) can be written as

where

Equation (19) is an analytical solution model under Laplacian space. Its applicable conditions are an anisotropic homogeneous rectangular reservoir with impermeable region and outer boundary closed by fracturing. By using Stehfest numerical inversion method [21], a numerical solution of the seepage pressure can be obtained. By changing the center of the fracture, the seepage pressure of the fracture at different locations can be obtained.

3.2. Model Validation

The mathematical model established in this paper can solve the pressure of multifracture system. For the sake of simplification, it is assumed that the number of fracture is 10. The numerical simulation is used to calculate the seepage field, and then the pressure value (simulated solution) at different time and different positions is output and compared with the numerical solution of the seepage pressure calculated by Stehfest numerical inversion (this work).

The reservoir E300 module in Eclipse 2011 is developed for fractured heterogeneous reservoirs. E300 is used to simulate the pressure variation around a fractured vertical well in a rectangular heterogeneous reservoir. In order to meet the assumptions of (1), the numerical model is set as follows.

The width and the length of the rectangular heterogeneous reservoir are 1 km, and there is an oil production well in the center of the reservoir, which is showed as in Figure 3(a). A five-point well pattern is used to simulate the production, that is, one production well in the reservoir center and four injectors in the four corners to ensure the production of the production well. By adjusting injection volume, the oil well production under different displacement pressure can be obtained.

(a)

(b)

Figure 3: The geometry information representation of the reservoir.

In order to describe the formation fluid heterogeneity, the triangular network of grid system is used to ensure that each crack at least has 3 grids, which is showed in Figure 3(b). Therefore, the plane is divided into meshes, and the average grid step is 50 meters. As (1) describing a single-phase fluid seepage, only 1 simulation layer is divided in the vertical direction of the reservoir according to the seepage characteristics of single-phase fluid seepage. The simulation of the total number of grid . The required parameters for the numerical simulation are shown in Table 1, and the results of two methods are shown in Table 2. From Table 2, it can be seen that the relative error is under the basic control of 5%, which is in consistent with the allowable error range, suggesting that this method we offered is reliable.

Table 1: Basic data of the system.

parameters
value

saturation pressure
25 MPa

oil viscosity
1.21 mPa·s

oil density
0.79 g/cm³

water compressibility
4.9 × 10-4 MPa⁻¹

oil volume coefficient
1.21 m³/m³

porosity
0.12

injection pressure
46 MPa

effective thickness
5 m

formation temperature
158°F

water viscosity
1.6 mPa·s

dissolved gas and oil ratio
22.31 m³/m³

oil compressibility
8.1 × 10-4 MPa⁻¹

rock Compressibility
4.5 × 10-4 MPa⁻¹

permeability
1.2 mD

original formation pressure
27 MPa

injecting water intensity
0.044 m³/(d·MPa·m)

Table 2: The results of comparative table.

production time

Oil production
(, m³/d)
Pressure of this work
(, MPa)
Pressure of E300
(, MPa)
relative error
(%)

30
17.492
22.984
22.249
3.201

60
17.202
22.603
21.870
3.244

90
16.951
22.274
21.539
3.298

120
16.729
21.982
21.251
3.323

150
16.532
21.723
20.993
3.363

180
16.355
21.490
20.767
3.369

210
16.194
21.279
20.558
3.390

240
16.048
21.087
20.366
3.421

270
15.913
20.910
20.196
3.412

300
15.770
20.723
19.912
3.913

330
15.627
20.538
19.731
3.929

360
15.486
20.353
19.450
4.436

3.3. Flow Period Division

The fracture center is located on the centerline of the wellbore axis, and the partially penetrating degree is 50%. Moreover, and . There is a fracture in the center of the reservoir. In order to study the trend of the bottom seepage pressure and pressure derivative, the flow division schematic of partial penetration fractured vertical wells is drawn, as shown in Figure 4. The partially penetrating fracturing straight well seepage pressure (pressure derivative) model curve can be divided into four flow periods: A (early linear flow), B (medium radial flow), C (late spherical flow), and D (boundary control flow). The early linear flow is affected by epidermal effect and wellbore storage effect, and progressive analysis shows that, during stage A, the slope of the pressure derivative curve is about 1, and the reservoir fluid is continuously infiltrated to the wellbore, and then the reservoir fluid enters the medium-term radial flow. The higher the degree of partially penetrating is, the longer the medium radial flow is. With the pressure gradually spreading out, before the pressure transmits to the boundary, it is mainly the late spherical flow, and the pressure derivative gradient becomes smaller until it tends to be stable. When the pressure propagates to the boundary, the boundary control flow occurs, and the pressure derivative value changes faster. Finally, the reservoir fluid flows to the wellbore in the form of quasi-steady flow.

Figure 4: The flow division schematic of partial penetration fractured vertical wells.

4. Sensitivity Analysis

Based on the control variable method, the parameters affecting pressure and pressure derivative (template curve) such as fracture orientation, fracture scale, the degree of penetration in the reservoir, permeability anisotropy, reservoir boundary condition, and reservoir scale were analyzed by using the parameters of Table 1.

4.1. Fracture Orientation

In this paper, the orientation of fractures is divided into two aspects, namely, in the reservoir center (as shown in Figure 5(a)) and not in the reservoir center (Figure 5(b)). By setting the coordinates of different fractures, the relationship between the reciprocal pressure, the reciprocal of pressure, and the producing time is shown as Figure 6.

(a)

(b)

Figure 5: Schematic diagram of fracture orientation.

Figure 6: The effect of fracture orientation on template curves.

Figure 6 shows that the fracture orientation mainly affects the late spherical flow stage of the template curve. When the fracture is located in the reservoir center, the pressure is easier to propagate outwards, because in the late spherical flow stage, the reservoir fluid flows to the wellbore in the form of space sinks. The fractures of hydraulic fracturing are the “flow network” established in the reservoir and connect the seepage channel of the reservoir, which increases the seepage area exposed to the reservoir. The more symmetrical the fracture is to the center, the more simple the fluid of the flow channel is, and the easier the fluid flow and pressure transmission are.

4.2. Fracture Length

With different fracturing scale, the length of the fractures is not the same. With different fracture length, the seepage area exposed to the reservoir is not the same, so the pressure transmission trend is different, as shown in Figure 7. Figure 7 shows that the length of the fracture mainly affects the early linear and medium radial flow, especially the medium radial flow. At the same production time, as the length of the fracture increases, the pressure drop becomes slower, and the pressure propagation becomes faster, and finally, the pressure increases significantly. In the early linear flow stage, small fractures tend to produce more pressure drop, which is because under the condition of the same production rate, the larger the size of the fractures is, the bigger the seepage area is. So the wellbore fluid supply capacity is strong, and the time of the early linear flow stage is longer. When the fracture length is the same, the pressure drop velocity tends to be consistent, so the fracture length only affects the duration of the early linear flow, and its effect on the pressure drop is not obvious. Using the progressive analysis method, we can see that the slope of the straight line in this stage is 0.5.

Figure 7: The effect of fracture length on template curves.

4.3. The Degree of Penetration in the Reservoir

Reasonable degree of penetration in the reservoir not only can save the cost of perforation, but also can get the maximum yield. The effect of the degree of penetration in the reservoir on fluid pressure is obvious, as shown in Figure 8. Figure 8 shows that the degree of penetration mainly affects the end of the late spherical flow and the beginning of the boundary control flow. When the height of the reservoir increases, the time of the late spherical flow becomes shorter, and the fluid seepage enters the boundary control flow stage earlier. When the reservoir is completely penetrated, the fluid flows into the boundary control flow stage without going through the spherical flow stage.

Figure 8: The effect of opening shot degree on template curves.

4.4. Reservoir Anisotropy

Since the permeability varies little in the horizontal direction, the permeability in and direction is considered as the same value, and the effect of the vertical permeability anisotropy on the template curve is studied, as shown in Figure 9. Figure 9 shows that the permeability of the vertical anisotropy mainly affects the medium radial flow stage. It is mainly because that, in the actual formation, the greater the vertical permeability is, the larger the probability of fluid flow between the fractures is. On the contrary, it is smaller. Therefore, when the ratio of horizontal permeability to vertical permeability increases, the time of the medium radial flow is longer. Because when the horizontal permeability plays a dominant role, the reservoir fluid mainly flows to the wellbore from the horizontal direction.

Figure 9: The effect of reservoir anisotropy on template curves.

4.5. Reservoir Boundary Conditions

The combination of different reservoir boundary conditions has different effect on the template curve, as shown in Figure 10. Figure 10 shows that the boundary condition mainly affects the boundary control flow stage. For the closed boundary, the value of the pressure derivative tends to be stable, and the progressive analysis shows that the pressure derivative is about 0.01. For the constant pressure boundary, the pressure is continuously propagating and the pressure derivative is gradually reducing due to sufficient supply of external energy until the pressure is equal to the boundary pressure. The progressive analysis shows that the pressure derivative slope is about 1, and the reservoir fluid makes quasi-steady state seepage. When the upper (lower) boundary is closed and the lower (upper) boundary is constant pressure, then the pressure change is between the above two cases, and only the pressure derivative slope is not a constant.

Figure 10: The effect of reservoir boundary conditions on template curves.

4.6. Reservoir Width

By setting different reservoir widths, the influence of reservoir width on the template curve is studied, as shown in Figure 11. Figure 11 shows that the reservoir width mainly affects the time when the fluid enters the boundary control flow stage and the speed of the pressure derivative curve decreasing at this stage. With the reservoir width decreasing, the seepage pressure enters the boundary control flow stage more quickly, and the pressure drop is smaller. The smaller the width of the reservoir is, the more quickly the pressure spreads to the boundary under the condition of constant production producing. The seepage pressure must be into the boundary control flow stage earlier, and finally the pressure increases earlier and faster, and the production pressure increases faster to maintain fluid flowing to the wellbore which is constant in the unit time.

Figure 11: The effect of reservoir width on template curves.

5. Conclusion

The mathematical model of the unstable seepage flow in the three-dimensional anisotropic rectangular reservoir is deduced by establishing a physical model which is consistent with the actual formation of the nonhomogeneous reservoir. The model considers the impermeable top, bottom and constant pressure, bottom boundary, and other different boundary conditions combined with each other. The numerical solution of the pressure in real domain of the model is obtained by using Laplace transform, Fourier cosine transform and Stephenson numerical inversion method. The calculation results are in good agreement with the numerical simulation, which proves the correctness of the model and the practicability of the method.

The pressure dynamic model curve can be divided into four flow periods: early linear flow, medium radial flow, late spherical flow, and boundary control flow. Different reservoir physical properties and different fracturing construction scale are in varying degrees affecting the seepage pressure. The depth of the fracture mainly affects the early linear flow stage. The permeability anisotropy mainly affects the medium radial flow stage. The reservoir firing degree and the fracture orientation mainly affect the late spherical flow stage. The boundary conditions and the reservoir width mainly affect the boundary control flow stage.

In the early linear flow stage, the pressure and pressure derivative curve are a straight line, the scale of the fracture mainly affects the early linear flow, and small fractures tend to produce larger pressure drop. In the medium radial flow stage, pressure and pressure derivative curve show radial flow characteristics of finite extending wells in infinite extension system, and the pressure derivative is approximately parallel to the abscissa, which is a constant. The size of the impermeable region and the orientation of the fracture determine the time when the medium radial flow occurs and its duration. In the late spherical flow stage, the degree of penetration and the location of the fracture determine the time when the late spherical flow occurs and its duration. The reservoir fluid flow is infinitely close to the quasi-steady state seepage flow, and the rate of the pressure derivative is about 0.5. In the boundary control flow stage, the pressure is affected by the boundary condition type and the reservoir width.

This method can determine the parameters such as optimal degree of penetration and vertical permeability and provide theoretical guidance for reservoir engineering analysis and fracturing process design.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

References

S. J. Al Rbeawi and D. Tiab, “Effect of penetrating ratio on pressure behavior of horizontal wells with multiple-inclined hydraulic fractures,” in Proceedings of the SPE Western Regional Meeting, Bakersfield, Calif, USA, 2012.

S. J. Al Rbeawi and D. Tiab, “Partially penetrating hydraulic fractures: pressure responses and flow dynamics,” in Proceedings of the SPE Production and Operations Symposium, Oklahoma City, Okla, USA, 2013.

O. Alpheus and D. Tiab, “Pressure transient analysis in partially penetrating infinite conductivity hydraulic fractures in naturally fractured reservoirs,” in Proceedings of the SPE Annual Technical Conference and Exhibition, Denver, Colo, USA, 2008.

M. Mills and M. W. Clegg, “Study of behavior of partially penetrating wells,” Society of Petroleum Engineers Journal (SPE), vol. 2054, 1969.

H. Ramey and A. C. Gringarten, “The use of source and Green’s function in solving unsteady-flow problem in reservoir,” Society of Petroleum Engineers Journal, vol. 13, no. 5, 1973.

A. C. Gringarten and H. J. Ramey, “Unsteady pressure distribution created by a single horizontal fracture, and partial penetration, or restricted entry,” Society of Petroleum Engineers Journal (SPE), vol. 14, no. 4, pp. 413–426, 1974.

M. Buhidmal and R. Raghavan, “Transient pressure of partially penetrating wells subject to bottom-water drive,” Journal of Petroleum Technology, vol. 32, no. 7, 1980.

F. J. Kuchuk and P. A. Kirwan, “New skin and wellbore storage type curves for partially penetrated wells,” SPE Formation Evaluation, vol. 2, no. 4, pp. 546–554, 1987.

M. Abbaszadeh and P. S. Hegeman, “Pressure-transient analysis for a slanted well in a reservoir with vertical pressure support,” SPE Formation Evaluation, vol. 5, no. 3, pp. 277–284, 1990.

M. Onur, A. Satman, and A. Reynolds, “New type curves for analyzing the transition time data from naturally fractured reservoirs,” in Proceedings of the Low Permeability Reservoirs Symposium, Denver, Colo, USA, 1993.

T. D. Bui, D. D. Mamora, and W. J. Lee, “Transient pressure analysis for partially penetrating wells in naturally fractured reservoirs,” in Proceedings of the SPE Rocky Mountain Regional/Low Permeability Reservoirs Symposium and Exbibition, SPE Paper No. 60289, pp. 1–8, Denver, Colo, USA, March 2000.

F. Wenguang and G. Jiali, “The non-Darcy flow problem of unsteady state in a single media or dual media,” Petroleum Exploration and Development, vol. 12, no. 1, pp. 56–62, 1985.

L. Fanhua and L. Ciqun, “Pressure transient analysis for unsteady porous flow with start-up pressure derivative,” Well Testing, vol. 6, no. 1, pp. 1–4, 1997.

C. Shiqing, L. Gongquan, L. Tao et al., “Mathematical model and typical curve for calculating effective hole diameter in the low velocity non-darcy flow testing of dual-media reservoir,” NGI, vol. 17, no. 2, pp. 35–37, 1997.

S. Fuquan and L. Ciqun, “Analasis of pressure and production in the deformable porous media,” Petroleum Exploration and Development, vol. 27, no. 1, pp. 57–59, 2000.

L. Qiguo, L. Xiaoping, and W. Xiaoqing, “Analysis of pressure transient behaviors in arbitrarily shaped reservoirs by the boundary element method,” Journal of Southwest Petroleum Institute, vol. 23, no. 2, pp. 40–43, 2001.

L. Qingshan, D. Yonggang, W. Chen et al., “Application of boundary element in unsteady state flow,” Petroleum Exploration and Development, vol. 23, no. 2, pp. 36-37, 2004.

C. Shiqing, X. Lunyun, and Z. Dechao, “Type curve matching of well test data for non-darcy flow at low velocity,” Petroleum Exploration and Development, vol. 23, no. 4, pp. 50–53, 1996.

W. Jianping, W. Xiaodong, and M. Shidong, “Unstable filtration theory on straight well penetrating in various parts,” Petroleum Exploration and Development, vol. 26, no. 3, pp. 65–71, 2007.

W. Xiaodong, The Foundation of Seepage Mechanics, China University of Geosciences Press, Beijing, China, 2006.

H. Stehfest, “Numerical inversion of Laplace transforms,” Communications of the ACM, vol. 20, no. 1, pp. 47-48, 1970.