MOD510- Project 2: Transport in Porous Media Solved

Transport in Porous Media - From

Laboratory Scale to Field Scale

Learning objectives.                 By completing this project, the student will:

Learn to use matrix equation to do linear regression, and estimate model parameters from experimental data
Solve the diffusivity equation using an implicit solution method
Compare sparse matrix solvers with standard dense solvers
Estimate reservoir properties, such as reservoir permeability, and reservoir volume using a model for the well pressure and data from a well test
Abstract
In this project we want to show you how we extract information from both lab and field data. As many of you do not have a proper background in the study of transport in porous media, we have included thorough introduction and derivation of the necessary equations. If you want to get started quickly with the coding, without having to read and understand all the theory first, we suggest that you jump straight to equations (19), (20), and (21). Show how you can derive these equations from the earlier equations involving the r-coordinate, equations (11), (12), and (13). Next, implement a solver that solves the system iteratively. Start from a known initial pressure pi, and then compute the pressure distribution at a sequence of times ∆t, 2∆t, 3∆t etc. After you have finished implementing the solver, you can go back and do exercise 1, and then 3 and 4 (or 3, 4 and then 1).

If you want to copy and paste some of the provided example code, it is best to use the html version of the project. If you copy from the pdf, some of the symbols might not be pasted correctly.

1           Introduction
Reservoirs are large pieces of porous rocks (usually sediments) that are buried underground. The size of a reservoir may vary a lot: typically by as much as 100-1000 m horizontally, and 10-100 m vertically. A sedimentary rock consists of many small grains, like sand on a beach. During geological time the sand grains have become cemented together, but still there is plenty of room for oil, water and gas to occupy the pores between the grains. Usually 20-30 % of a rock is void space. The absolute permeability is a measure of how easily fluids flow through the rock. It plays the same role as does the inverse of resistance in electromagnetism. For example, a high permeability means that fluids flow easily through the rock, just as a low resistance conductor transmits electric current easily. What actually drives the fluid flow are pressure differences, and if we continue the analogy with electricity, these pressure differences play the same role as voltage gradients inside a conductor.

When a reservoir is opened for production, the pores close to the producer might clog due to the mobilization of particles, and because of chemical reactions taking place etc. [2]. If this happens, the permeability will be reduced. To monitor changes in reservoir permeability, a key tool is well pressure testing [3], which is the topic of this project. It turns out that by studying the timedependence of well pressures during production, one can actually learn a lot about the properties of a reservoir. This information is clearly extremely relevant for an oil company, but it of equal importance in geothermal applications, as well as during groundwater monitoring. Reservoirs are huge (∼ km size), and the well is only about half a foot, but still the information coming from a single producing well can provide valuable information about, e.g., the size of the reservoir, or about formation damage close to well. If one finds that the reservoir has become clogged, one has to implement some kind of intervention, after which one can measure the effect of the intervention by repeating the well test.

Another factor that may be important to consider is reservoir compaction. As fluids are produced, the pressure in the formation drops, and this could lead to large increases in the effective stress in the reservoir, and to seabed subsidence. An example of this is Ekofisk, Norway’s largest oil field, where platforms were observed to sink by several meters [4]; even if the reservoir itself is located more than 3 km below the seabed! A similar process has been taking place in Venice, which is currently sinking. Part of the recorded subsidence of the city has been attributed to groundwater pumping operations [5].

2           Exercise 1
Darcy’s law [1] is frequently used to calculate the relationship between flow rate and pressure when a fluid phase is moving inside a permeable medium. If we neglect gravity and if we use a consistent set of units (e.g., SI-units), Darcy’s law can be expressed in coordinate-free form as

k

u = −    ∇p,         (1) µ

where u[L/T] is the Darcy velocity, k[L2] is the absolute permeability tensor, µ[ML/T] is the fluid viscosity, and p[M/LT2] is fluid pressure.

Part 1. Consider a homogeneous cylindrical core sample of length L = 7 cm, diameter d = 3.8 cm. Imagine that we have first saturated the sample with water, and that we subsequently perform a sequence of injection tests in which water is pumped through the core at a constant volumetric flow rate, Q. For each flow rate, we record the (steady-state) pressure difference across the core, ∆p. Darcy’s law then implies that:

kA      ∆P

Q = · , (2) µ              L

where A = πr2 is the cross-sectional area of the core.

Read the data stored in the file core_flood.dat. Make a scatter plot of ∆P versus Q.
Part 2. At room temperature, the viscosity of water is approximately 1 cP (or 10−3 Pa · s in SI-units).

Use linear regression to estimate the slope of the data, and estimate the permeability from equation (2). If you use the given units (atm, ml/s, cm, cP) the permeability, k, will have unit of Darcy (1 D ' 10−12 m2). When doing the linear regression formulate the problem as a linear problem (see example in the lecture notes), and use a matrix solver do find the regression coefficients.
3           Theoretical background for the rest of the project
The time-dependent pressure development in a reservoir can be studied by setting up a fluid mass balance:

,                                                    (3)

where ρ[M/L3] is the fluid density, and φ is the porosity; the void fraction of the reservoir that permits fluid flow. We are going to assume that

Darcy’s law accurately describes fluid motion,
gravity and capillary forces can be neglected,
viscosity is independent of pressure,
the permeability tensor is isotropic and homogeneous, so that it can be represented by a single constant value, and
fluid compressibility is low (i.e., we do not consider gases)
Inserting Darcy’s law, equation (1), and then expanding both sides of (3), yields:

φ∂ρ + ρ∂φ = k ∇ρ · ∇p + ρ∇2p . (4) ∂t ∂t µ

The left side contains the fluid density, ρ. We want to replace the density with pressure. This can be done by using the fluid compressibility:

1     ∂V        1    ∂ρ

cf = −        ·          =      ·        .                                               (5)

V      ∂p        ρ    ∂p

Similarly, we want to replace the derivative of porosity by using the rock compressibility:

cr .                                                            (6)

Both cf and cr are measured in units of inverse pressure (LT2/M). Equation (4) can therefore be written as

ρφct  cf∇p · ∇p + ∇p2 .                                          (7)

where ct = cf + cr is the total compressibility. Since we have assumed our fluid to be relatively incompressible (this is the case for water), it is often a very good approximation to drop the term proportional to cf on the right-hand side:

ρ cf∇p · ∇p ≈ 0. µ

If we use this approximation, and if we further define
(8)
k

η =   ,
(9)
k µφct

we end up with the diffusivity equation.

3.1          Diffusivity equation in radial flow

If we use cylindrical coordinates, and if we assume perfectly radially symmetric flow, the diffusivity equation reduces to

= η             r          .                                                 (11)

∂t           r ∂r        ∂r

To solve this equation we need some boundary conditions. Later when doing numerical calculations, we shall assume that there is a vertical well placed in the middle of the reservoir, with well radius r = rw. The well is set to produce fluids at a constant flow rate, Q > 0, in other words:

2πhkr ∂p

|r=rw = Q.     (12) µ ∂r

Far away from the well, at a distance r = re, we shall assume that there is no flow of fluids:

.                                                         (13)

Note that equation (12) is a special case of Darcy’s law, equation (1); this can be seen by inserting ∇ = ∂/∂r, and Q = −uA = −u · 2πrh [1].

3.2             Analytical solutions to the radial diffusivity equation?

The diffusivity equation is extremely hard to solve exactly, and explicit formulas can be found only in idealized cases. One such solution is called the line source solution. It employs slightly different boundary conditions than formulated above:

It considers the well to be infinitely small (i.e. a line), thus neglecting the influence of a finite well radius, and
The reservoir is assumed to be infinite. Moreover, rather than assuming no flow far away from the well, the line source solution assumes a constant pressure pi, the initial reservoir pressure, at infinity.
In mathematical terms, the boundary conditions for the line source solution are:

2πhk ∂p

lim r = Q, (14) r→0 µ ∂r

and

lim p(r,t) = pi .         (15) r→∞

The exact solution to equations (11), (14), and (15) is[2]

Qµ  r2  p(r,t) = pi +  · W −  , (16)

4πkh                 4ηt

where W is the exponential function, defined in terms of an integral:

Z x        eu

W(x) =               du.                                                 (17)

−∞ u

In Python this function can easily be computed by using the sp.special.expi function.

4           A Crash Course on Classes
You are free to code the project as you see fit, but regardless of how you do it you should strive to a) write clear, concise code and b) reduce unnecessary code repetition. One way to achieve that is to use classes, because classes provide a mechanism for wrapping parts of your code into objects, which you can afterwards reuse in several places.

If you want to use classes for this project, there are really only a couple of things you need to know. First, all of your classes should include a special function called __init__, in which you declare the variables (attributes) you wish an instance / object of the class to keep track of. Second, all of these attributes should start with the keyword self, followed by a dot and then the actual attribute name. Third, functions inside a class should have self as the first function argument (when defining the function). All of this is best understood via an example:

class Person:

"""

A class representing custom objects of the type Person.

""" def __init__(self, name, age):

self.name = name self.age = age self.species = ’Human’

def greeting(self):

print(’Hello, my name is {}.’.format(self.name))
This first few lines tells Python that in order to create a new object of the type Person, you have to provide the name and age of the person as input arguments, for example:

p1 = Person(’John’, 36)

When executing the above code, Python automatically invokes the __init__ method, and sets the object attributes self.name and self.age based on the user input. In addition, we see that the object p1 has a third attribute, species, which is hard-coded inside __init__, and is therefore not needed as input (this is a good way to set constants that you expect to re-use across parts of your program).

Once we have created a particular Person object, we can call the other function we have defined, namely greeting:

p1.greeting()

Note that we do not include self as an input argument when calling the function; only when we define it [3] If you want to access variables belonging to an instance of the class, you can also access it via "dot-notation", e.g.:

print(p1.name) print(p1.age)

Why we suggest using classes.

The main advantage of using classes in this project is that you can define all the numerical model parameters once and for all within __init__, and then simply re-use them inside all your other functions.

However, you are not required to use classes to complete this project. If it seems too advanced at the moment, just code using ordinary functions. You should still try to make your functions re-usable, though!

5           Exercise 2
You are going to implement a numerical solver for the radial diffusivity equation. To avoid unnecessary complications with units, you may use the following class definition as a starting point:

class RadialDiffusivityEquationSolver: """

A finite difference solver for the radial diffusivity equation. We use the coordinate transformation y = ln(r/rw) to set up and solve the pressure equation.

The solver uses SI units internally, while "practical field units" are required as input.

Except for the number of grid blocks / points to use, all class instance attributes are provided with reasonable default values.

Input arguments:

name                                                                   symbol              unit

--------------------------------------------------------------

Number of grid points                                     N                          dimensionless

well radius              rw ft extent of reservoir re ft height of reservoir h ft permeability k mD

porosity phi dimensionless fluid viscosity mu mPas (cP) total compressibility ct 1 / psi constant flow rate at well Q bbl / day initial reservoir pressure pi psi

constant time step dt days maximal simulation time max_time days

-------------------------------------------------------------

"""

def __init__(self, N, rw=0.328, re = 100000., h=8.0, k=500, phi=0.1, mu=1.0, ct=17.7e-6, Q=1000, pi=2000, dt=0.1, max_time=10):

# Unit conversion factors (input units --> SI)

self.ft_to_m_ = 0.3048 self.psi_to_pa_ = 6894.75729 self.day_to_sec_ = 24*60*60 self.bbl_to_m3_ = 0.1589873

# Grid self.N_ = N

self.rw_ = rw*self.ft_to_m_ self.re_ = re*self.ft_to_m_ self.h_ = h*self.ft_to_m_

# Rock and fluid properties self.k_ = k*1e-15 / 1.01325 self.phi_ = phi self.mu_ = mu*1e-3

self.ct_ = ct / self.psi_to_pa_

# Initial and boundary conditions self.Q_ = Q*self.bbl_to_m3_ / self.day_to_sec_ self.pi_ = pi*self.psi_to_pa_
# Time control for simulation self.dt_ = dt*self.day_to_sec_ self.max_time_ = max_time*self.day_to_sec_ self.current_time_ = 0.

READ THIS FIRST: A note on default model parameters.

Whether you code with classes or not, we shall assume that the default model input parameters are those shown in the above class definition. For example, the default value for viscosity is 1 cP, the default well radius is

0.328 ft, and so on.

Part 1.

Implement a Python function called line_source_solution that returns the analytical line source solution, equation (16), as a function of r and t. The function should work both in the case when r and t are floating-point numbers, and when one of them is a float and the other is a NumPy array.
Implement another function called plot_line_source_solution that makes a plot of the reservoir pressure versus radial distance at a specific point in time. This second function should take the time in question as an input argument, and it should call line_source_solution in order to compute the pressures.
For those of you who choose to code using classes, the functions should be member functions of the solver class; otherwise, use standalone functions with appropriate input arguments and default values.

The following code snippet shows how the above functions might work in practice (if following the class approach):

# Create a solver object. Use only default arguments: solver = RadialDiffusivityEquationSolver(N=100)

# We can calculate the analytical solution at a specific point # in space and time: analytical_sol = solver.line_source_solution(10.0, 42.0)

# ...or at a specific position, but at several times:

time_array = np.linspace(0.0, 50, 51)

analytical_sol = solver.line_source_solution(10.0, time_array)

# ...or at a specific time, but at multiple locations:

position_array = np.linspace(0.1, 10, 11)

analytical_sol = solver.line_source_solution(position_array, 42.0)
Part 2.

Plot the line source solution as a function of r assuming default model input parameters (i.e., at t = 10 days).
What do you see? Does the curve approach the expected limits?
Part 3.             Next, we are going to use a clever coordinate transformation. Define

r

y(r) = ln   ,                                                        (18)

rw

where, as before, rw is the radius of the well.

Show that in the new y-coordinate, equations (11), (12), and (13) can be rewritten as:
∂p           e−2y ∂2p

= η    , ∂t          rw2 ∂y2
(19)
Qµ

 ,                                                        (20)

πhk

.                                                                  (21)

Part 4.   We are now going to discretize equation (19). We start by dividing the total flow domain into N equally sized grid blocks in the y-coordinates, see figure 1 for an illustration. The numerical pressure solution will be calculated at the midpoints (centers) of the y-blocks; note that this does not correspond to the center of the r-blocks!

Make a NumPy array of size N where the elements are the midpoint coordinates of each y-interval. You should give the array a meaningful name, for example: y_centres_ or y_midpoints_.
Create another array holding the r-coordinates of the grid points. Again, this will be a NumPy array of size N (hint: invert equation (18).)
Plot the line solution using your new r-coordinates. Use a logarithmic scale on the x-axis. Comment on the plot. Are the r-points spaced out as they should?
If you code your solver as a class, the arrays should be attributes that are initialized within the __init__ function. However, the plotting should not be included inside __init__.

Part 5. Let i denote an arbitrary grid point inside our simulation domain. To discretize the left-hand side of (19) for point i, we use the following approximation:

∂p = p(yi,t + ∆t) − p(yi,t) + O(∆t) ≡ pin +1 − pni + O(∆t).                       (22)

∂t                          ∆t                                                    ∆t

We have introduced the following short-hand notation: pni +1 ≡ p(yi,t+∆t), and pni ≡ p(yi,t). Next, we are going do discretize the spatial (right-hand side) term

Figure 1: Sketch of the coordinate transformation r → y. Note that the pressure is always evaluated at the center of the y-blocks.

of Eq. (19). However, before we can approximate partial derivatives with respect to y, we have to decide what value to use for the time; after all, the pressure is a function of both space and time, p = p(r,t).

If we select the time t, corresponding to the pressure solutions pni , we say that we use an explicit scheme. On the other hand, if we use t + ∆t, we say that we have an implicit scheme. It turns out that the implicit scheme is much more numerically stable, hence we choose that option here:

∂2p       pni                                                +1

=(23)

∂y2                                     ∆y2

Suppose that we use N = 4 grid points. Write down equation (24) for all four points, i = 0,1,2,
Note that some of the equations contains p−1 and p4, both of which are outside of the physical simulation domain!

Part 6.

|
{z

A
}| {z } pn+1
|
{z pn−d
}
where we have defined

and
 
ηe−2yi∆t

ξi ≡           rw2 ∆y2 ,

Qµ∆y β ≡            .

2πkh
 
 
 
Again, let N = 4. Replace the derivatives in the boundary conditions equation (20) and (21) by the (second order) central finite difference approximation to eliminate p−1 and p4. Show that we now can formulate the equations as the following matrix equation
1 + ξ0

 −ξ1

 0



0
−ξ0

2 + ξ1

−ξ2

0
0

−ξ1

2 + ξ2

−ξ3
0 pn0+1 pn0 − ξ0β

02 pn12+1 =        ppn1n2  ,

−ξ pn+1 

1 + ξ3                pn3+1                              pn3
(25)
If we know the initial pressure, p0, in the reservoir at t = 0, the pressure distribution at t = ∆t , p1, is found by solving equation (25). Schematically:

p1 = A−1 p0 − d , p2 = A−1 p1 − d , p3 = A−1 p2 − d , etc.

6           Exercise 3 Implement the full numerical solution
Once you have understood what goes on in special case N = 4, it is time to tackle the general situation (any choice of N). We assume the pressure is the same everywhere at time zero:

p00 = p01 = p02 ... = pi .

After setting the initial condition, we need to implement a time loop and solve the discretized equations (25) repeatedly for each time step: ∆t, 2∆t, 3∆t, etc. As before, we need to take special care when coding the equations for the first and last grid points, i = 0 and i = N − 1.

Part 1.            Start by making three different Python functions:

A function setup_matrix, in which you store the coefficient matrix as a regular NumPy matrix[4].
Another function setup_sparse_matrix, where only the non-zero elements of the matrix are stored. For ideas on how to do this, see the official SciPy documentation, e.g.: sparse.diags.
A function rhs, which calculates the right-hand side of the linear equation system.
Part 2.   Note that as long as we use a constant ∆t, we only need to set up the matrices once; when initializing the model. On the other hand, the right-hand side vector must be updated each time step.

Finish implementing the numerical solver for the radial diffusivity equation:

Write a "top-level" function that executes all the necessary steps. The simulator should be able to select at run-time whether the dense and sparse matrix solver is going to be used.
Solve the diffusivity equation numerically assuming default model input parameters and N = 100.
Compare your numerical solutions with the analytical line source solution. Check that you get the same result regardless of whether you use a dense or sparse matrix solver.
Part 3.

Increase the number of grid points to N = 1000. Compare the efficiency of the two linear solvers.
7           Exercise 4 Compare with well data
Although we have measured properties of a rock sample in the lab, the reality in the reservoir might be very different. In this exercise we are going to study data from a well test. During a well test, the production engineer starts to produce from the reservoir, while at the same time monitoring how the well pressure changes as a function of time.

Part 1.

Read the data from well_test.dat into a Pandas DataFrame.
Plot the experimental well pressures versus time. Describe what you see.
Part 2. So far we have only calculated the pressure distribution inside the reservoir; the actual observable well pressure is missing! However, we can estimate the well pressure by discretizing equation (21):

Use a first-order finite difference approximation to determine the well pressure
(Hint: Use Taylor’s formula with step-size ∆y/2.)

Part 3.

Further extend your program by calculating, and storing, the estimated well pressure at each time step in your simulation.
Part 4.

Try to manually match your numerical model to the well test data (Hint: use log scale on the x-axis)
To be able to do so, you might have to modify one or both of the following input parameters: k, pi.
Make a plot in which you show 1) the experimental well test data, 2) your numerical well pressure solution, and 3) the corresponding line source solution (i.e., at r = rw)
Part 5. Notice that towards the end of the well test, the pressure starts to rapidly decline. This is actually a sign that the pressure wave has reached the end of the reservoir.

Use the permeability, k, and initial pressure, pi, you found in the previous task, but change re to to match the start of the late decline period.
Use the value you found for re to calculate the volume of fluids in the reservoir.
Part 6.   (OPTIONAL) Use a non linear solver in Python to estimate pi, k, re as a best fit to the data from the well pressure test.

[1] We use the sign convention that Q> 0 for production from the well.

[2] You do not have to show this.

[3] This is to let Python know that the function should be called on a particular object.

[4] An N×N NumPy array.

References
King Hubbert. Darcy’s law and the field equations of the flow of underground fluids. Transactions of the AIME, 207(01):222–239, 1956.
Roland F. Krueger. An overview of formation damage and well productivity in oilfield operations: an update. In SPE California Regional Meeting. Society of Petroleum Engineers, 1988.
John Lee, John B. Rollins, and John Paul Spivey. Pressure Transient Testing. Richardson, Tex.: Henry L. Doherty Memorial Fund of AIME, Society of Petroleum Engineers, 2003.
AM Sulak and J. Danielsen. Reservoir aspects of Ekofisk subsidence. In Offshore Technology Conference. Offshore Technology Conference, 1988.
Luigi Tosi, Pietro Teatini, and Tazio Strozzi. Natural versus anthropogenic subsidence of venice. Scientific reports, 3(1):1–9, 2013.