[SOLVED] MOD510 Project 2-Transport in Porous Media

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).

A final tip.

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 time- dependence 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

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

u = − k ∇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:

Q = kA · ∆P , (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).

3

• 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.

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:

∂ (φρ)=−∇·(ρu), (3) ∂t

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

3

• 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:

  cf = − 1 · ∂V = 1 · ∂ρ . (5) V ∂p ρ ∂p

  Similarly, we want to replace the derivative of porosity by using the rock com- pressibility:

  cr = 1 · ∂φ . (6) φ ∂p

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

  ρφct · ∂p = ρ k 􏰅cf ∇p · ∇p + ∇p2􏰆 . (7) ∂t μ

  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:

  ρ k cf ∇p · ∇p ≈ 0 . (8) μ

If we use this approximation, and if we further define
η= k , (9)

μφct we end up with the diffusivity equation.

4

The diffusivity equation.

The diffusivity equation describes how pressure waves travel in a porous

rock.

∂p =η∇2p. (10) ∂t

The factor η can be consider a diffusion constant, in analogy with molecular diffusion. Whenever you have a diffusion equation, you can use the relation

length2 diffusion constant ≈ 4 · time .

to give an order of magnitude estimate the speed of the “diffusion”. For example, if η = 2.5 m2/s and the reservoir radius is 1 km, it will take approximately

(1000 m)2
4 · 2.5 m2/s ≃ 1.2 days

for the pressure wave to reach the outer boundaries of the reservoir.

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

∂p 1 ∂ 􏰀 ∂p􏰁
∂t=ηr∂r r∂r . (11)

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:

∂p|r=re =0. (13) ∂r

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

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

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:

and

􏰀2πhk ∂p􏰁
lim r =Q, (14) r→0 μ ∂r

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

The exact solution to equations (11), (14), and (15) is2 Qμ 􏰀r2􏰁

p(r,t)=pi+4πkh·W −4ηt , (16) where W is the exponential function, defined in terms of an integral:

􏰇x eu

W(x)=
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:

−∞

u du. (17)

2You do not have to show this.

6

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.:

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!

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

7

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

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:

8

# 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 spe- cific 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):

Part 2.

• Plot the line source solution as a function of r assuming default model input parameters (i.e., at t = 10 days).

9

• What do you see? Does the curve approach the expected limits?
Part 3. Next, we are going to use a clever coordinate transformation. Define

y(r) = ln r , (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 , (19) ∂t rw2∂y2

∂p|y=yw = Qμ , (20) ∂y 2πhk

∂p|y=ye =0. (21) ∂y

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:

Part 4.

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

∂t ∆t ∆t
We have introduced the following short-hand notation: pn+1 ≡ p(yi, t + ∆t), and

i
pni ≡ p(yi , t). Next, we are going do discretize the spatial (right-hand side) term

10

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 pn+1 + pn+1 − 2pn+1
= i+1 i−1 i +O(∆y2) (23)

11

∂y2 ∆y2

Numerical scheme for interior grid points.

By neglecting higher order terms, the finite difference discretization of equation (19) becomes

0 0

􏰊

−ξ3

1 + ξ3

pn+1 pn3 3

􏰋􏰊 􏰉􏰈 􏰋 􏰊 􏰉􏰈 􏰋

pn+1 − pn
i i =η·

pn+1 + pn+1 − 2pn+1
· i+1 i−1 i (24)

e−2yi rw2

four points, i = 0, 1, 2, 3.
Note that some of the equations contains p−1 and p4, both of which are outside

of the physical simulation domain!

Part 6.

• 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 0 pn+1 pn −ξ β 00000

−ξ 2+ξ −ξ 0 pn+1 pn
1 1 1 1=1,(25)

22222  0 −ξ 2+ξ −ξ pn+1  pn 

∆y2
• Suppose that we use N = 4 grid points. Write down equation (24) for all

∆t

where we have defined

and

􏰉􏰈

A

pn+1 pn−d ξi ≡ ηe−2yi ∆t ,

r w2 ∆ y 2

β≡ Qμ∆y. 2πkh

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.

12

6 Exercise 3 Implement the full numerical solu- tion

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:

p0 =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:

1. A function setup_matrix, in which you store the coefficient matrix as a

   regular NumPy matrix4.

2. 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.: scipy.sparse.diags.
3. 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.

4An N × N NumPy array.

13

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.

  14

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.

8 Guidelines for project submission

The project has to be handed in as a notebook, though you can submit an additional PDF if you want. You should bear the following points in mind when working on the project:

• Start your notebook by providing a short introduction in which you outline the nature of the problem(s) to be investigated.
• End your notebook with a brief summary of what you feel you learned from the project (if anything). Also, if you have any general comments or suggestions for what could be improved in future assignments, this is the place to do it.
• All code that you make use of should be present in the notebook, and it should ideally execute without any errors (especially run-time errors). If you are not able to fix everything before the deadline, you should give your best understanding of what is not working, and how you might go about fixing it.
• Functions that you code should be reusable, never copy and paste code. Copy and paste is a strong indication that you instead should define function.
• If you use an algorithm that is not fully described in the assignment text, you should try to explain it in your own words. This also applies if the method is described elsewhere in the course material.
• In some cases it may suffice to explain your work via comments in the code itself, but other times you might want to include a more elaborate explanation in terms of, e.g., mathematics and/or pseudocode.
• In general, it is a good habit to comment your code (though it can be overdone).
• When working with approximate solutions to equations, it is always useful to check your results against known exact (analytical) solutions, should they be available.
• It is also a good test of a model implementation to study what happens at known ’edge cases’.
• Any figures you include should be easily understandable. You should label axes appropriately, and depending on the problem, include other legends etc. Also, you should discuss your figures in the main text.

  15

• It is always good if you can reflect a little bit around why you see what you see.

References

1. [1]  M. King Hubbert. Darcy’s law and the field equations of the flow of under- ground fluids. Transactions of the AIME, 207(01):222–239, 1956.
2. [2]  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.
3. [3]  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.
4. [4]  AM Sulak and J. Danielsen. Reservoir aspects of Ekofisk subsidence. In Offshore Technology Conference. Offshore Technology Conference, 1988.

[5] Luigi Tosi, Pietro Teatini, and Tazio Strozzi. Natural versus anthropogenic subsidence of venice. Scientific reports, 3(1):1–9, 2013.

16