Python Series – Basic plotting and curve fitting

//Python Series – Basic plotting and curve fitting

Python Series – Basic plotting and curve fitting

Basic Gas MBE with Python

Written and Edited by Debora Martogi

One of the main approaches in estimating reserves as we produce hydrocarbons is to use the Material Balance Equation. The material balance equation (MBE) is simply a method to track the volumetric balance into, out, and stored (select cases) in the reservoir. The general MBE relates the initial pressure/ volume of hydrocarbons in the reservoir to total production volumes and current pressure conditions. There are several types of MBE depending on the types of fluid(s) in the reservoir (i.e., dry gas, wet gas, gas condensate, volatile oil, black oil) and the primary driving mechanism(s) (i.e., solution gas drive, gascap drive, natural water drive, compaction drive) for hydrocarbon production (Dake, 2010). Two of the main MBEs are designed for gas reservoirs (i.e., dry gas, wet gas, gas condensate) and oil reservoirs (i.e., volatile oil, black oil). The MBE for gas reservoir with no significant water influx nor other additional driving mechanisms (i.e., water drive), also known as volumetric depletion reservoir, can simply be defined as follows (Dake, 2010)


Or in terms of PVT (Pressure, Volume, and Temperature), relations can be written as

Where pi and Zi are the initial pressure and gas compressibility factor, p and Z are the current/ considered pressure and gas compressibility factor,  Gp being the cumulative gas production, and G is the initial gas in-place (GIIP). The ratio between Gp and G is called the fractional gas recovery during depletion. When evaluated at abandonment pressure, the ratio is known as the recovery factor. Provided the PVT history during production, we can determine initial gas in-place from Equation 2 by plotting the history of cumulative gas production against the Pressure/ Compressibility factor ratio as shown in Figure 1.


Figure 1 – Abnormally pressured gas reservoir GIIP approximation

Plotting Gas MBE



Plotting the figure above can be done quickly in Python. The data plotted in Figure 1 is based on the Louisiana offshore gas reservoir production data published in SPE 10125 (Ramagost and Farshad, 1981). To start programming, open Anaconda on your device and launch Sypder or any other Python Compiler you are using. When you launched Spyder (see Figure 2), you should have an empty script editor window on the left-hand side and console windows on the right-hand side (i.e., Help console, Variable explorer, Plots, IPPython console, etc.). We will mainly use the Script editor to write our program, the plots console to view the plots we made, and the IPPython console window to see the results we printed.

The first step to creating the plot in Figure 1 is to initialize the library. Two of the libraries we will need are NumPy and matplotlib.pyplot and are written into the editor window on lines 12 and 13, as shown in Figure 2. For convenience, we rename the NumPy library as np and the matplotlib.pyplot library as plt.

Figure 2 – Spyder interface

The next step we need to do is to load the data we need. I simply copied the data from SPE 10125 to a notepad. The information that we will need for this analysis is the p/z and pressure history data. Following is the script to load the data.

Figure 3 – Loading data stored in a text file named SPE_10125_MBE_Gas into Python

Np.loadtxt is a package under the NumPy library to load text – files. Usecols=(1,2) simply means we are extracting the second and third columns of data stored in the txt – file. Python numbering starts from zero, so to extract the first column, you will need to specify 0 (i.e., usecols=(0,2)). Data extracted will be stored as G_z (from the second column) and pz (from the third column) in the order it is called under Usecols. Unpack = true is to transpose the data we have.  

The next step we can do is to plot the data using the matplotlib. To plot the data as a line, we can use the command plt.plot, while for a scatter plot, we can use plt.scatter. We can then simply plot using the lines shown in Figure 4. The resulting plot was shown in Figure 1 as the blue dots.

Figure 4 – Plotting data as a scatter plot

We can format our plot to show the grid, axes label, and adjust the axes limits using the grid, _label, and  _lim packages. The complete formatting of the plot is shown in Figure 5. Line 22 defines a major grid as a continuous line with a line width of 0.5 and red color. Lines 23-24 defines the minor grid as a dashed line with a line width of 0.5 and black in color. We can also specify the horizontal axis label (i.e., line 25 – xlabel) and the vertical axis label (i.e., line 26 – ylabel). Notice, to define Gp with a subscript, we use the $ sign and the _ so that we scripted $G_p$ to show the formatting correctly. The horizontal axis limit (i.e., xlim) and vertical axis limit (i.e., ylim) are provided in lines 27-28.

Figure 5 – Plot formatting script

GIIP Approximation


The volumetric depletion plot, such that shown in Figure 1, is one sample of “history matching” and production “prediction”/”forecasting.” To predict the initial gas in place, the simplest way to do this is to use the simple gas reservoir MBE as shown in Equation 2 and draw a linear (1-D) trendline that fits the production history data shown in blue dots in Figure 1. However, upon a quick glimpse, if we fit a linear line through all the data as shown, we may end up with an overestimated GIIP. This issue can be avoided by looking back into the reservoir’s initial pressure gradient suggested by Ramagost and Farshad in SPE 10125. They reported that the reservoir’s initial pressure gradient is much higher than what is typically observed in the region, making the reservoir abnormally pressured. This anomaly is due to the overpressure condition sustained by the reservoir fluids. To prevent overestimating the GIIP in volumetric depletion reservoirs, we must monitor the pressure gradient and check whether the reservoir is overpressured. We can then determine the GIIP when the production data shows a trend of a normally pressured reservoir.

We can segment our analysis into two segments: abnormal pressure (early data) and normal pressure (later data) to analyze this data. We can identify where the transition is observed, and in our case, it is around Gp = 120 BSCF. We can then determine the apparent GIP from the abnormally pressured data segment and the actual GIIP from the normal pressure data segment. To do so, we can fit several lines to data within each segment. This can be done by hand, too but relies on the user (i.e., user-dependent). Utilizing Python, we can determine the best-fit linear line that plots through the abnormal production history data in Figure 1 using the polyfit package available under the NumPy library.

Figure 6 – Example script to use polyfit and plot trendlines

In this case, we can select to analyze the abnormal pressure segment up to index 8 or dataset 8 (or 7 in Python, see line 31 in Figure 6). We then define the associated Gp and p/z data using the np.array package listed on lines 32-33 (Figure 6). Using the polyfit package (line 34), we can get the trendline information for our data by specifying the following command np.polyfit(x-data, y-data, degree of fitting polynomial). Since we are looking to get a linear trend, we use a degree of a fitting polynomial of 1. The polyfit function’s output is the coefficients of a linear line equation, as shown in Equation 3. Equation 3 is obtained after rearranging Equation 2 with b being the intercept related to our initial pressure/ compressibility factor ratio, m is the slope related to how we will predict our apparent GIP, x is the cumulative production data, and y is the p/z ratio history. We can then replot the data using the linear coefficients obtained through polyfit. The first coefficient relates to the slope and can be called using poly_1[0], and the second coefficient relates to the intercept. We can then use this information to define new ydata, as shown in lines 38 and 43. Plotting can then be performed similarly to describe previously, but this time we will plot a line, so we are using plt.plot instead (lines 39 and 44). The plot of the apparent pressure prediction is labeled as “Apparent” in Figure 1. Finally, we can determine the apparent GIP, G, utilizing the definition of the b-intercept and m-slope (which includes G) and determined as shown in line 47 (i.e., dividing the intercept-b by the slope-m).

You can then follow a similar procedure to determine the GIIP from the normal pressure history data (i.e., data beyond index 8). In this case, the apparent GIP is 654 BSCF, and the GIIP is 504 BSCF. The error observed is about 29.74% if we were to use the apparent GIP.

Earlier studies had also considered a quadratic formulation to approximate the GIIP for abnormally pressured reservoirs (Duggan, 1992). An example of this is provided in Figure 7 using the same data we use previously. We can obtain the quadratic coefficients and the GIIP similarly, but this time by defining a fitting polynomial degree of 2 in the polyfit command. The GIIP we got is 451.91 BSCF, which is about a 10.37% error.

Figure 7 – GIIP approximation using quadratic fitting line

So here are some basic skills on how to perform a basic curve-fitting approach to “history matching”. Curve fitting is not only been applied to MBE, but also other concepts such as wellbore storage, skin calculations, etc. Thus knowing this basic technique might be beneficial to perform some basic estimations.

Next week, we will look into some more advanced analysis on plotting a dimensionless constant rate solution of an infinite acting reservoir. Until then, stay well everyone!


Disclaimer: The Well Log is a non-profit publication aimed “purely” to educate students at Texas A&M University and beyond on information pertinent to the petroleum engineering industry. All articles are written by student volunteers based on information obtained through online sources and SPE publications. If you are the owner of any materials we cited and would like us to remove it from our publications, please contact The Well Log Editors at


Dake, L. P. The Practice of Reservoir Engineering. Rev. Ed. L.P. Dake. Rev. ed., Elsevier, 2001. 

Dr. Blasingame’s PETE 324 Lecture Notes

Ramagost, B. P., & Farshad, F. F. (1981, January 1). P/Z Abnormally Pressured Gas Reservoirs. Society of Petroleum Engineers. doi:10.2118/10125-MS

Duggan, J. O. (1972, February 1). The Anderson “L” -An Abnormally Pressured Gas Reservoir in South Texas. Society of Petroleum Engineers. doi:10.2118/2938-PA

Commands :       Spyder     Grid     subscript and symbol     Axis labels     Array      Polyfit



  1. Keaton Upshaw November 14, 2020 at 4:40 am - Reply

    Thank you for starting this Python series! This is a great beginner’s example – using Python makes performing regressions so much easier… I am curious about what your code was to compute the 2nd order polynomial regression to obtain the GIIP.

    Here is the code I used for the 2nd order polynomial (GIIP was ~ 450 BSCF):

    index_int = 14
    Gp1 = np.array(Gp[0:index_int])
    Pz1 = np.array(Pz[0:index_int])

    mymodel = np.poly1d(np.polyfit(Gp1, Pz1, 2))
    myline = np.linspace(0, 600)
    plt.plot(myline, mymodel(myline))

    Thanks again!

    • The Well Log November 16, 2020 at 5:05 am - Reply

      Howdy, Keaton!

      Thank you for reading and trying it out. Your code is correct. The slight difference in the answers you obtained and the one posted is due to the data source. I obtained a raw txt file from the source website. I noticed the data in the txt file is provided up to the third decimals while the data from the actual paper is rounded to the nearest integer.

      Thank you for checking this. I hope this clears up the confusion.



Leave A Comment

Welcome Ags!

Welcome to the 2020-2021 academic year. Have a blast!