UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Thermo-viscous fingering in porous media and in-situ soil remediation Moyles, Iain 2011

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2011_fall_moyles_iain.pdf [ 994.86kB ]
JSON: 24-1.0080593.json
JSON-LD: 24-1.0080593-ld.json
RDF/XML (Pretty): 24-1.0080593-rdf.xml
RDF/JSON: 24-1.0080593-rdf.json
Turtle: 24-1.0080593-turtle.txt
N-Triples: 24-1.0080593-rdf-ntriples.txt
Original Record: 24-1.0080593-source.json
Full Text

Full Text

Thermo-Viscous Fingering in Porous Media and In-Situ Soil Remediation by Iain Moyles B.Sc., The University of Ontario Institute of Technology, 2009 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Mathematics)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2011 c Iain Moyles 2011  Abstract The removal of artificially and naturally placed oils found in the ground is of paramount concern for environmental reasons and for the extraction of crude oil for energy. While many companies are involved in this soil remediation process, the physics of the underlying problem are not very well understood. We present a mathematical model and analysis of oil extraction in a porous soil based around the injection of water. We first consider a saturation analysis based on conservation of mass between oil and water in a one-dimensional setting and we find that based on certain parameter values, the water-oil interface goes unstable producing viscous finger patterns. We then include the effects of surface tension between oil and water to determine how this affects the growth of such fingers. We conclude that in the limit of small surface tension effects, the results generalize to the original problem, but more importantly, we deduce that there is a scaling which places the formulation into a setting that is invariant with respect to the magnitude of surface tension effects. With this scaling, we notice that the effect of surface tension is to limit the growth of fingers to a maximal wave number and to prevent their formation entirely beyond a certain critical wave number. Finally with the inclusion of temperature, via heated water injection, we see the formation of a dual fingering pattern: one associated with the mass conservation analysis  ii  Abstract of the oil-water interface and one associated with the conservation of energy across an interface where thermal gradients occur. We see that the thermal gradients across the interface where temperature drops induce unstable viscous patterns with a higher wave number than would occur for an equivalent isothermal interface where there was solely a change in viscosity. The thermal gradients also promote fingering development downstream across the classical viscosity differential driven interface but at the expense of lowering the interfacial velocity. It is interesting to note that the change in saturation that occurs across the energy interface is a result of a pseudo-free boundary created by the thermal problem.  iii  Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iv  List of Tables  vi  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Introduction 2 Model  xi  . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  1  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  7  2.1  Volume Averaging Method  . . . . . . . . . . . . . . . . . . .  7  2.2  Water Saturation . . . . . . . . . . . . . . . . . . . . . . . . .  9  2.3  2.2.1  Capillary Pressure . . . . . . . . . . . . . . . . . . . . 15  2.2.2  Non-Dimensionalization . . . . . . . . . . . . . . . . . 19  2.2.3  Cory-Type Permeabilities . . . . . . . . . . . . . . . . 21  Temperature  . . . . . . . . . . . . . . . . . . . . . . . . . . . 22  2.3.1  Volume Averaged Temperature Equation  . . . . . . . 23  2.3.2  Temperature Dependent Viscosity  2.3.3  Non-Dimensionalization . . . . . . . . . . . . . . . . . 37  . . . . . . . . . . . 36  iv  Table of Contents 3 Results 3.1  1D Model in the Absence of Capillary Pressure 3.1.1  3.2  3.3  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42  Linear Stability Analysis  . . . . . . . . . . . . . . . . 50  1D Model With Capillary Pressure . . . . . . . . . . . . . . . 61 3.2.1  Linear Stability Analysis  3.2.2  Numerical Solution . . . . . . . . . . . . . . . . . . . . 75  3.2.3  Asymptotic Analysis . . . . . . . . . . . . . . . . . . . 83  1D Temperature Model  . . . . . . . . . . . . . . . . 70  . . . . . . . . . . . . . . . . . . . . . 99  3.3.1  Stationary Boundary Characteristics . . . . . . . . . . 106  3.3.2  Internal Interface Characteristics . . . . . . . . . . . . 107  3.3.3  Other Cases  3.3.4  Solving The Saturation Problem  3.3.5  Linear Stability Analysis  . . . . . . . . . . . . . . . . . . . . . . . 116  4 Conclusions and Future Work 4.1  . . . . . . . . 42  Future Work  . . . . . . . . . . . . 117  . . . . . . . . . . . . . . . . 121  . . . . . . . . . . . . . . . . . . 138  . . . . . . . . . . . . . . . . . . . . . . . . . . . 144  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148  v  List of Tables 2.1  Parameter Values for Model. . . . . . . . . . . . . . . . . . . . 20  2.2  Parameter Values for Model with Temperature. . . . . . . . . 40  vi  List of Figures 1.1  Porous sand with fluids filling the void space. Image from [39].  2.1  Here we consider a control volume V which is broken up into  3  its volume fractions. The porosity φ makes up the portion of the volume that allows flow and then that is further subdivided into fractions s and o which make up the compartments where water and oil flow respectively. . . . . . . . . . . . . . . . . . . 10 2.2  A 1D model interface contact, Γ, between oil and water. . . . 15  3.1  Characteristic diagram with expansion fan and shock for the 1D problem (3.3) in the absence of capillary pressure. . . . . . 46  3.2  Saturation profile from solution (3.8) with s∗ = 0.1443 and t = 1/8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49  3.3  Saturation steady-state solution for the inner boundary layer centered around η = 0. It is a smooth heteroclinic orbit connecting the steady-state values of the outer hyperbolic problem. 68  3.4  Pressure derivative steady-state solution for the inner boundary layer centered around η = 0. The derivative has the appropriate limiting behaviour to the outer hyperbolic problem.  69  vii  List of Figures 3.5  An example of a discretization of the saturation and pressure eigenfunctions. The thick line represents the actual discretization with dark circles representing nodes and the white circle representing cell-centers. The flow chart above the line shows the symmetry of the cell-centered discretization and how information flows to describe the i + 1/2 cell. . . . . . . . . . . . 76  3.6  Numerical solution of uˆ1 in (3.21) for n = 1. To compute the solution we used 201 evenly spaced mesh points. . . . . . . . . 80  3.7  Numerical solution of w1 in (3.21) for n = 1. To compute the solution we used 201 evenly spaced mesh points. . . . . . . . . 81  3.8  This figure shows the largest value of σ for any transverse wave number, n. All wave numbers after the terminus n on the plot are negative and thus stabilize. . . . . . . . . . . . . . . . . . 82  3.9  Saturation eigenfunction computed numerically (solid line), asymptotically to first order (dot-dash line) and to second order (dash-line) for n = 1. 201 mesh points were used in the solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95  3.10 Pressure eigenfunction computed numerically (solid line), asymptotically to first order (dot-dash line) and to second order (dash-line) for n = 1. 201 mesh points were used in the solution. 96 3.11 Eigenvalue plot for various values of n ˆ = n. . . . . . . . . . . 97 3.12 Visualization of the area change via a moving internal interface in the water saturation problem. . . . . . . . . . . . . . . . . . 104  viii  List of Figures 3.13 Characteristic diagram for saturation when internal interface characteristics emerge into the hot region. Here, solid thin lines represent the constant solution characteristics from the stationary boundaries, dotted lines represent rarefaction characteristics, dashed lines represent the constant solution characteristics from the internal interface, and the dot-dashed line represents the shock. Notice that the internal interface at x = t corresponds to the jump in temperature while the shock solution at x = ξ corresponds to the traditional BuckleyLeverett behaviour. . . . . . . . . . . . . . . . . . . . . . . . . 111 3.14 This is an artificial example showing how the hot emergent characteristic scenario holds. The solid line is the analytical solution while the dashed line is the numerical one. The dotdashed line is the temperature front. . . . . . . . . . . . . . . 113 3.15 Characteristic diagram for saturation when internal interface characteristics emerge into the cold region. The thin solid lines represent the constant solution characteristics from the stationary boundaries. The dotted lines represent the characteristics from the rarefaction wave, the dashed lines represent the constant solution characteristics from the internal interface and the dot-dashed line represents the shock. Like in the hot region, the internal interface at x = t corresponds to the jump in temperature while the shock solution at x = ξ corresponds to the traditional Buckley-Leverett behaviour. . . 115  ix  List of Figures 3.16 Analytic and numerical results for saturation as it depends on temperature. The solid curve is the analytic solution while the dashed curve is the numeric solution. The temperature is plotted as dot-dash in black. . . . . . . . . . . . . . . . . . . . 119 3.17 Linear stability eigenvalue in the absence of capillary pressure with and without temperature effects. The solid line is the eigenvalue without temperature (3.15) and the dashed line is the eigenvalue with temperature (3.49). . . . . . . . . . . . . . 125 3.18 Figure from [19] showing thermo-viscous fingering in glycerin. The snapshots are of different instances in time. . . . . . . . . 135 3.19 Fingering growth rate for various values of non-dimensional temperature gradient, ∆T . The solid line is ∆T = 1, the dashed line is ∆T = 1/3, the dot-dashed line is ∆T = 1/30 and the dotted line is ∆T = 1/300. . . . . . . . . . . . . . . . 136  x  Acknowledgements Firstly, I’d like to thank my supervisor, Dr. Brian Wetton for his endless guidance, support, and patience through the completion of this thesis. He provided invaluable advice both as a leading academic in applied mathematics, and as a trusted representative for academic-industrial relations, promoting the strong application-centered goals of this work. I appreciate the often longer than scheduled meetings so that we had the opportunity to discuss everything in the appropriate detail it needed. I’d like to thank Lee Yupitun, the math graduate secretary, for all her administrative help.  I appreciate the endless support my family provides from across the country and the enthusiasm with which they discuss my accomplishments. They provide a constant source of motivation to adhere to the tasks that I set out to accomplish and a strong faith in my abilities to complete those tasks.  No one has had to be more patient through this process than my loving girlfriend Sonia Soares whom has constantly motivated me to work hard and strive for greatness. I thank her for her love and understanding through this research process and for her extra assistance and time offered towards any task that provided me the opportunity to increase my research efficiency. She  xi  Acknowledgements has supported me through all the positive and negative endeavours associated with every aspect of this research and my life. I thank her profusely for her love and belief in me.  Finally I’d like to thank many of the friends who have provided welcome distractions through various points of the research process. I’d like to thank my research group consisting of Mark Willoughby, Michael Lindstrom, and Ricardo Alves-Martins for the many laughs and conversations, Kai Rothauge for his strong abilities, providing a constant positive motivation to continue to strengthen my own, my friend Laura Peterson for her years of support, and of course the Physics Council, Ryan McVey, Mike Busuttil, and Matt Hennessy, for the brotherhood they provide. Specifically I’d like to thank Matt Hennessy for providing both a strong pride and admiration in his academic endeavours that inspires me to be successful, for reading and editing this thesis, and also for the much needed vacation that allowed me to maximize my efficiency.  Funding for my graduate work was provided through a one-year NSERC CGS-M fellowship, a UBC affiliated award and the Westcoast Energy Inc Jack Davis Scholarship in Energy studies. All of the financial support provided to me is greatly appreciated.  xii  Chapter 1 Introduction The process of soil remediation is designed, through various means, to remove unwanted contaminants present in the ground. These contaminants, often oil and its byproducts, may enter the soil through various methods such as intentional dumping or seepage through pipelines that run deep beneath the earth [45]. Quite often the remediation is needed to satisfy certain mandated environmental conditions or to rezone industrial land for commercial or residential purposes. However, the extraction of oil from the ground is not limited exclusively to environmental concern or rezoning but is also becoming of paramount importance when obtaining crude oil from tar sands. The issue of oil recovery is of such great interest to industry right now that a conference was organized specifically for CEOs and vice presidents of companies associated with the soil remediation process entitled In-Situ Oil Sands Summit 2011. It was a collaborative effort with businesses and engineers to understand and optimize many of the processes and technologies associated with soil remediation [1]. It is this strong industrial application that motivated our study of the remediation problem. The word in-situ refers to the most common oil extraction process whereby a fluid that is immiscible to oil (typically water) is forced into the system at a given flow rate in order to move the oil to a collection site where it may be extracted [27]. It  1  Chapter 1. Introduction is also fairly common to heat the water that is being injected into the ground.  The extraction process via fluid injection is not trivial. Due to spatial variations in the soil and other factors, the contact surface at which oil and water meet may not be planar due to small inhomogeneities or roughness and hence have oscillatory behaviour. Under certain conditions, this behaviour persists throughout time and this phenomena is known as viscous fingering, a process that cannot be neglected in soil remediation studies due to the strong fluid displacement it causes [38]. Fingering is not a mechanism solely associated with soil remediation; it is known to occur in many applications where there is interfacial growth such as the formation of liquid snowflakes in superheated ice ([16]) and it has been studied extensively for many decades experimentally ([32]), analytically ([46]), and numerically ([35]). Homsy in [20] provides a thorough treatment on the history and progress of viscous fingering as it relates to porous media.  In the remediation problem, the solid component, soil, is very porous and the two fluids fill the surrounding void space. Figure 1.1 shows the porous structure of sand with the fluids filling the void space.  2  Chapter 1. Introduction  Figure 1.1: Porous sand with fluids filling the void space. Image from [39]. Due to the oil displacement associated with the viscous fingering phenomena, we postulate that the fingering mechanism in soil remediation has an adverse effect in the isothermal case. Fingering would hinder the movement of oil towards the collection site and thus lower the remediation efficiency. However, we expect that fingering in the non-isothermal case actually promotes soil remediation. As the heated water comes in to contact with the oil, it lowers the viscosity of oil, thus decreasing the effective ratio of the two viscosities. As the viscosity ratio comes closer to unity, the water is able to drive the oil toward the collection site because viscous fingering has been reduced. Therefore, just as the displacing oil in the isothermal case poses a problem, the displacing hot water in the non-isothermal case allows for greater transfer of heat to the oil, ultimately lowering the viscosity in the most efficient manner thus reducing viscous fingering overall. The analytical work in this thesis is a step towards understanding the exact role of heating in the remediation 3  Chapter 1. Introduction process.  Mathematical modelling of the soil remediation problem arguably began with the publication of [6] by S.E. Buckley and M.C. Leverett in 1941, entitled Mechanism of Fluid Displacement in Sands. They, through mass balance equations, derived a simple one-dimensional hyperbolic partial differential equation, now famously known as the Buckley-Leverett equation, which was considered to produce solutions of the movement of water and oil saturation in the soil. Since then, many engineers and scientists have studied the porous displacement problem to varied degrees of success ([32], [47], [46], [34]). However, most studies neglect the application towards soil remediation. Therefore, it is our aim to develop a model whereby the parameters are consistent with the physics involved in the remediation process and to carefully analyze that model by taking advantage of parameter relationships. In what follows we use the porous structure problem foundation and derive models for water saturation and temperature inside the sand based on parameter values that are commonly associated with soil remediation sites. We begin with the classic Buckley-Leverett equation and the viscous fingering that originates from such a formulation, carefully analyzing the requirements for which the fingering is persistent in time. From there, we aim to understand the basic stability of viscous fingering in the presence of small interface dynamics (i.e. surface tension effects) that occur due to the immiscibility of oil and water. Our goal is to establish a more robust result than in previous literature whereby we exploit small parameters through the method of matched asymptotic expansions ([18], [43]) to clearly demonstrate the tran-  4  Chapter 1. Introduction sition in fingering properties between the cases with and without interface dynamics. Through our analysis we are able to determine the correct scaling between the eigenfunctions as well as recover, to first order, the linear eigenvalue without interfacial dynamics. As well, we are able to cast the problem numerically into one that holds regardless of the relative effect of surface tension. This is actually a fairly novel solution formulation as it produces a localized structure that can be exploited to include all possible magnitudes of surface tension.  Next, we introduce temperature into the model via fixed heated water injection to begin to understand how the changes in thermal energy affect fingering stabilities. This is of important concern for two reasons. Firstly, many of the in-situ processes for removing oil, as stated, use a heated flushing system and secondly, viscosities can vary greatly with temperature and thus, the thermal inclusion is an additional mechanism for viscosity gradients. We postulate that the inclusion of thermal effects will promote fingering. This is due to the fact that with a temperature dependent viscosity, oil becomes less viscous in the presence of a high temperature and thus in regards to the cold oil, it acts like the water, in that it is a less viscous fluid trying to displace a more viscous one. With the inclusion of temperature, we show that a dual fingering pattern occurs, one associated with the natural viscosity gradients and a second pattern due to viscosity gradients caused by thermal changes. The stability of these fingering patterns have a strong interdependence that makes them different than the isothermal case. In fact, it is a pseudo-free boundary in the saturation problem created from the solution of the thermal  5  Chapter 1. Introduction problem that drives the saturation changes in the temperature induced fingering pattern.  The thesis is presented as follows, in Chapter 2 we derive our model through the use of partial differential equations and conservation laws. We then present parameters consistent with those used in the remediation process and determine appropriate order of magnitude estimates and dominating physical processes. In Chapter 3 we analyze one-dimensional variants of our model and subject them to transverse perturbations in order to understand the growth of instabilities or viscous fingers. We end with a discussion of the context provided by the results as well as insight into the future improvements of the model.  6  Chapter 2 Model The soil remediation problem involves the elimination of a stagnant viscous fluid (oil) via a possibly heated front of a much less viscous fluid (water) that is immiscible with it. In order to create a tractable model for this problem, we will assume continuum modeling through the use of partial differential equations in order to conserve mass and energy of all constituents in the soil.  2.1  Volume Averaging Method  The idea behind the volume averaging method is that by exploiting structural symmetry in a domain, we can infer macroscopic properties via averages of conservation laws on the microscopic level. It is worth noting here that we use the term microscopic in a rather liberal sense. We mean that in a multicomponent system, we can look at continuous models in small subdomains of the problem each of which contain a single component only and then average over the entire domain to extract global information. If we consider Vi to be a subdomain volume of some domain V then the local volume average of some quantity (vector or scalar) w contained entirely in Vi is defined as, w =  1 V  w dV. Vi  7  2.1. Volume Averaging Method Note this is just a standard volumetric average where we are able to integrate solely over Vi under the assumption that the variable is only defined in that region. We also define another useful quantity known as the intrinsic volume average, w  ∗  =  1 Vi  w dV. Vi  The volume averaging method is essentially a formulation of homogenization theory as both methods average microscopic details into a macroscopic framework: however standard homogenization theory usually capitalizes on the unification of many small length scales into one macroscopic scale. It is important to note that we assume the domain control volume V to be periodic in structure so that each cell is a microscopic blueprint for the global problem. Making this assumption avoids the need for empirical relationships between adjacent control volumes that define unique properties and allows us to study the behaviour of one representative volume and extract information regarding the global problem. Before we derive the effective temperature equation we will require an appropriate divergence theorem for averaged variables. It is not difficult to show ([21],[41],[40]) that the following identities hold, ∇w = ∇ w +  1 V  wn dA,  (2.1a)  Aint  ∇·w =∇· w +  1 V  w · n dA  (2.1b)  Aint  8  2.2. Water Saturation where n is the standard outward normal vector and Aint is the area of interfacial contact surfaces between adjacent subdomains. It will also be important that [21], w ∇w  = w =∇ w  which states that the average value (and gradient of the average value) are constant with respect to the domain they are averaged over.  2.2  Water Saturation  We consider an immiscible two-liquid species interaction in porous media whereby a viscous fluid (oil) is displaced by a relatively non-viscous fluid (water) in soil. To model the equations for saturation of either fluid, we must consider a control element volume V (Figure 2.1) which will be divided into a porous portion containing the two fluids as well as the solid rock matrix.  9  2.2. Water Saturation  s  (1 − φ)  φ  o  Figure 2.1: Here we consider a control volume V which is broken up into its volume fractions. The porosity φ makes up the portion of the volume that allows flow and then that is further subdivided into fractions s and o which make up the compartments where water and oil flow respectively. We define the fluid phase distribution function as,  a(x) =    1, x is in fluid  0, x is in solid  so that we may define the fraction of a porous media that accounts for void  10  2.2. Water Saturation space, the porosity, as φ=  1 V  a(x) dV =  Vf V  V  where Vf is the fluid volume made up of oil and water. This volume fraction φ, is volume available for flow while the remaining portion (1 − φ) is the solid matrix or non-porous portion. We must further divide φ into subdomain volumes which constitute the fraction of space that will be filled with water (s) and that which will be filled with oil (o). By using a similar phase distribution function as above, we could provide a more rigorous definition of the volume fractions, s and o. Because we have defined our problem so that only water and oil are admissible then for the volume that makes up the porous region, it must be true that s+o=1  (2.2)  and thus it will be satisfactory to derive models for water saturation only. It should be clear from (2.2) that 0 ≤ s ≤ 1. If we consider the volume that only contains water to be Vs and for liquid transport to only be caused by convection1 then we can assume that the continuity equation holds in Vs , ∂ρs + ∇ · (ρs qs ) = 0 ∂t where ρ is the density, q is the volumetric flux (velocity) and the subscript 1  More specifically we assume that any small scale diffusive fluxes are ignored due to the dominance of convection.  11  2.2. Water Saturation s denotes each variable’s association with the water phase. We now take a local volume average of this equation,     ∂ 1 ∂t V  ρs dV  +  1 V  Vs  ∇ · (ρs qs ) dV = 0. Vs  Now there is a subtle interchange of operations here that may seem invalid since Vs = Vs (t). However, we assume that each individual control volume, V , is isochoric; the loss of water must be compensated by a gain in oil. This is clear from (2.2). Therefore when averaging the time derivative term, we may write more carefully, 1 V  ∂ρs dV ∂t V  and by the previous statement that V is constant in time we can remove the derivative from the integrand,     ∂ 1 ∂t V  ρs dV  . V  We complete the argument by noting that,  ρs =    ρ s  x ∈ Vs   0  x∈ / Vs  .  We may now use the divergence theorem for averages (2.1b) to simplify the  12  2.2. Water Saturation continuity equation further, φ  1 ∂s + ∇ · qs + ∂t V  qs · n dA = 0 Aint  where we assume the density is constant and thus it no longer appears in the equation. We assume that the soil is uniformly homogeneous and thus that the porosity is constant so we have removed porosity from the derivative. We have also used the fact that, φs ≡  Vs . V  Note that here the interfacial area is the total integral from both the water-oil interface (so) and the water-rock interface (sr), Aint = Aso + Asr . We can make the surface integral in the continuity equation vanish by imposing the boundary condition that the normal component of the velocity along the interfacial boundaries vanishes, nij · qsi = 0,  on Aij .  This is valid since water is localized entirely in Vs and thus it does not have an outward flux into the other regions. By using this we then get for our final mass conservation differential equation, φ  ∂s + ∇ · qs = 0 ∂t  (2.3)  13  2.2. Water Saturation where we have removed the average brackets for simplicity and assume hereafter that all quantities are referring to volume averaged ones, as is standard in porous media analysis. We can generate a similar expression for the oil and if we combine that expression with equation (2.2), we see that, ∇ · (qs + qo ) = 0  (2.4)  which translates to conservation of volume, the natural assumption we made above. Now, for our momentum equations we take the standard porous media flow version of the Stokes equation which is Darcy’s Law [21], κs ∇Ps , µs κo qo = − ∇Po µo  qs = −  (2.5) (2.6)  where q is the water flux like before, κ is the permeability, µ is the viscosity and P is the pressure. We note that the subscript denotes association with the water or oil phase. There is a subscript on the pressure because of the immiscibility of the fluids. The pressure difference across the interface of the liquids will create a contact surface. We denote this possible pressure difference as the capillary pressure Pc , P c = Po − Ps .  (2.7)  We will also allow the permeability, κ to be a function of the saturation. This is due to the fact species have a greater mobility when occupying larger volumes. 14  2.2. Water Saturation  2.2.1  Capillary Pressure  In order to derive an expression for the capillary pressure, we will consider a standard Euler-Lagrange method for deriving the surface tension between two immiscible substances in contact. For derivation purposes we will consider a 1-dimensional interface (a curve) but the analysis that follows is readily carried to higher dimensions. Consider a droplet of oil with height h(x) and arc length Γ on a substrate with a surrounding media of water such that the length of contact between the oil and the substrate is L as depicted in Figure 2.2. water  oil  x=0  h(x)  Γ  x=L  Figure 2.2: A 1D model interface contact, Γ, between oil and water. There are three contact surfaces of interest in Figure 2.2. First there is contact between water and the substrate, and we will let γss be the energy per unit length of this surface. Secondly, there is the contact surface between oil and the substrate and γos will denote the energy per unit length of that surface. Finally there is the droplet interface between the oil and water. Consider the energy per unit arc length of the droplet surface to be γ such  15  2.2. Water Saturation that the energy of the entire surface is2 , | dΓ| + γos L + γss Lss  E=γ  where Lss is the length of the substrate that makes contact with the water. Now the forces acting will be such that this energy is a minimum. If we parameterize the droplet curve we see that, 1 + [h (x)]2 dx  | dΓ| =  but in order to obtain reasonable functions h(x), we need a constraint. We assume that the area (or volume in higher dimensions) of the oil is fixed, i.e. regardless of h(x) we have that, L  h(x) dx = A, 0  a constant. With this in mind the problem of minimizing energy reduces to finding the h(x) that minimizes the surface energy of the droplet, L  1 + [h (x)]2 dx  Edroplet = γ 0  subject to the constraint. Let u = h(x) and v = h (x) and use the standard Euler-Lagrange equations to minimize this integral, i.e. find a solution to, d Fv = Fu dx 2  We assume that this derivation occurs in the absence of a gravitational field and hence there is no need to overcome deformation due to gravity.  16  2.2. Water Saturation where the subscript indicates differentiation. Here, F is the functional, √ F = γ 1 + v 2 + λu with λ a Lagrange multiplier. Upon carrying out these operations we get that, γ  d v √ = λ. dx 1 + v 2  Now, we notice that the curve y = h(x) = u has normal vector, 1 ˆ=√ v, −1 n 1 + v2 and so we can write our Euler-Lagrange expression as ˆ = λ, γ∇ · n γν = λ where we define the curvature ν as the standard divergence of the normal vector, ˆ. ν =∇·n This derivation can easily be extended to higher dimensions. In our case we have a 2-dimensional surface energy to minimize and the parameter γ is known as the surface tension coefficient. It is known ([17],[25]) that the Lagrange multiplier represents the pressure difference across the interface, or Pc in our case. This dimensionally agrees with the constraint (i.e. the Lagrange multiplier has the appropriate units for pressure) and also produces  17  2.2. Water Saturation the Young-Laplace equation, γν = Pc . Now Leverett in [25] determined an expression for the curvature in terms of the saturation s, ν=  φ J(s) κ  where J(s) has since become to be known as the Leverett function, which is a dimensionless function that includes saturation dependence. The lack of subscript on the permeability κ is due to the fact that this is the total permeability. Leverett’s motivation for creating such an expression was simply to derive a non-dimensional form of a curvature that could be used in his experiments to determine the dependency of capillary pressure solely on saturation. He derived it by assuming that the contact surface area must be a fraction of the surface of the sand itself or by noting that it must depend on the voids in the sand which are determined through the porosity and permeability. His experiments determined that the curvature varied quite substantially over different ranges of saturation. We will treat J(s) as any monotonically decreasing function that incorporates the coupling of water saturation to the pressure difference whether it truly conforms with Leverett’s experiments or not. The reason for the monotonicity restriction is so that the pressure difference drives the fluid left to right. With this Leverett formulation we have an expression for the capillary pressure,  Pc = γ  φ J(s). κ  (2.8)  18  2.2. Water Saturation  2.2.2  Non-Dimensionalization  Consider the rescaling, x¯ =  x , L  t t¯ = , τ  ¯i = q  qi , Q  αi =  κi , κ  P¯i = Pi  κ µs QL  where L is a characteristic length, τ is a characteristic time, κ is the total permeability, αi is the relative permeability of species i, and Q is a characteristic velocity. Using these non-dimensional variables, we get for our equations (after removing the overbars), ∇ · (qs + qo ) = 0,  (2.9a)  ∂s + ∇ · qs = 0, ∂t  (2.9b)  φFl  qs = −αs (s)∇Ps ,  (2.9c)  qo = −δαo (s)∇Po ,  (2.9d)  Po − Ps = Ca−1 J(s),  (2.9e)  where here we introduce the non-dimensional numbers, µs , Viscosity ratio, µo µs QL Ca = √ , Capillary Number, γ κφ L Fl = , Flow Number. Qτ δ=  (2.10a) (2.10b) (2.10c)  19  2.2. Water Saturation Now, we wish for all of the terms in (2.9b) to be present with equal magnitude so we will choose our flow rate, Q, such that φFl = 1, i.e., Q=  Lφ . τ  (2.11)  We choose as our two free parameters a typical site length scale of 20m (400 square-metre site) and, as is consistent in literature ([21],[36],[13],[8]), water flow rates are based on a per day scale so we choose a time scale of 1 day. Table 2.1 defines all of the parameter values used as well as references where required.  Symbol φ L τ Q κ µs µo γ δ Ca  Table 2.1: Parameter Values for Model. Value Description Reference 0.49 Porosity [21] (average value) Characteristic 20m Chosen length scale Characteristic 1 day Chosen time scale 9.8m/day Characteristic velocity Eqn (2.11) 7.1 × 10−12 m2 Total permeability [21] (average value) 1.00 × 10−3 Pa·s Water viscosity (300K) [23] 0.523Pa·s Oil viscosity(300K) Chosen Interfacial tension 32.7mN/m [22] of water & benzene 1/523 Viscosity Ratio Eqn (2.10a) 37.2 Capillary number Eqn (2.10b)  In [36], Ritter states that ambient groundwater flow rates are about 15m/day which is the same order of magnitude as predicted by (2.11). The values in Table 2.1 are consistent and necessary for porous media considerations and the use of Darcy’s Law, i.e. the appropriate length and time scales are present 20  2.2. Water Saturation [21]. For the viscosity of oil, we choose our value based on the knowledge that we are aiming to extract crude oil. Motor oil has a viscosity of 0.319Pa·s ([11]) and crude oil has a higher viscosity than this ([10]) hence we choose a slightly higher value. However due to the complex structure of crude oil it is difficult to determine other parameters for it and so we choose benzene to represent our oil for the other parameters. This is due to the fact that benzene is a common component in crude oil and involved in most industrial processes that cause oil seepage into the earth [45]. Now, in (2.9e) it is not the capillary number but rather the inverse Capillary number that is present. By looking at Table 2.1 we see the inverse Capillary number is small and we will let, = Ca−1 = 0.03. We also note that by our choice of parameters, δ  1 as well, though in  general it does not have to be.  2.2.3  Cory-Type Permeabilities  We now have a fairly tractable model but we see the presence of unknown functions in the form of the relative permeabilities, αi (s). We follow, [5],[34], and choose a polynomial fit with respect to s which is known as a Cory-Type permeability. Following [5] we choose a cubic exponent to get that, αs (s) = s3 ,  (2.12a)  αo (s) = (1 − s)3 ,  (2.12b)  21  2.3. Temperature which is always positive since 0 ≤ s ≤ 1. It is worth noting here that physically, there are upper and lower bounds to water saturation because in standard conditions, soil is rarely fully saturated or depleted of water. However, we can always choose a suitable rescaling of s such that we map those bounds to be s = 0 and s = 1. With this choice of permeabilities we do not exhibit any percolation thresholds [46].  2.3  Temperature  We also wish to include the effects of temperature into the soil model for two main reasons. First, many modern in-situ methods involve a mechanism for heating the water, oil, or both which indicates that its inclusion should significantly effect the model [27] and second, in real applications, viscosities can be strongly temperature dependent as temperature affects the kinetic energy of the molecules, lowering their resistance to motion, which viscosity is a measure of. We will first introduce a macroscopic energy equation for temperature effects via the volume averaging technique previously used and then identify a temperature dependent model for the viscosities to completely couple the problem. For completeness, we will introduce how the volume averaging method can lead to the inclusion of anisotropic effects such as structural deformities or tortuosity in the porous structure but ultimately our model will ignore such effects.  22  2.3. Temperature  2.3.1  Volume Averaged Temperature Equation  To begin the analysis we consider a standard “point” heat equation which is a continuous model valid in each subdomain volume. We have, ∂Tr = −∇ · JT r ∂t ∂Ts (ρcp )s = −∇ · JT s ∂t ∂To = −∇ · JT o (ρcp )o ∂t (ρcp )r  for rock, water, and oil respectively. We denote JT for thermal fluxes in the problem, cp as specific heat capacities and ρ as the density. We assume that ρcp is constant, hence its exclusion from the time derivative. We take the volume average of these equations to get, ∂Tr ∂t ∂Ts (ρcp )s ∂t ∂To (ρcp )o ∂t (ρcp )r  = − ∇ · JT r  (2.13a)  = − ∇ · JT s  (2.13b)  = − ∇ · JT o .  (2.13c)  23  2.3. Temperature We will deal with the various fluxes individually. First consider the left-hand side of (2.13a), (ρcp )r  ∂Tr ∂t   = (ρcp )r    Vr ∂  1 V ∂t Vr  Tr dV  Vr  = (1 − φ)(ρcp )r  ∂ Tr ∂t  ∗  where we have exploited the fact that Vr does not change in time and have rewritten in terms of intrinsic volume averaging. We can do a similar simplification for (2.13b) with slightly more caution because here the volume of water, Vs = Vs (t). However, due to the fact that we disallow any net changes in volume we can exploit the fact that the control volume V is time independent as was the case with the derivation for water saturation, ∂Ts ∂t  1 ∂ = (ρcp )s  ∂t V (ρcp )s   Ts dV  Vs  = φ(ρcp )s  ∂ (s Ts ∗ ) ∂t  24  2.3. Temperature where we have removed φ from the derivative under the assumption the porosity is constant. We can use the same methodology for (2.13c), ∂To ∂t ∂ = φ(ρcp )o ((1 − s) To ∗ ) . ∂t (ρcp )o  We will now analyze the various fluxes that contribute to thermal changes in the system and derive their local averages. Convective Fluxes We allow for the transport of water or oil throughout the soil and as such, there will be an associated heat flux with that motion. The convective heat fluxes are given by, Jconv(s) = (ρcp )s qs Ts Jconv(o) = (ρcp )o qo To where qi are the velocities of each fluid. We note from (2.4) that, ∇ · (qs + qo ) = 0 which is the same conservation of volume required to state that the control volume V was time independent. We do not have a convective flux for the solid matrix (rock) because it is stationary. We will do the analysis for the water and that will trivially carry over to the oil convective flux. We can  25  2.3. Temperature write the volume average of the divergence of the convective flux as, ∇ · ((ρcp )s qs Ts ) = (ρcp )s ∇ · qs Ts + (ρcp )s  1 V  qs Ts · n dA  (2.14)  Aint  using identity (2.1b). Now before proceeding we need some extra definitions. Firstly, it is unlikely that information about the true microscopic variables will be known yet they appear in the surface integral of (2.14). To avoid this, it is customary to write for any quantity w, wi = wi  ∗  +w ˜i  which says that the true microscopic quantity is the intrinsic volume average of that quantity plus any spatial variations from the average denoted w˜i . Notice if we take the intrinsic volume average of that definition, wi  ∗  = wi  w˜i  ∗  = 0.  ∗  + w˜i  where we have used the fact that over Vi , wi  ∗  ∗  is constant. Using these two  facts we can write, qs Ts = qs  ∗  Ts  ∗  ˜s + q ˜ s T˜s . + qs ∗ T˜s + Ts ∗ q  If we then take the intrinsic volume average of this, use the identities for w˜i and the fact that a volume averaged quantity is constant over its intrinsic 26  2.3. Temperature volume we get, qs Ts  ∗  = qs  ∗  qs Ts = φs  Ts qs  ∗  ∗  ˜ s T˜s + q Ts  ∗  ∗  ˜ s T˜s + q  ∗  .  We may substitute this expression into (2.14) to get Jconv(s) = ∇ · ((ρcp )s qs Ts ) = (ρcp )s ∇ · (φs qs + (ρcp )s  1 V  ∗  ˜ s T˜s Ts ∗ ) + (ρcp )s ∇ · φs q  ∗  qs Ts · n dA. Aint  We can define a similar convective flux for oil, Jconv(o) = ∇ · ((ρcp )o qo To ) = (ρcp )o ∇ · (φ(1 − s) qo + (ρcp )o  1 V  ∗  ˜ o T˜o To ∗ ) + (ρcp )o ∇ · φ(1 − s) q  ∗  qo To · n dA. Aint  We now need appropriate boundary conditions for the problem. We assume that at any interface, Aij , Ti = Tj enforcing continuity of temperature. We will also force the normal velocity across any interface to be zero to disallow for outgoing flux of any species,  27  2.3. Temperature consistent with our derivations of water flux, nij · qi = 0, on Aij . Now we also want to have a locally average thermal equilibrium so that we only need to consider one temperature. This is a fair assumption to make since we have no heat sources inside any of the subdomain volumes ([21]). This means that, Ts  ∗  = To  ∗  = T .  In order to produce closure arguments, it is typical ([21],[30]) to rewrite the quantity T˜i as, T˜i = bi · ∇ T with bi = bi (x) a vector in order to allow for anisotropy. This allows us to write the variations, T˜i as a projection of the average temperature gradient. The bi quantities effectively represent the inherent property of the material being studied such as structure and chemical composition. This is explained in more detail in [21] as well as methods for approximating their values mathematically or obtaining them empirically. With this representation, we can write, 1 Vi  1 ˜ i T˜i dV = q Vi Vi  ˜ i bi · ∇ T dV = q ˜ i bi q  ∗  ·∇ T  Vi  28  2.3. Temperature where the product of two vectors u and v is to be understood as a dyadic product,     uv uv uv  1 1 1 2 1 3   uv = u2 v1 u2 v2 u2 v3  .   u3 v1 u3 v2 u3 v3 Here, the subscript indicates the position within the vector. Notice that we have removed the gradient of the average temperature from the integral as this is constant over its subdomain volume. We denote the tensor as, ˜ s bs D = (ρcp )s φs q  ∗  ˜ o bo + (ρcp )o φ(1 − s) q  ∗  and understand it as a dispersion tensor since it accounts for the temperature deviations due to a spatially non-uniform velocity. Summing the convective flux terms and using the appropriate boundary conditions noting that nij = −nji , Jconv = ∇ · ([(ρcp )s φs qs  ∗  + (ρcp )o φ(1 − s) qo ∗ ] T ) + ∇ · (D · ∇ T ) . (2.15)  29  2.3. Temperature Diffusive Fluxes Even in the absence of motion, there will be a natural conduction of heat due to molecular contact, which we assume takes the form of Fourier’s Law, Jcond(s) = −ks ∇Ts , Jcond(o) = −ko ∇To , Jcond(r) = −kr ∇Tr , where we notice here that the stationary solid matrix is able to absorb or dissipate heat through conduction. Here we define ki as the thermal conductivity of species i and treat it to be constant. As with the convective fluxes we will do the full analysis for water and then generalize to the other two fluxes. We can write the volume average of the divergence of the conductive flux as, ∇ · (ks ∇Ts ) = ∇ · ks ∇Ts +  ks V  ∇Ts · n dA  (2.16)  Aint  30  2.3. Temperature using identity (2.1b). Since we assume that ks is constant we can use identity (2.1a) to write, ks ∇Ts = ks ∇ Ts + ks  1 V  Ts n dA Aint  = ks ∇(φs Ts ∗ ) + ks  1 V  Ts n dA Aint  = ks φs∇ Ts  ∗  + ks Ts ∗ ∇(φs) + ks  1 V  Ts n dA  (2.17)  Aint  where we have rewritten the volume averages in terms of intrinsic ones. Next we once again define the true temperature in terms of the average and deviation via, Ts = Ts  ∗  + T˜s  and do a surface integral over the interfacial area, 1 V  Ts n dA = Aint  1 V  Ts ∗ n dA + Aint  1 V  T˜s n dA.  (2.18)  Aint  Now for the second term in this expression we can reverse identity (2.1a) to get that, 1 V  Ts ∗ n dA = ∇ Ts  ∗  − ∇ Ts  ∗  .  Aint  We can write, ∇ Ts  ∗  = ∇ (φs Ts ∗ ) = Ts ∗ ∇(φs) + φs∇ Ts  ∗  31  2.3. Temperature by utilizing the fact that the intrinsic volume average is constant over its subdomain volume. Now we may write, ∇ Ts  ∗  = φs∇ Ts  ∗  by utilizing that the gradient of the average is constant over its subdomain volume. Finally we can use these simplifications in (2.18) to get that, 1 V  Ts n dA = Aint  1 V  T˜s n dA − Ts ∗ ∇(φs). Aint  This integral simplification appears in [7] (for the solid matrix) using order of magnitude arguments. Now if we substitute this into (2.17) and simplify we get, ks ∇Ts = ks φs∇ Ts  ∗  +  ks V  T˜s n dA Aint  so that the conductive flux for water is,  Jcond(s) = −∇ · (ks φs∇ Ts ∗ ) − ∇ ·    ks T˜s n dA − V  ks V Aint  ∇Ts · n dA. Aint  32  2.3. Temperature In a similar fashion we can derive conductive fluxes for oil and the rock,  Jcond(o) = −∇ · (ko φ(1 − s)∇ To ∗ ) − ∇ ·    ko V  T˜o n dA Aint  −  ko V  ∇To · n dA, Aint     Jcond(r) = −∇ · (kr (1 − φ)∇ Tr ∗ ) − ∇ ·   kr V  T˜r n dA Aint  kr − V  ∇Tr · n dA. Aint  Now as for the convective fluxes we will assume that we have thermal equilibrium so that all the intrinsic volume average temperatures are equal to the global volume average temperature. For boundary conditions we will impose that at each interfacial surface we will have continuity of temperature and temperature flux, i.e. on surface Aij , Ti = Tj nij · (ki ∇Ti ) = nij · (kj ∇Tj ). Notice that the first boundary condition along with thermal equilibrium implies that T˜i = T˜j . Now like with convection, we will rewrite the temperature variations by, T˜i = bi · ∇ T  33  2.3. Temperature so that, T˜i n dA = Aint  bi n dA∇ T Aint  where once again the product of vectors is understood to be a dyad product. By summing the three fluxes, using the thermal equilibrium and boundary conditions we get that the total conductive flux can be expressed as, Jcond = −∇ · ((ks φs + ko φ(1 − s) + kr (1 − φ))∇ T + K · ∇ T )  (2.19)  where, K=  ks − ko V  bs nso dA + Aso  ks − kr V  bs nsr dA + Asr  ko − kr V  bo nor dA Aor  is the conduction tensor. Notice we could incorporate all of the terms multiplying the gradient into K but we separate them to visualize the presence of the standard Fourier conductive terms. It is clear now as to why the reformulation using ∇ T was used for the variations in the convective fluxes that lead to the dispersion tensor D as it can now be incorporated with K as a total tensor that multiplies the temperature gradient, measuring the property effects such as structure and tortuosity. Fully Averaged Temperature Equation We know have averaged expressions for both the conductive (2.19) and convective (2.15) fluxes so that we can finally substitute all the information into  34  2.3. Temperature (2.13) and write the averaged heat equation as, φ(ρcp )s  ∂ ∂ T ∂ (s T ) + φ(ρcp )o ((1 − s) T ) + (1 − φ)(ρcp )r ∂t ∂t ∂t  + ∇ · ([(ρcp )s φs qs  ∗  + (ρcp )o φ(1 − s) qo ∗ ] T )  = ∇ · ((ks φs + ko φ(1 − s) + kr (1 − φ))∇ T ) + ∇ · ((K − D) · ∇ T ) . Now if we expand the time derivative out and use (2.3) for the following, ∂s = −(ρcp )s ∇ · qs ∂t ∂s φ(ρcp )o = (ρcp )o ∇ · qo ∂t φ(ρcp )s  and take advantage of the fact that, φs qs  ∗  = qs  φ(1 − s) qo  ∗  = qo  then we can drastically simplify the problem, (φ(ρcp )s s + φ(ρcp )o (1 − s) + (1 − φ)(ρcp )r )  ∂T ∂t  + (φ(ρcp )s qs + φ(ρcp )o qo ) · ∇T = ∇ · ((ks φs + ko φ(1 − s) + kr (1 − φ))∇T ) + ∇ · ((K − D) · ∇T ) . (2.20) Notice we have dropped the average brackets for convenience and understand that the relevant variables are the same as introduced in the water saturation derivation, i.e. the velocities qi are understood to be volume averages and  35  2.3. Temperature not intrinsic volume averages. As mentioned we included the derivation encompassing microscopic effects for completeness only and so in what follows we will assume that there is no spatial variations from any average values and thus set the tensors K and D to be zero. Of course, this in general, is inappropriate to ignore but we wish to simplify the analysis.  2.3.2  Temperature Dependent Viscosity  We currently have an expression for the water saturation (2.3) and an expression for the temperature (2.20) but as it stands the systems are not coupled because while s appears in the temperature equation, the temperature, T , is absent from the water saturation. The most natural way to couple the system is to impose that the viscosity is a function of temperature, µi = µi (T ). This will make the velocity and thus the saturation depend on T . We will assume that water remains a constant viscosity even with the introduction of temperature on the basis that its changes will be minimal in relation to that of oil. For oil we will assume an exponential type behaviour whereby the viscosity takes on a known value at a reference temperature Ta and then decays exponentially to be of equal viscosity as water. Therefore we assume, µs (T ) = µs0 µo (T ) = (µo0 − µs0 ) exp −  (T − Ta ) Ta  + µs0  36  2.3. Temperature where µi0 is the reference viscosity at temperature Ta . These forms for the viscosities have not been argued valid on a physical basis, they are merely introduced as a simplified way of incorporating the thermo-dependence to investigate the solution structure in this case. However, the exponential dependence does match significantly well with experimental reports of viscosities [2].  2.3.3  Non-Dimensionalization  First we consider the same rescaling as before, x¯ =  x , L  t t¯ = , τ  ¯i = q  qi , Q  αi =  κi κ  where all parameters are as defined in Table 2.1. We will rescale temperature by an ambient value, Ta , T T¯ = Ta consistent with the introduction of a temperature dependent viscosity. We will also scale the viscosities by their respective reference values, µ ¯i =  µi . µi0  With the use of this rescaling we may utilize the same scaling for pressure as in the temperature independent model, P¯i = Pi  κ µs0 QL  37  2.3. Temperature where we now consider the reference values since viscosity is no longer constant. Our non-dimensional equations can then be written as (removing overbars), ∇ · (qs + qo ) = 0,  (2.21a)  ∂s + ∇ · qs = 0, ∂t  (2.21b)  µs (T ) = 1,  (2.21c)  µo (T ) = (1 − δ) exp (1 − T ) + δ,  (2.21d)  φFl  qs = −αs (s)∇Ps , qo = −δ  (2.21e)  αo (s) ∇Po , µo (T )  Po − Ps = Ca−1 J(s), ∂T ˆ + (Cs Cˆs (s)qs + Co Cˆo (s)qo ) · ∇T = D∇ · (D(s)∇T ) ∂t  (2.21f) (2.21g) (2.21h)  where the nondimensional numbers δ, Ca and Fl are defined as before but with reference viscosity values used. We define the convection numbers Ci as the average ratio of convective energy transport to total heat energy, Cs = Co =  φ(ρcp )s Fl  1 0  (φ(ρcp )s s + φ(ρcp )o (1 − s) + (1 − φ)(ρcp )r ) ds φ(ρcp )o  Fl  1 0  (φ(ρcp )s s + φ(ρcp )o (1 − s) + (1 − φ)(ρcp )r ) ds  (2.22) (2.23)  38  2.3. Temperature and the diffusion number D as the average ratio of diffusive energy transport to total heat energy,  D=  1 0  QLFl  1 0  (ks φs + ko φ(1 − s) + kr (1 − φ)) ds  (φ(ρcp )s s + φ(ρcp )o (1 − s) + (1 − φ)(ρcp )r ) ds  .  (2.24)  Notice that the diffusion number acts like an averaged inverse Peclet number. ˆ are linear order one functions that contain the s The functions Cˆi and D dependence. They are order one because 0 ≤ s ≤ 1. These functions are what remain after the extraction of the convection or diffusion numbers. We do a rescaling based on those numbers so that we can observe the average effect each process has on the energy equation. New parameter values are listed in Table 2.2 which is meant to complement the previous list of parameters in Table 2.1.  39  2.3. Temperature  Table 2.2: Parameter Values for Model with Temperature. Value Reference Symbol Description s o r s o r Heat Capacity cp 4.18×103 1.76×103 794 [15] [28] [4] (Jkg−1 K−1 ) Density ρ 1000 877 2.65×103 [14] [26] [3] (kgm−3 ) Thermal k 0.610 0.151 1.04 Conductivity [33] [26] [44] (Wm−1 K−1 ) Reference Ta 300K Chosen Temperature µs0 Reference −3 1.00×10 Pa·s Water [23] Viscosity µo0 Reference 0.523Pa·s Oil Chosen Viscosity Cs Water 0.405 Convection Eqn (2.22) Number Co Oil 0.150 Convection Eqn (2.23) Number D Diffusion 7.24×10−10 Eqn (2.24) Number  We have chosen reference values for temperature and viscosity such that when T = 1 in a non-dimensional sense the non-isothermic problem resolves to the model absent of temperature. As stated before, all of the parameters (except for viscosity) match the literature values for benzene. The density of oil was chosen to match the corresponding thermal conductivity. For soil values we have chosen to use coarse quartz sand, a common element in soils. Based on  40  2.3. Temperature the convection and diffusion numbers, the thermal transport is dominated by convective flow.  Now that we have a full model we will begin to analyze it from a basic one-dimensional isothermic flow in the absence of capillary pressure. From there we will then incorporate capillary pressure before finally incorporating pressure effects in the non-capillary setting.  41  Chapter 3 Results 3.1  1D Model in the Absence of Capillary Pressure  In order to begin analyzing the interaction at the porous level between oil and water, we consider the isothermal model (2.9) in a 1-D setting in the absence of capillary pressure. This is a reasonable first approximation since 1. In this case, the water and oil pressures are the same and we denote them as P to get a 1D model, ∂ (qs + qo ) = 0, ∂x ∂s ∂qs + = 0, ∂t ∂x dP , dx dP qo = −δαo (s) . dx qs = −αs (s)  (3.1a) (3.1b) (3.1c) (3.1d)  The benefit of the 1D model is that we see the divergence condition becomes a simple 1-dimensional conservation law (3.1a). We can then quickly solve  42  3.1. 1D Model in the Absence of Capillary Pressure that to get, qs + qo = 1  (3.2)  where we see that the non-dimensional Q from Table 2.1 represents the total flow in a dimensional sense. By summing (3.1c) and (3.1d) while using (3.2) we can simplify the model to get, ∂ ∂s + f (s) = 0, ∂t ∂x  − ∞ < x < ∞, t > 0  (3.3)  where, for the relative permeability functions assumed here, f (s) =  αs s3 = 3 3 s + δ(1 − s) λ  (3.4)  with, λ = αs + δαo = s3 + δ(1 − s)3 being denoted as the total permeability. Note that this is the famous BuckleyLeverett formulation as postulated in [6]. We will choose as initial conditions a standard Riemann problem with a planar front separating water and oil, i.e.,  s(x, 0) =    1 x < 0  .   0 x > 0  43  3.1. 1D Model in the Absence of Capillary Pressure We also implicitly define that the Riemann problem holds in the far-field for all time, i.e., lim s(x, t) = 0  x→∞  lim s(x, t) = 1.  x→−∞  Now, from standard PDE theory ([12],[31]) we know that in a hyperbolic problem, if the flux function (f (s)) has a strictly positive second derivative then there is a shock; conversely, if f (s) is strictly negative for all s then there is a rarefaction. Here, f (s) as written, changes convexity and so we anticipate both types of solution. By the method of characteristics, the equation (3.3) reduces to, ds =0 dt along characteristics that satisfy, dx = f (s). dt Now since along these characteristics, s does not depend on t, we can easily integrate and use the initial conditions to obtain that the characteristics are, x=β  44  3.1. 1D Model in the Absence of Capillary Pressure where x(0) = β and the solution is then,  s(β) =    1 β < 0  .   0 β > 0 All the characteristics that satisfy β < 0 have s = 1 while β > 0 have s = 0 and thus we have a jump in the solution. Therefore, we fit a rarefaction fan that smoothly connects the two solutions. The problem point is β = 0 and so, knowing that the characteristic lines that form the rarefaction must emerge from the origin, we seek a solution, x = s(χ) t  s=s  there. Substituting this into (3.3) we get that, ds (f (s) − χ) = 0. dχ Since we expect that s = s(χ) we assume that the derivative is non-zero and thus we get for our rarefaction fan, f (s) =  x , t  0 ≤ s ≤ 1.  A full picture of the characteristic diagram with the expansion fan is in Figure 3.1.  45  3.1. 1D Model in the Absence of Capillary Pressure s=1  x=0  x=ξ  s=0  t=0  Figure 3.1: Characteristic diagram with expansion fan and shock for the 1D problem (3.3) in the absence of capillary pressure. Now as the expansion fan meets with the s = 0 solution, a shock x = ξ(t) forms (Figure 3.1). We know that the conditions across the shock must satisfy a Rankine-Hugoniot condition ([31]), v=  f (s+ ) − f (s− ) dξ = dt s+ − s−  where v is the shock velocity, s− is the saturation value just before the shock and s+ is the saturation value just after the shock. Since the shock begins to form at the first intersection of characteristics (the origin), the shock path will be a straight line from the origin which means that the shock is also a characteristic of the rarefaction wave. Therefore, v has the form, v = f (s). The shock profile must be linear because there will always be a uniform spacing between characteristic intersections due to the fact that there are an 46  3.1. 1D Model in the Absence of Capillary Pressure infinite number of vertical characteristics with which to annihilate the infinite number of rarefaction characteristics. We know that the solution after the shock is s = s+ = 0 but we do not know the exact value of s for which the rarefaction characteristic becomes a shock. To find this value we can take advantage of the form of v to find s− = s∗ < 1 that satisfies, f (s∗ ) =  f (s∗ ) . s∗  Now using the expression for f (s) we get that a non-trivial s∗ must satisfy (δ − 1)s∗3 − 3δs∗ + 2δ = 0. In the case where δ  (3.5)  1 then (δ − 1) ≈ − 1 and so, s∗3 + 3δs∗ − 2δ ≈ 0.  To leading order in δ, this only produces trivial roots but if we scale s∗ = δ 1/3 k ∗ and expand, k ∗ = k0 + δ 1/3 k1 then we get that the asymptotic shock saturation value s∗a is, s∗a ∼ 21/3 δ 1/3 − 2−1/3 δ 2/3 .  (3.6)  For the specific choice of δ from Table 2.1, we get that s∗a ≈ 0.1442 and by numerically solving (3.5) we get the actual shock saturation is s∗ ≈ 0.1443,  47  3.1. 1D Model in the Absence of Capillary Pressure a good agreement. Now for the shock velocity, v=  f (s∗ ) s∗2 = s∗ s∗3 + δ (1 − s∗ )3  (3.7)  which we can approximate using the asymptotic shock saturation value. For simplification, we will take only the first term in s∗a which will provide accuracy up to o(δ 1/3 ). By substituting this into the expression for shock velocity (3.7) and expanding for δ  1 we get that the asymptotic shock velocity va  is, va ∼  22/3 2 + . 1/3 3δ 3  For our choice of parameters, we get that va ≈ 4.93 while numerically we calculate from (3.7) that v ≈ 4.95, another good agreement. By finding s∗ and v, the problem is complete and has solution,     1    s(x, t) = [f ]−1      0  x≤0 x t  0 < x < f (s∗ )t .  (3.8)  x > f (s∗ )t  For the δ we chose, the solution is plotted in Figure 3.2 for t = 1/8 which in dimensional time is 3 hours.  48  3.1. 1D Model in the Absence of Capillary Pressure Saturation Profile at t=0.125 1 0.9 0.8 0.7  s  0.6 0.5 0.4 0.3 0.2 0.1 0 -2  -1  0  1  2  3  4  5  x  Figure 3.2: Saturation profile from solution (3.8) with s∗ = 0.1443 and t = 1/8. Now, notice that if we consider a moving coordinate frame, η = x − vt then (3.3) becomes, ∂ ∂s + (f (s) − vs) = 0 ∂t ∂η and if we seek steady state solutions then we set the time derivative to zero. If we integrate this expression then we see a steady state solution s0 exists if d (f (s0 ) − vs0 ) dη = f (s0 ) − vs0 = C, dη  49  3.1. 1D Model in the Absence of Capillary Pressure a constant. When C = 0 then f (s∗ ) =v s∗ satisfies this condition. These are precisely the conditions that describe the shock and thus the shock front,  s0 =    s ∗ ,  η<0 (3.9a)   s − , η > 0 is a steady-state solution in this frame of reference [47]. What this means is that, the steady state solution is a step profile with s = s∗ to the left and s = 0 to the right of the shock. We can also use (3.2) and (3.1c) to define a steady state pressure as, dP0 1 1 =− 3 3 = − , dη λ s0 + δ (1 − s0 )  η = 0.  (3.9b)  We will now analyze the stability of this travelling wave solution under perturbations in higher dimensions.  3.1.1  Linear Stability Analysis  We consider the steady-state shock solution (3.9a) to be a 1D cross section of a 2-dimensional problem. We wish to see how the shock front behaves in the presence of transverse wave disturbances, i.e. we form a linear expansion  50  3.1. 1D Model in the Absence of Capillary Pressure as3 , s ∼ s0 + εs1 (η) exp (iny + σt) , P ∼ P0 + εP1 (η) exp (iny + σt) , s1 ,  ds1 dP1 , P1 , → 0, dη dη  |η| → ∞  η=0 η=0  (3.10a) (3.10b) (3.10c)  where s0 and P0 are the steady state saturations and pressures respectively given by (3.9) and n > 0 are the transverse wave modes. The sign of the real part of σ then determines the stability of the shock to disturbances of wave number n. The parameter ε is an artificial parameter introduced just to control the balance of terms, it is not to be confused with , the inverse capillary number. The form of the perturbation is also mathematically valid since the full two-dimensional PDE decouples the variables x, y, and t from one another (i.e. produces constant coefficient problems). Hence the expansion is simply a Fourier transform eigenfunction expansion in the y direction followed with an exponential ansatz (or Laplace Transform) for the time dependence. We choose our far-field conditions on s1 and P1 such that saturation and pressure decay smoothly to the steady state solution. If we substitute the ansatz (3.10) into the full isothermal PDE, (2.9) (while still 3  It may seem peculiar that velocity is not being perturbed but indeed it is due to the perturbation of pressure and the explicit form for velocity in terms of pressure as Darcy’s Law.  51  3.1. 1D Model in the Absence of Capillary Pressure excluding capillary pressure)4 , ∂s + ∇ · (qs − vsˆ x) = 0 ∂t qs = −αs ∇P qo = −δαo ∇P ∇ · (qs + qo ) = 0 ˆ is a unit vector in the η direction, then, keeping only linear terms where x (i.e. terms of O(ε)), we get for our eigenvalue problem, dP0 dP1 d + α¯s s1 + vs1 − n2 α¯s P1 = σs1 α¯s dη dη dη d ¯ dP1 ¯ dP0 ¯ 1 = 0. λ +λ s1 − n2 λP dη dη dη  (3.11a) (3.11b)  Here, we let over-bars indicate a function evaluated at the steady state, and prime indicates differentiation with respect to s. Now since we have localized around the shock front, the steady-state away from η = 0 is simply constant ¯ and its derivatives are all constants and can be removed from (3.9) so, α¯s , λ, the derivatives. We can rearrange (3.11) solely for s1 , 1− 4  f (s0 ) v  ds1 σ = s1 dη v  We have redefined the gradient operator as, ∇=  ∂ ∂ , ∂η ∂y  .  52  3.1. 1D Model in the Absence of Capillary Pressure by noting that the identity, f (s0 ) = vs0 holds for both steady state saturation values. Now if η > 0 then s0 = 0 and we get that, s1 = D+ exp  σ η v  but for any unstable modes5 (σ > 0) this cannot vanish as η → ∞ and so we conclude that D+ = 0. If η < 0 then s0 = s∗ < 1 and f (s∗ ) = v so therefore, 0=  σ s1 . v  Thus for general σ, s1 = 0. Therefore, s1 = 0,  η = 0.  (3.12a)  Now, knowing that s1 = 0 away from η = 0 means that we get a simplified expression for P1 ,  d 2 P1 = n 2 P1 2 dη  5  Fingering is a phenomena due to instability and thus we are not interested in any modes that stabilize.  53  3.1. 1D Model in the Absence of Capillary Pressure which has solution,  P1 =    C+ exp (−nη) , η > 0  C− exp (nη) ,  (3.12b)  η<0  in order to satisfy the far-field conditions in both cases. This analysis provides us with the solution that s1 is zero everywhere, except possibly at η = 0. In [47] Yortsos and Huang analyze a similar problem but for a true front (s = 1 transitions to s = 0) and claim that an appropriate choice of s1 is s1 = Aδ(η)  (3.13)  in order to satisfy the appropriate conditions and conclusion about s1 away from η = 0. Upon doing so, they generate a series of jump conditions across the shock front that are solvability conditions for the eigenvalue, ¯ dP1 = 0, λ dη dP1 = s∗ Aσ, α¯s dη dP0 [P1 ] = −A , dη  (3.14a) (3.14b) (3.14c) (3.14d)  which we have generalized to suit our problem. Note that we let [·] denote the jump from the left state to the right state. While the delta function is equivalent to imposing continuity of physical variables across the shock interface, we instead take the more intuitive approach of directly matching 54  3.1. 1D Model in the Absence of Capillary Pressure these continuities. To do this, we instead treat the problem expansion as before by expanding away from η = 0 but also by perturbing the interface to get, s ∼ s0 + εs1 exp (iny + σt) P ∼ P0 + εP1 exp (iny + σt) ηˆ ∼ 0 + εA˜ exp (iny + σt) where once again ε is a small parameter introduced to control term balance, then we will still conclude that s1 = 0 and that P1 is exponential away from η = ηˆ ≈ 0 and thus the eigenfunctions are still given by (3.12). However at the boundary, ηˆ, we impose proper physical conditions in that the total pressure and normal fluid flux must be continuous. In the travelling wave coordinate frame, the fluid fluxes are given by, qs − vsˆ x qo − v(1 − s)ˆ x for water and oil respectively. We find the normal fluxes by taking an inˆ . However at the ner product with the normal direction of the interface, n interface η = ηˆ an we have that the normal vector is, ˆ= n  1, εAiny exp (iny + σt) 1 + (εAiny exp (iny + σt))2  55  3.1. 1D Model in the Absence of Capillary Pressure which we can expand as, ˆ≈x ˆ + εAiny exp (iny + σt) y ˆ + O(ε2 ) n ˆ is a unit vector in the y-direction. However, to first order, the flux where y terms only have components in the x-direction, qi = qi0 , 0 + ε qi1x , qi1y ˆ =x ˆ. then the only contribution for normal fluxes up to O(ε) is to have n ˆ direction, the total jump in Therefore by taking the inner product in the x water across the interface ηˆ = 0 is, [ˆ x · qs − vs] = total flux across η = ηˆ+ − total flux across η = ηˆ− ∼ total flux across η = 0+ − total flux across η = 0− = s(0+ )  d ηˆ dt  − s(0− ) η=0+  d ηˆ dt  η=0−  with, d ηˆ = εσA exp (iny + σt) dt being the speed of the interface. However we also know that, [ˆ x · qs − vs] = −α¯s  dP0 dP1 − v [s0 ] + ε −α¯s exp (iny + σt) + O(ε2 ) dη dη  where we have explicitly used the fact that, s1 =  ds1 d2 P0 = = 0. dη dη 2 56  3.1. 1D Model in the Absence of Capillary Pressure These two expressions for the water flux must be equal and so we get (after substituting in known values and expressions), [ˆ x · qs − vs] = −α¯s  dP1 dP0 + ε −α¯s exp (iny + σt) dη dη  = −vs∗ − εs∗ σ A˜ exp (iny + σt) which states that dP0 = vs∗ dη dP1 ˜ = s∗ σ A. α¯s dη α¯s  The first of these expressions is satisfied by applying the known expression for the steady-state pressure and the second equation provides the first solvability condition for σ. Notice physically this means that a change in flux must be compensated by the motion of the interface. Similarly we could define a coordinate frame centered around η = ηˆ in which case we must impose a strict continuity of flux, however, the motion of the interface will be included in the relative velocity term and the same conclusions will apply. If we were to apply the same analysis for the oil flux across the boundary we would get that, [ˆ x · qo − v(1 − s)] = −δ α¯o  dP0 dP1 + ε −δ α¯o exp (iny + σt) dη dη  = vs∗ + εs∗ σ A˜ exp (iny + σt) .  57  3.1. 1D Model in the Absence of Capillary Pressure If we sum the expressions for the jump in individual fluxes we get the jump in total flux which satisfies, ¯ dP1 exp (iny + σt) ¯ dP0 + ε λ [ˆ x · (qs + qo ) − v] = λ dη dη = 0. Here we see that the jump in total fluid flux across the boundary is conserved as we expected. From this condition we generate another solvability condition, ¯ dP1 = 0. λ dη Finally we look at the continuity of pressure. On either side of the interface we can expand P (ˆ η k ) with ηˆ ≈ 0 where k = + or − depending on the left or right near ηˆ to get, P (ˆ η k ) ∼ P0k (0) + ε  dPo dη  A˜ exp (iny + σt) + εP1k (0) exp (iny + σt) η=0k  so that we can write the jump in pressure as, dP0 [P ] = [P0 ] + ε [P1 ] + A˜ dη  exp (iny + σt) = 0.  This generates the conditions that, [P0 ] = 0 dP0 [P1 ] = −A˜ dη  58  3.1. 1D Model in the Absence of Capillary Pressure and therefore we get for our three solvability conditions, ¯ dP1 = 0 λ dη dP1 α¯s = s∗ σA dη [P1 ] = −A  dP0 dη  which are precisely the jump conditions (3.14) predicted by [47]. Now we can rewrite these jump conditions in terms of a matrix problem as,       s∗ nσ 1 ¯ λ(0)  −  1 ¯ ∗) λ(s  0     α¯s (s∗ ) α¯s (0) A 0       1 −1  C−  0    ∗ ¯ ¯ C+ 0 λ(s ) λ(0)  and if we wish to have non-trivial eigenfunctions then the determinant of this matrix must be zero. This means that, ¯ ∗ ) − λ(0) ¯ ¯ α¯s (s∗ ) − λ(s ¯ ∗ )α¯s (0) λ(s λ(0) σ = ¯ λ(s ¯ ∗ )(λ(0) ¯ ¯ ∗ )) n + λ(s s∗ λ(0) which we can rearrange to σ =v n  ¯ ∗ ) − λ(0) ¯ λ(s ¯ ∗ ) + λ(0) ¯ λ(s  (3.15)  by recalling that, f (s) =  αs λ  59  3.1. 1D Model in the Absence of Capillary Pressure and v=  f (s∗ ) f (s∗ ) − f (0) = s∗ s∗ − 0  since f (0) = 0. This result states that the front will be unstable whenever the total relative permeability at the shock saturation value is larger than the final state (s = 0). Whenever δ is sufficiently small then this will always be true. Consider the numerator for the eigenvalue, λ(s∗ ) − λ(0) = s∗3 + δ(1 − s∗ )3 − δ. Recall from (3.6) that for δ  1 the solution for the shock saturation value  looks like, s∗ ∼ 21/3 δ 1/3 where we have taken only the first term which will provide us accuracy up to o(δ) but this is all that is required. With this in consideration we get that, s∗3 + δ(1 − s∗ )3 − δ ≈ 2δ > 0 and thus for large viscosity differences (small δ), fingering will always occur. As mentioned, this eigenvalue result is a generalized analysis of a front solution which transitions from s = 1 to s = 0 carried out in [47], however it yet generalizes even further when αs (s∗ ) = αo (0) = 1 to the Taylor-Saffman instability [32]. A unique difference between the Taylor-Saffman instability and our instability problem is that in the Taylor-Saffman problem, an instability occurs whenever δ < 1; however in this problem we find values of δ < 1 for which the system is stable. For example, consider δ = 1/5, if we repeat 60  3.2. 1D Model With Capillary Pressure the analysis we get that s∗ = 1/2. We then get that, λ(s∗ ) − λ(0) = −  1 < 0. 20  Therefore we see in this quick example that there is a δ < 1 producing negative eigenvalues whereby the Taylor-Saffman result would still produce positive eigenvalues. The real inherent difference between this model and the Taylor-Saffman result is that in the latter, there is no dual phase, i.e. the upstream value is s = 1 and the downstream value is s = 0 so the permeabilities are simply take as the permeabilities of water and oil. However, in our case, we have an upstream saturation value different from one or zero and hence there is dual phase, i.e. a certain percentage of pores are filled with water while others with oil, which requires the use of the total permeability. This mixed volume can act to stabilize fingers because at the sharp interface the saturation gradient is smaller than the Taylor-Saffman case so there is less work required in moving the front.  3.2  1D Model With Capillary Pressure  Now we consider the same 1D analogue of (2.9) considered in the previous section but with a non-zero Capillary pressure, Pc = Po − Ps = J(s) where J(s) is the Leverett function. As mentioned, it is chosen so that it is monotonically decreasing (J (s) < 0), i.e. there is a pressure gradient in the 61  3.2. 1D Model With Capillary Pressure direction of the flow field (left to right). With this in mind we shall write, J (s) = −|J (s)| for sign consistency. Having a form for the capillary pressure means we can still write the model solely in terms of one of the pressure terms, P = Po (although the choice is irrelevant, we will shortly justify the use of Po over Ps ) to get that, αs αo ∂s αs − δ |J (s)| , λ λ ∂x δαo δαo αs ∂s qo = + |J (s)| . λ λ ∂x qs =  Now we expect that since  1, we will maintain the behaviour of the  solution to the hyperbolic equation even though the shock is completely regularized. With this consideration we once again rescale, η = x − vt where v is the hyperbolic shock velocity from before. We are then able to write the now parabolic problem for the water saturation as, ∂s ∂ ∂ + (f (s) − vs) = δ ∂t ∂η ∂η  αs αo ∂s |J (s)| λ ∂η  (3.16)  of which we will seek steady-state solutions. Firstly we note that, for η of size O(1), if we were to expand s in powers of  then to first order we just  62  3.2. 1D Model With Capillary Pressure get for the steady state, ∂ (f (s) − vs) = 0, ∂η ∂P 1 =− , ∂η λ  (3.17a) (3.17b)  which is exactly the hyperbolic problem with no capillary pressure and thus we get the shock wave steady state profile (3.9). Since the first order steady state solution satisfies, ds =0 dη then if we continue the expansion we would see that the first order solution would persist at all orders of  and hence we have recovered the full outer  steady-state solution. However, we see that there is a singular perturbation problem or internal boundary layer as η → 0 and so if we rescale our space variable, η z= ,  −∞<z <∞  u(z) = s( z) w(z) = P ( z) then we get that the diffusive terms created from the capillary pressure are of O(1) in this region and we define a new steady-state problem, d d (f (u) − vu) = δ dz dz  αs α0 du |J (u)| . λ dz  63  3.2. 1D Model With Capillary Pressure This boundary layer reformulation is a standard technique for weak viscous shock profiles [42]. Now we can integrate this once (and take the integration constant to be zero so that we match to the outer constant solution at z → ∞) to get a first order problem, du = dz  λ αs αo  f (u) − vu δ|J (u)|  .  (3.18a)  We may also define a steady-state pressure, dw =− dz  1 + αs |J (u)| du dz λ  by using the expressions for qs , qo and noting that qs + qo = 1. We may simplify the pressure expression by utilizing the saturation derivative to get, dw vu − 1 . = dz δαo  (3.18b)  While the choice of the pressure variable which we use is irrelevant because of the capillary relationship, if we had chosen instead to look at the water pressure P = Ps instead of P = Po we would get, dw vu =− dz αs which diverges as s → 0. Since the oil pressure only diverges as s → 1, which is outside the domain of interest, the function is analytic where we are concerned and thus it is more suitable for analysis. We now, in principle, have the solution for water saturation when capillary pressure is included. For a  64  3.2. 1D Model With Capillary Pressure non-specific J(s) we notice that our parabolic problem (3.16) degenerates to a hyperbolic one for all  since αs (0) = 0. This is noted in that the solu-  tion (3.18a) has a divergent derivative at s = 0. It is known ([29],[9]) that parabolic equations which degenerate to hyperbolic equations have solutions with finite wave speed propagation or solutions with compact support and also that they may exhibit “waiting-time” solutions [24]. This means that the diffusive problem will not reach the lower far-field steady-state of s = 0 in an exponential sense but rather in a finite distance. This naturally creates an infinite slope at the point where the travelling wave intersects s = 0. Physically this corresponds to the system attempting to preserve the hyperbolicity and hence shock like behaviour for small s. We can seek the form of this degenerate solution by asymptotically expanding (3.18a) for small u.  If u  1 (more specifically, we require u  δ 1/3 , a reasonable assumption  since we are near u = 0) then, αs (u)αo (u) u3 ∼ λ(u) δ and f (u) ∼  u3 δ  as well and therefore, du 1 ∼ dz |J (0)|  1 v − 2 δ u  where we assume J(u) is smooth enough that J (0) is defined. Since u then it is reasonable to assume that u2  δ 1/3  δ and so we can approximate  65  3.2. 1D Model With Capillary Pressure further, v 1 du ≈ − dz |J (0)| u2 which has solution, u∼  3v (z − a) − |J (0)|  1/3  .  Here a is the arbitrary point at which the travelling wave intersects s = 0. This arbitrary value is due to the translational invariance of travelling waves and thus can be determined uniquely by fixing a point on the wave. Notice that this indeed has the property that u(a) = 0 but that the derivative is singular there. If we choose a specific J(s) to make the degeneracy vanish6 , i.e., |J (s)| =  1 αs  or some multiple there of then we still get for small u that, f (u) ∼  u3 δ  but now, αo (u) 1 ∼ λ(u) δ so that when u  δ 1/3  1 that, du ∼ u3 − δvu. dz  6  While there could be other choices to make degeneracy vanish, we choose the simplest form for our analysis  66  3.2. 1D Model With Capillary Pressure However since u  δ 1/3 then it is reasonable to assume that u3  δu and so  we can approximate further to finally get that, du ∼ −δvu dz so that without degeneracy we have the standard linear decay rate which corresponds to exponential decay, u(z) ∼ exp (−δvz) near u = 0. For completeness we will look how the solution decays to the other constant state, u = s∗ . In this case, 1 f (u) − vu ≈ (f (s∗ ) − vs∗ ) + (f (s∗ ) − v) (u − s∗ ) + f (s∗ ) (u − s∗ )2 2 =0  =0  1 = f (s∗ ) (u − s∗ )2 2 so therefore, ∂u λ(s∗ ) f (s∗ ) ∼ (u − s∗ )2 ∂z αs (s∗ )αo (s∗ ) 2δ|J (s∗ )| for u ≈ s∗ where once again we require J(u) smooth enough so that J (s∗ ) is defined. Note that in this case, u(z) is not degenerate as the upstream decay is algebraic so that the function is analytic with respect to u. If the steady-state solution were a generic Riemann problem shock solution (i.e. it did not originate from a rarefaction wave) then we would expect the standard exponential decay, akin to the solution for u near zero in the non-degenerate case. This is because the linear term in the Taylor expansion would no longer 67  3.2. 1D Model With Capillary Pressure vanish. We conclude this section by noticing that in both cases (ignoring degeneracy), the limiting behaviour matches the outer solution, i.e., du ds = 0 = lim z→±∞ dz η→0 dη dP 1 dw lim = − = lim . η→0 dη z→±∞ dz λ lim  In the degenerate case, we still match to the outer solution but at a finite z. Figure 3.3 plots the non-degenerate saturation steady state while Figure 3.4 plots the corresponding pressure derivative. Saturation Steady-State:Boundary Layer 0.16  0.14  0.12  0.1  u  0  0.08  0.06  0.04  0.02  0  -0.02 -10  -8  -6  -4  -2  0 z  2  4  6  8  10  Figure 3.3: Saturation steady-state solution for the inner boundary layer centered around η = 0. It is a smooth heteroclinic orbit connecting the steady-state values of the outer hyperbolic problem.  68  3.2. 1D Model With Capillary Pressure Pressure Derivative Steady−State:Boundary Layer −2  −2.5  −3  dw0 dz  −3.5  −4  −4.5  −5  −5.5 −10  −8  −6  −4  −2  0 z  2  4  6  8  10  Figure 3.4: Pressure derivative steady-state solution for the inner boundary layer centered around η = 0. The derivative has the appropriate limiting behaviour to the outer hyperbolic problem. The saturation is solved numerically using ODE45 in MATLAB and the pressure derivative is solved using this solution along with (3.18b). We see from the figures that we do indeed have the correct matching behaviour. Since we now have an understanding of the behaviour of the steady-state solution with capillary pressure, we will perform a stability analysis as was the case when there was no capillary pressure to determine if the steady state is stable to transverse disturbances.  69  3.2. 1D Model With Capillary Pressure  3.2.1  Linear Stability Analysis  We consider the steady-state solutions for the outer problem s0 and P0 (3.17) as well as u0 and w0 for the inner problem (3.18) to be a 1D cross section of a 2-dimensional problem like with the steady-state shock solution. Therefore, like before, we wish to see how the wave front stabilizes in the presence of transverse wave disturbances. We once again form a stability expansion as, s ∼ s0 + εs1 (η) exp (iny + σt) , P ∼ P0 + εP1 (η) exp (iny + σt) , s1 ,  dP1 ds1 , P1 , → 0, dη dη  η=0 η=0  |η| → ∞  where as before, the conditions are chosen to have smooth decay to the base state for large |η| and ε is used to control term balance. We now wish to substitute this into the full PDE (2.9), ∂s − ∇ · (qo + vsˆ x) = 0 ∂t qs = −αs ∇(P − J(s)) qo = −δαo ∇P ∇ · (qs + qo ) = 0 where once again we use the travelling wave gradient and are using the oil flux in the hyperbolic equation for convenience since we did our analysis using the oil pressure. If we do this substitution then we get to linear order  70  3.2. 1D Model With Capillary Pressure in ε, d dη  ¯ dP1 + λ ¯ dP0 s1 λ dη dη  +  d2 (α¯s |J (s0 )|s1 ) dη 2  ¯ 1 + α¯s |J (s0 )|s1 = 0 −n2 λP d dη  δ α¯o  dP0 dP1 + δ α¯o s1 − vs1 dη dη  − n2 δ α¯o P1 = −σs1  (3.19a) (3.19b)  where once again, overbars indicate a function evaluated at the steady state and prime indicates differentiate with respect to s. Notice that when  = 0,  with the use of the identity δαo = λ − αs , we indeed recover the growth rate for small disturbances in the absence of capillary pressure (3.11) as expected. Outer Solution Now, if we consider η to be O(1) then in this region we take for the steady state (3.17), ds0 =0 dη 1 dP0 = −¯ dη λ and if we proceed in an expansion sense for s1 , P1 , and σ as, s1 ∼ s10 + s11 + . . . P1 ∼ P10 + P11 + . . . σ ∼ σ0 + σ1 + . . .  71  3.2. 1D Model With Capillary Pressure then we get to first order the eigenvalue problem for = 0 and thus, s10 = 0   C0+ exp (−nη) , η > 0 . P10 =  C0− exp (nη) , η<0 Now in the = 0 case we proceeded to find the eigenvalue via jump conditions across the shock interface. However, with  = 0 we no longer have a true  sharp interface and so we cannot strictly impose the same conditions for the eigenvalue. Because s10 = 0, the problem at O( ) is the same problem as at O(1) and thus the solution is unchanged. This type of recursion occurs so that in the outer solution, the eigenvalue problem when  = 0 holds at all  orders. Thus we get for the outer solution, s1 = 0  (3.20a) ∞ p  P1 =  Cpk  exp (−knη)  (3.20b)  p=0  where k = + or k = − to differentiate between the solutions for η > 0 and η < 0. Inner Solution, Non-Degenerate Case Like with the steady-state solution we are motivated to pose an inner variable expansion for the linear stability analysis to analyze the region where the capillary-driven front dominates the steady-state profile. Therefore, we  72  3.2. 1D Model With Capillary Pressure expand around the front η = 0 and let, η z= ,  −∞<z <∞  u(z) = s( z) w(z) = P ( z). We consider the problem for the non-degenerate steady-state solution so that we have analyticity throughout the domain. We therefore will assume that, α¯s |J (u0 )| = 1 and we take as our steady-state solution in this region, (3.18), ¯ du0 λ (f (u0 ) − vu0 ) = dz δ α¯o dw0 vu0 − 1 d w˜0 = = dz δαo dz where we write the pressure derivative as such so that we can explicitly see the order of all terms. While this assumption on the Leverett function is motivated for simplification of analysis, it is not entirely without its physical merits. The experimental Leverett function in [25] has an asymptote as s = 0 consistent with our choice. With this inner coordinate transformation we are  73  3.2. 1D Model With Capillary Pressure able to write our localized eigenvalue problem as, d dz  δ α¯o  dw1 d w˜0 + δ α¯o u1 − vu1 − dz dz d dz  2 2  n δ α¯o w1 = − 2 u1 σ  d w˜0 dw1 d2 ¯ ¯ λ + λ u 1 + 2 ( u1 ) dz dz dz ¯ 1 + u1 = 0 − 2 n2 λw  where once again overbars indicate functions evaluated at the base state. Now in order to comply with matching into the outer region, we wish for the eigenfunctions to be of the same order and so we scale, u1 =  uˆ1  so that we can write as our final inner variable eigenvalue problem, d dz  δ α¯o  dw1 d w˜0 + δ α¯o uˆ1 − vˆ u1 dz dz d dz  −  2 2  n δ α¯o w1 = − uˆ1 σ  (3.21a)  dw1 ¯ d w˜0 d2 ¯ λ +λ uˆ1 + 2 (ˆ u1 ) dz dz dz ¯ 1 + uˆ1 = 0. − 2 n2 λw  (3.21b)  Another motivation for this scaling is that based on the outer solution s1 we seek u1 to be a homoclinic orbit. For any properly normalized homoclinic orbit, h(x), we can define the Dirac delta function as, 1 x δ(x) = lim h →0  74  3.2. 1D Model With Capillary Pressure which is exactly the scaling utilized for our problem. This could also be expected because in the situation with no capillary pressure we saw that the saturation eigenfunction could be expressed as a delta function (3.13) which is essentially the same as imposing the correct jump conditions across the layer [47].  3.2.2  Numerical Solution  We are interested in the numerical solution to (3.21) so that we have a foundation for which to compare any analytic results as well as to have a qualitative understanding of the full solution. We use a cell centered finite difference scheme in order to discretize the problem, i.e. for a set of nodes with z values, we place all information relating to uˆ1 and w1 on the spaces between the nodes and then all derivative or flux information is placed on the nodes. For example, to discretize, d dz  ¯ dw1 λ dz  we denote, ¯ dw1 J =λ dz as a flux term and place it on the nodes as in Figure 3.5.  75  3.2. 1D Model With Capillary Pressure  uˆ1(i−1/2) w1(i−1/2) ¯ i−1/2 λ  xi Ji  uˆ1(i+1/2) w1(i+1/2) ¯ i+1/2 λ  xi+1 Ji+1  α¯oi+1/2  α¯oi−1/2  uˆ1(i+3/2) w1(i+3/2) ¯ i+3/2 λ α¯oi+3/2  Figure 3.5: An example of a discretization of the saturation and pressure eigenfunctions. The thick line represents the actual discretization with dark circles representing nodes and the white circle representing cell-centers. The flow chart above the line shows the symmetry of the cell-centered discretization and how information flows to describe the i + 1/2 cell. The term we are discretizing then is, dJ dz which we can discretize so that on the ith node we get, Ji+1 − Ji dJ = dz h where h is the distance between adjacent nodal points (or adjacent cell points). Since we are using a cell-centered discretization, we use a short difference operator on the nodes so that the two flux terms share a cell-center. By definition, ¯ i dw1(i) Ji = λ dz but since λ and w1 are function of uˆ1 and w1 then they are only defined on ¯ we the cell-centers and thus we must extrapolate them to the nodes. For λ  76  3.2. 1D Model With Capillary Pressure simply take an average, ¯ ¯ ¯ i = λi−1/2 + λi+1/2 λ 2 and for the derivative term we once again use a short difference across the two adjacent cell-centers to extract information solely for Ji , dw1(i) w1(i+1/2) − w1(i−1/2) = . dz h The differences are always set-up so that there is a cell-centered symmetry with each cell (Figure 3.5) and such a discretization is second order accurate, i.e. it has errors that are O(h2 ). We can use our information about the discretizations to get that, d ¯ dw1 λ = dz dz ¯ i+3/2 + λ ¯ i+1/2 )(w1(i+3/2) − w1(i+1/2) ) − (λ ¯ i+1/2 + λ ¯ i−1/2 )(w1(i+1/2) − w1(i−1/2) ) (λ . 2h2 ¯ = 1 we recover the standard second-order centered operator, Notice if λ w1(i+3/2) − 2w1(i+1/2) + w1(1−1/2) d2 w1 = . 2 dz h2 We could discretize the other differential operators in a similar fashion to generate a generalized numerical eigenvalue problem, Au = σBu  77  3.2. 1D Model With Capillary Pressure where the vectors u are, u = [uˆ1 , w1 ]T . Boundary Conditions Now in order to solve the eigenvalue problem we require homogeneous boundary conditions so that we do not add any inhomogeneous terms to our righthand side, eliminating the eigenvalue problem. To accomplish this goal for the pressure we notice that we seek an exponentially decaying solution in the far-field and so we impose a Robin condition there, dw1 = − nw1 , dz dw1 = nw1 , dz  z→∞ z → −∞.  Notice that since the domain endpoints are nodal points, our far-field conditions apply on the nodes and so to resolve this we create what are known as ghost points, artificial cell centers that are introduced in order to use the far-field conditions. If we create an artificial point x1−1/2 to the left of the first node x1 then the far-field condition tells us that,  w1(1−1/2) =  1− 1+  nh 2 nh 2  w1(1+1/2) .  To determine this, we once again used short differences for derivative terms and averaging for non-derivative terms. We can create a ghost point on the right after the final node, the N th node, xN +1/2 and the far field condition  78  3.2. 1D Model With Capillary Pressure imposes that, w1(N +1/2) =  1− 1+  nh 2 nh 2  w1(N −1/2) .  With these values determined, the original equations will hold on the first and last cell center by employing the use of these ghost points. The saturation eigenfunction decays to zero in the far-field and so we can simply impose that7 , uˆ1(1+1/2) = 0 uˆ1(N −1/2) = 0. With the boundary conditions complete we use the eig command in Matlab to compute eigenvalues and eigenfunctions and we notice that all of the numerically computed eigenvalues are real. We take the largest eigenvalue for analysis. This is motivated for two reasons; first, the most unstable mode is always the one observed experimentally and second, in the  = 0 case,  when σ > 0 we had a unique eigenvalue and this is the eigenvalue we would like to extract since, as stated, highly unstable modes are the ones observed in nature. Figure 3.6 and 3.7 show the eigenfunctions for n = 1. 7 While technically it is necessary to also use the ghost point ideology for this eigenfunction as well, we expect that the saturation eigenfunction decays to zero quickly enough that we can state with impunity that the far-field condition holds on the cell center itself rather than the node.  79  3.2. 1D Model With Capillary Pressure Saturation Eigenfunction for n=1 1.4 1.2 1  uˆ1  0.8 0.6 0.4 0.2 0 −0.2 −10  −8  −6  −4  −2  0 z  2  4  6  8  10  Figure 3.6: Numerical solution of uˆ1 in (3.21) for n = 1. To compute the solution we used 201 evenly spaced mesh points.  80  3.2. 1D Model With Capillary Pressure Pressure Eigenfunction for n=1 1400  1200  1000  800  w1  600  400  200  0  -200  -400  -600 -10  -8  -6  -4  -2  0 z  2  4  6  8  10  Figure 3.7: Numerical solution of w1 in (3.21) for n = 1. To compute the solution we used 201 evenly spaced mesh points. When n = 1 our numerical code determines that σ = 1.81. Figure 3.8 plots the value of σ for various values of n.  81  3.2. 1D Model With Capillary Pressure Largest Eigenvalue vs. n  60  50  σ  40  30  20  10  0  0  20  40  60  80  100  120  140  160  180  n  Figure 3.8: This figure shows the largest value of σ for any transverse wave number, n. All wave numbers after the terminus n on the plot are negative and thus stabilize. From the figure we notice that, unlike the linear case in the previous section, there is a defined maximum as well as a critical nc , for which all fingering perturbations of higher wave-number stabilize. The maximum wavenumber is what we would expect to see in actual fingers since it is the most unstable, and therefore most dominant mode. The usefulness of the scaling argument with the inner and outer solution is very apparent when computing the numerical problem. While the properties of σ in the presence of capillary effects has been noted, often in the literature ([34]) a highly robust scheme is required to resolve sharp changes around the saturation interface. However, via our asymptotic scaling argument, the inner problem is always smooth 82  3.2. 1D Model With Capillary Pressure regardless how small  is and thus is handled quite easily through standard  techniques. We conducted a numerical convergence study to test robustness by first letting h → 0, i.e. by increasing the number of grid points over a fixed domain and secondly by extending the far-field conditions for large |z| values. In both cases, we noticed that the results remained unchanged within 1% of the largest computed value and thus we have confidence in our results.  3.2.3  Asymptotic Analysis  We now wish to asymptotically determine the eigenvalue σ and thus we proceed in the usual sense by expanding the eigenfunctions and eigenvalue of (3.21) as, uˆ1 ∼ uˆ10 + uˆ11 + . . . w1 ∼ w10 + w11 + . . . σ ∼ σ0 + . . . . To first order we get that, d dz  δ α¯o  dw10 d w˜10 + δ α¯o uˆ10 − vˆ u10 dz dz  d dz  2 ¯ dw10 + λ ¯ d w˜0 uˆ10 + d (ˆ λ u10 ) = 0. dz dz dz 2  =0  (3.22a) (3.22b)  This is valid as long as n is not sufficiently large (more specifically as long as n is not O(1/ )). In order to be able to find a solution for (3.22) we need to furnish matching conditions to the outer problem. Recall that the outer  83  3.2. 1D Model With Capillary Pressure problem admits solutions, s1 = 0 ∞  p  Cp+ exp (−nη) , η > 0  p=0 P1 = ∞  p   Cp− exp (nη) , η<0 p=0  and so if we Taylor expand around η = 0 we generate the following matching conditions, uˆ10 → 0, |z| → ∞  (3.23a)  w10 → C0+ , z → ∞  (3.23b)  w10 → C0− , z → −∞.  (3.23c)  Since the outer solution for the saturation eigenfunction vanishes smoothly, we assume that the derivatives of uˆ10 vanish in the outer limit as well. Before proceeding we notice immediately from (3.22) and its matching conditions (3.23) that we furnish two solvability equations for the system. Firstly by integrating (3.22a) over the entire z domain and applying the matching conditions we get, dw10 δ α¯o dz  ∞  = 0.  (3.24a)  −∞  If we integrate (3.22b) and use the matching conditions we get that, ¯ dw10 λ dz  ∞  = 0.  (3.24b)  −∞  84  3.2. 1D Model With Capillary Pressure Now to solve the problem we take a first integral of (3.22) to get, ¯ d w˜0 uˆ10 + d uˆ1 = B0 ¯ dw10 + λ λ dz dz dz dw10 d w˜0 δ α¯o + δ α¯o uˆ10 − vˆ u10 = D0 dz dz where B0 and D0 are constants of integration. We can rearrange the second equation to get, dw10 = dz  v δ α¯o D0 − (vu0 − 1) uˆ10 + 2 δ α¯o (δ α¯o ) δ α¯o 1 d vu0 − 1 D0 = du0 uˆ10 + dz δ α¯o δ α¯o dz  1 d2 w˜0 D0 = du0 uˆ10 + 2 dz δ α¯o  (3.25)  dz  where we have used the knowledge that for any g(u(z)), g =  dg 1 dg = du0 . du0 dz dz  Notice firstly we can integrate this expression to get another solvability condition, w10 |∞ −∞  ∞  = −∞  1 d2 w˜0 D0 u ˆ + 10 du0 dz 2 δ α¯o  dz  (3.26)  dz  85  3.2. 1D Model With Capillary Pressure but also that we can solve solely for the saturation eigenfunction, 1 du0 dz  ¯ d w˜0 d2 w˜0 d λ ¯ λ 2 + dz dz dz  uˆ10 +  ¯ d uˆ10 λ = B0 − D0 . dz δ α¯o  (3.27)  After applying the matching conditions we see that, B0 = D0 = 0 is the only valid possibility8 . Therefore we can rewrite (3.27) in terms of a simple first order differential equation, d uˆ10 1 d = − du0 dz dz dz  ¯ d w˜0 λ dz  uˆ10 .  In order to proceed further we notice that we can relate the base state derivatives via, du0 ¯ d w˜0 + 1 =− λ dz dz so that we get, d uˆ10 1 d2 u0 d = du0 = ln dz dz 2 dz dz  du0 dz  8 It may seem audacious to state this result since u ˆ10 is dividing by the derivative of the base-state which itself decays to zero. We then merely assume that the saturation eigenfunction must decay as quickly as the base state derivative. By not doing so we either have an unbounded eigenfunction or a trivial one, neither of which is permitted.  86  3.2. 1D Model With Capillary Pressure which has solution, uˆ10 = A0  du0 . dz  (3.28)  We can substitute this back into (3.25) to solve for w10 , w10 = A0  d w˜0 + E0 dz  (3.29)  where we can use the matching conditions to get that the constant is,  E0 =  ¯ λ(0) 1− ¯ ∗ λ(s )  −1  ¯ λ(0) C− − ¯ ∗ C+ . λ(s )  We can now simplify solvability conditions (3.24) and (3.26), δ α¯o  dw10 dz  ¯ dw10 λ dz  ∞  =0  (3.30a)  =0  (3.30b)  −∞ ∞ −∞  w10 |∞ −∞ = A0  d w˜0 dz  ∞  (3.30c) −∞  which are indeed satisfied with our solution for any constants A0 , C0+ and C0− .  87  3.2. 1D Model With Capillary Pressure Determining σ0 We are now able to look at the next order, O( ) to get that our eigenvalue problem is, d w˜0 dw11 + δ α¯o uˆ11 − vˆ u11 dz dz  d dz  δ α¯o  d dz  dw11 ¯ d w˜0 d2 ¯ λ +λ uˆ11 + 2 (ˆ u11 ) = 0. dz dz dz  = −σ0 uˆ10  (3.31a) (3.31b)  In order to generate matching conditions at this order we must once again Taylor expand our outer solution, (3.20) to get that, uˆ11 → 0,  |z| → ∞  (3.32a)  w11 → C1+ − nC0+ z,  z→∞  (3.32b)  w11 → C1− + nC0− z,  z → −∞.  (3.32c)  However, in what follows we will set C1+ = C1− = 0 since, like with the first order problem, any arbitrary values will suffice until we analyze the next order problem. Like in the first order case, we immediately generate two solvability conditions from integrating (3.31) over the domain and using the matching conditions (3.32), dw11 δ α¯o dz ¯ dw11 λ dz  ∞  ∞  uˆ10 dz = A0 σ0 s∗  = −σ0  (3.33a)  −∞  −∞ ∞  = 0.  (3.33b)  −∞  88  3.2. 1D Model With Capillary Pressure To solve (3.31) we take a first integral to get, δ α¯o  d w˜0 dw11 + δ α¯o uˆ11 − vˆ u11 = −A0 σ0 u0 + D1 dz dz ¯ d w˜0 uˆ11 + d uˆ11 = B1 ¯ dw11 + λ λ dz dz dz  where B1 and D1 are constants of integration. We can rearrange the first equation to get, dw11 1 d2 w˜0 A0 σ0 D1 = du0 uˆ11 − u0 + 2 dz dz δ α¯o δ α¯o  (3.34)  dz  which we can use to solve entirely for the saturation eigenfunction, d uˆ11 1 − du0 dz dz  d2 u0 dz 2  uˆ11 = B1 +  ¯ ¯ 0 σ0 λ λA u0 − D1 . δ α¯o δ α¯o  (3.35)  Notice in both of these simplifications we have recycled ideas from the previous section regarding relationships between steady states. Now when we apply the matching conditions we get that the constants must satisfy, B1 = D1 =  A0 σ0 v  where we have used the fact that, 1 s∗ = . f (s∗ ) v  89  3.2. 1D Model With Capillary Pressure We can then write the solution for (3.35) as, uˆ11 = A1  du0 A0 σ0 du0 + dz v dz  However, 1+  1+  ¯ λ (vu0 − 1) δ α¯o  1 du0 dz  dz.  ¯ λ du0 (vu0 − 1) = − δ α¯o dz  and so finally we get that, uˆ11 = A1  du0 A0 σ0 du0 − z. dz v dz  (3.36)  We can substitute this into (3.34) to get that, d2 w˜0 A0 σ0 dw11 = A1 − dz dz 2 v = A1  z  d2 w˜0 A0 σ0 d − dz 2 v dz  d2 w˜0 vu0 − 1 + dz 2 δ α¯o z  d w˜0 dz  where we have used that, d w˜0 vu0 − 1 . = δ α¯0 dz Finally then we get that, w11 = A1  d w˜0 A0 σ0 d w˜0 − z + E1 dz v dz  (3.37)  90  3.2. 1D Model With Capillary Pressure and after using the matching conditions (3.32), we conclude that9 , E1 = 0 A0 σ0 ¯ ¯ ∗ = −λ(0)nC + = λ(s )nC− . v Notice this last statement seems ill-posed unless the two constants are equal but this is resolved in the solvability condition (3.33b) where we have that, ∞  ¯ dw11 λ dz  =0 −∞  which after using the matching conditions (3.32) becomes, ¯ ¯ ∗ −λ(0)nC + = λ(s )nC− , precisely the condition we need. Having the solutions uˆ11 and w11 we can write the complete solvability conditions at this order as, δ α¯o  dw11 dz  = A0 σ0 s∗  (3.38a)  =0  (3.38b)  (w11 )∞ −∞ = 0.  (3.38c)  ¯ dw11 λ dz  9  ∞ −∞ ∞ −∞  This conclusion for E1 is based on the arbitrary choice of setting the constants C1k = 0.  91  3.2. 1D Model With Capillary Pressure Finally we can write the solution up to O( ) as, du0 A0 σ0 du0 du0 + A1 − z dz dz v dz d w˜0 A0 σ0 d w˜0 w1 = A 0 + E0 − z dz v dz uˆ1 = A0  (3.39a) (3.39b)  and we can combine the solvability conditions from the first order problem (3.30) and the second order problem (3.38) to get, dw1 δ α¯o dz ¯ dw1 λ dz  ∞  = A0 σ0 s∗ −∞ ∞  =0 −∞  (w1 )∞ −∞  = A0  d w˜0 dz  ∞  . −∞  If we then scale back to the outer variables, let A˜0 = −A0 , and use the identity λ = αs + δαo , we get, ¯ dP1 = 0 λ dη dP1 α¯s = A˜0 σ0 s∗ dη dP0 [P1 ] = −A˜0 dη  (3.40a) (3.40b) (3.40c)  where [·] indicates a jump across η = 0 like before. These are in fact the same jump conditions (3.14) for the  = 0 problem and thus we conclude  92  3.2. 1D Model With Capillary Pressure that, ¯ ∗ ) − λ(0) ¯ λ(s σ0 = nv ¯ ∗ ¯ λ(s ) + λ(0)  (3.41)  as we would have expected from the onset. Knowing σ0 we can use (3.40) to get that for any A˜0 , we require, ¯ λ(0) A˜0 1+ ¯ ∗ C+ = − ¯ λ(0) λ(s ) ¯ A˜0 λ(0) C− = ¯ ∗ 1 + ¯ ∗ λ(s ) λ(s )  −1  −1  ¯ λ(0) 1− ¯ ∗ λ(s ) ¯ λ(0) 1− ¯ ∗ . λ(s )  The fact that we are still left over with an arbitrary constant is not unsurprising since homoclinic orbits have a translational invariance associated with them and thus for unique solutions, the eigenfunction requires some sort of prescribed condition. While we could define a value for uˆ1 (0), we instead impose a normalizing condition, ∞  uˆ1 dz = 1. −∞  This supports the eigenfunction approaching a Dirac-delta function as → 0. If we impose such a normalization condition we get that, ∞  uˆ10 dz = −A0 s∗ = 1 −∞  and so, A0 = −  1 . s∗  93  3.2. 1D Model With Capillary Pressure At the next order then we get that, ∞  uˆ11 dz = 0 −∞  which requires that, ∞  A1 −∞  du0 du0 −z dz dz  dz = 0.  If we integrate the second term by parts and notice that since u0 is a travelling wave then for large arguments,  u0 ∼    s∗ z, z < 0  0,  z>0  then the second term vanishes and the above condition is indeed satisfied if A1 = 0. Notice again that the inclusion of non-zero constants C1+ and C1− in the far-field has no effect on the solvability condition or eigenvalue at this order. We can use our finite difference solver to compare our asymptotic solution with the numeric one. In order to impose that our numerical saturation eigenfunction obeys our solvability condition, we compute it and then numerically integrate over the domain. We then divide our entire eigenvector by this result and this ensures that the saturation eigenfunction has the required property that, ∞  uˆ1 dz = 1. −∞  There is a subtleness involved when attempting to compare the code and the asymptotic result. Our asymptotic assumption was that the solution in the 94  3.2. 1D Model With Capillary Pressure far-field is linear if we consider terms up to O( ). However, in order to have a well posed numeric problem, we required homogeneous conditions and thus we imposed that the solution was exponential in the far-field. Thus, automatically, the numeric solution includes all terms in the Taylor expansion of the outer solution (3.20). As n increases, the higher order Taylor series terms become significant and thus the solution is no longer linear. Based on our value of  we choose n = 1 since the next Taylor series term is small. Figure  3.9 shows the numerical saturation eigenfunction computed using the finite difference solver along with the uˆ10 and uˆ11 computed asymptotically. Figure 3.10 shows the same result but for the pressure saturation eigenfunction. Saturation Eigenfunction for n=1 1.4 1.2 1 Numeric u10  0.8 uˆ1  u10+εu11 0.6 0.4 0.2 0 −0.2 −10  −8  −6  −4  −2  0 z  2  4  6  8  10  Figure 3.9: Saturation eigenfunction computed numerically (solid line), asymptotically to first order (dot-dash line) and to second order (dash-line) for n = 1. 201 mesh points were used in the solution.  95  3.2. 1D Model With Capillary Pressure Pressure Eigenfunction for n=1 1500  1000 Numeric w10  500 w1  w10+εw11  0  −500  −1000 −10  −8  −6  −4  −2  0 z  2  4  6  8  10  Figure 3.10: Pressure eigenfunction computed numerically (solid line), asymptotically to first order (dot-dash line) and to second order (dash-line) for n = 1. 201 mesh points were used in the solution. Recall that numerically we computed that σ0 = 1.81 while asymptotically we estimate that σ0 = 1.85, a good agreement. Large n solution If n is large, i.e., O(1/ ) then terms previously neglected in (3.21) are no longer inadmissible. With this in mind we scale, n ˆ= n σ ˆ= σ  96  3.2. 1D Model With Capillary Pressure and we get a completely O(1) problem, free of , d dz  δ α¯o  dw1 d w˜0 + δ α¯o uˆ1 − vˆ u1 − n ˆ 2 δ α¯o w1 = −ˆ u1 σ ˆ dz dz d dz  ¯ dw1 + λ ¯ d w˜0 uˆ1 λ dz dz  (3.42a)  d2 + 2 (ˆ u1 ) dz  ¯ 1 + uˆ1 = 0. −ˆ n2 λw  (3.42b)  We can use our finite difference solver and compute the largest eigenvalue and eigenfunction for each n ˆ . This is plotted in Figure 3.11. Largest Eigenvalue vs. n ˆ  0.6  0.5  σˆ  0.4  0.3  0.2  0.1  0  0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  ˆn  Figure 3.11: Eigenvalue plot for various values of n ˆ = n. Notice this plot is the exact same as Figure 3.8 only scaled by the appropriate factors above. This should be the case since both plots are solutions to the same problem just with different scaling considerations. It was previously 97  3.2. 1D Model With Capillary Pressure mentioned that the inner solution scaling is advantageous because it ensures a smooth problem regardless of . This formulation extends that idea further; by considering a strictly O(1) problem;  has truly been removed and  the numerical methods will always be solving a smooth problem. Therefore but of course, for fixed n, n ˆ → 0 as  Figure 3.11 holds for all values of → 0. This means that as  → 0, the linear eigenvalue persists for larger  values of n as expected. The universality of this formulation being independent of  is that now, the analysis would hold for any value of the inverse  capillary number, large or small. As the inverse capillary number gets larger then for a fixed n ˆ , the actual value of the wave number, n, decreases. This means that for large capillary dominating effects, even small wave numbers can produce stable fingering patterns and as → ∞, then all wave-numbers produce stable fingering patterns, an interesting result although not entirely unsurprising since  → ∞ is equivalent to large surface tension which is a  strong stabilizing mechanism.  There is a caveat associated with this formulation. If the remaining asymptotic terms to (3.21) are algebraic with respect to n, i.e. they can be written in the form,  ∞ p  σ∼  σp np+1  p=0  and if we rescale using the large n formulation as, n ˆ= n σ ˆ= σ  98  3.3. 1D Temperature Model then,  ∞  σp np+1 .  σ ˆ∼ p=0  Therefore we see that the asymptoticness of the original series is lost and all of the terms are of equal importance. Unfortunately, this does not mean that having the full asymptotic series determines σ in the large n case nor viceversa because the scaled large n problem includes all terms in (3.21) where as the asymptotic analysis surrounding n not sufficiently large neglects certain terms. Thus the asymptotic series obtained, cannot be scaled to determine σ for large n.  An important note for both numerical eigenvalue problems is that, while we sought the most positive (unstable) eigenvalue on the basis that it would be the most commonly observed mode, only one positive eigenvalue ever persisted. This means that for all n, there really is only one unstable mode, an interesting result.  3.3  1D Temperature Model  With the inclusion of temperature we will once again restrict ourselves to a one-dimensional domain and we will neglect capillary pressure for simplicity. We begin by looking at (2.21h), ∂T ∂ ˆ ∂T ∂T + (Cs Cˆs (s)qs + Co Cˆo (s)qo ) = D (D(s) ) ∂t ∂x ∂x ∂x  99  3.3. 1D Temperature Model where, recalling from Table 2.2 we have that, Cs = 0.405 Co = 0.105 D = 7.24 × 10−10 . From here it is clear that diffusion is negligible and so we will discount it from the model10 . Now we observe that the convection numbers are of the same order and reasonably close to 1 (within an order one factor). With this in mind we will take for simplicity that, Cs Cˆs = Co Cˆo = 1, which essentially amounts to assuming the volumetric heat capacities, ρcp are equal. Now since we are analyzing the 1D model, the conservation of volume (2.4) once again amounts to, qs + q o = 1 so that we are able to simplify the temperature dependence to simply being the hyperbolic partial differential equation, ∂T ∂T + = 0, ∂t ∂x  − ∞ < x < ∞.  √ There will be a boundary layer of width O( Dt) which means that as time progresses the temperature front will become smooth but this is on a long time scale and we neglect it. 10  100  3.3. 1D Temperature Model We will take for initial conditions to this problem, Riemann data such that there is a “hot” front coming in from the left, x < 0 with T = 1 and then everything on the right x > 0 is “cold” with T = 0,  T (x, 0) =    1, x < 0  .   0, x > 0 With this consideration then we can easily solve the temperature profile to get that,  T (x, t) =    1, x < t  .   0, x > t We now have solved for the temperature and thus can easily compute the viscosity. We have from (2.21d) that, µo (T ) = (1 − δ) exp(1 − T ) + δ however upon the knowledge that for our particular problem δ  1 we will  simply ignore δ in the viscosity. Therefore using the temperature we can get that the oil viscosity is approximately,  µo (T ) =    1,  x<t   exp(1), x > t  101  3.3. 1D Temperature Model and of course we have that the water viscosity is µs (T ) ≡ 1. Having defined this we can now write our velocities as, dP dx αo dP qo = −δ µo (T ) dx qs = −αs  which when combined with, qs + qo = 1 can define the water (or oil) velocity solely in terms of s and T as, qs (s, T ) =  αs (s) Λ(s, T )  where we define Λ as the modified total mobility, Λ(s, T ) = αs (s) + δ  αo (s) . µo (T )  For our currently specified temperature we have that,  Λ=    αs + δαo = λ,  x<t   αs + δ exp(−1)αo , x > t and therefore we have that,  qs =      s3 s3 +δ(1−s)3  = f1 ,  x<t     s3 s3 +δ exp(−1)(1−s)3  = f2 , x > t  .  102  3.3. 1D Temperature Model We now have two equations that describe saturation, one for the left of the temperature interface and one for the right and therefore we must solve, ∂s1 ∂f1 + = 0, ∂t ∂x ∂s2 ∂f2 + = 0, ∂t ∂x  −∞<x<t  (3.43a)  t<x<∞  (3.43b)  where, for initial conditions, we will impose the Riemann data as in the scenario which did not include temperature,  s(x, 0) =    1, x < 0  .   0, x > 0 Notice that this is not a standard characteristics problem as we have a pseudo-free boundary at x = t and so we will have a Stefan condition that describes the dynamics across that boundary. We call it a pseudo-free boundary because it is a boundary that moves regardless of saturation dynamics and thus is free in the saturation context. However, it is not a natural growing boundary such as ice melting, and so it is restricted by the temperature dynamics. For this reason we will call it the internal interface condition. Consider a general interface at time t that moves with a normal velocity vn and then consider that same interface having moved a time t + ∆t as shown in Figure 3.12.  103  3.3. 1D Temperature Model  vn  H hot  cold  t  t + ∆t  Figure 3.12: Visualization of the area change via a moving internal interface in the water saturation problem. Notice that the shaded region in Figure 3.12 represents an area given by Hvn ∆t where H is the height of the shock. In this area there has been a water saturation differential as the cold water that was present in this region at time t is now replaced with hot water at time t + ∆t and thus we can write this saturation differential as, Sd = (s2 − s1 )Hvn ∆t. Now as the front moves there is also a net change in the flux of water over the time period t + ∆t and thus we get a flux differential, Sf = (f2 − f1 )H∆t. Conservation of mass/volume dictate that these differentials must balance  104  3.3. 1D Temperature Model and thus we get the Stefan condition, vn =  f2 − f1 . s2 − s1  (3.44)  Notice that the internal interface condition acts like a Rankine-Hugoniot condition which is not coincidental as the internal interface can be thought of as the shock that forms from intersecting characteristics of the two problems. However, this terminology is being avoided because we truly do have a free interface which characteristics can emanate from, something that is invalid for a standard shock. Now we have derived (3.44) for a general moving interface and typically this would provide information about the location of the interface however, we are in a unique situation in which we already know the velocity of the interface as it is prescribed by the velocity of the temperature front and therefore, vn = 1. We therefore see that the Stefan condition is used to determine discontinuities in the saturation profile.  105  3.3. 1D Temperature Model  3.3.1  Stationary Boundary Characteristics  We can apply the standard method of characteristics to (3.43) along the stationary boundary t = 0. If we consider the problem for s1 then we have, ∂s1 ∂s1 + f1 (s1 ) = 0, ∂t ∂x  −∞<x<t  s1 (x, 0) = 1 where once again prime notation denotes differentiation with respect to s. We are able to state with impunity that the Riemann initial condition will hold because when t = 0 our domain is −∞ < x < 0 which satisfies the requirements. After applying the method we get that along vertical characteristics x = β1 , the solution is, s1 = 1. If we applied the same procedure to the problem for s2 using the other Riemann data condition than we get that the solution is, s2 = 0 along vertical characteristics x = β2 . Note that the problem can no way be complete because we have not utilized any of the characteristic data from the internal interface.  106  3.3. 1D Temperature Model  3.3.2  Internal Interface Characteristics  Now consider solving either of the saturation problems, ∂si ∂si + fi =0 ∂t ∂x without any initial data. We see that the solution is, si = C, a constant, along any characteristics, x = fi (C)t + β. Therefore the solution that originates from the boundary x = t must be of the form si = C, however we do not know which region the boundary generates characteristics for nor what this value C is to be. We know that at the internal interface, dx =1 dt and so if the characteristics are to emerge into the hot region we require that the slope of the characteristics in the hot region to be less than one, i.e., f1 (C) < 1.  107  3.3. 1D Temperature Model Conversely if the characteristics emerge into the cold region then we require that the characteristic slopes exceed one, f2 (C) > 1. Each scenario produces a unique characteristic plot and hence a unique solution and we will investigate both cases. Hot Region Emergent Characteristics If characteristics emerge from the internal interface x = t into the hot region then they will be of the form, x = f1 (C)t + β with f1 (C) < 1. Along these characteristics, s = C. We will then need a rarefaction fan to connect the solution from the vertical characteristics in the hot region (s = 1 on x = β1 ) to these new linear ones. For the rarefaction fan the characteristics will emanate from the origin and thus be a function of χ = x/t. Using this variable transform and substituting it into the saturation equation for the hot region (3.43a), we get that the characteristics satisfy, x = f1 (s)t,  C < s < 1.  Now in the cold region, we do not wish for the x = t boundary to be a true shock since we have characteristics emerging from it and so we place another rarefaction fan to prevent the vertical characteristics in this region (s = 0 on  108  3.3. 1D Temperature Model x = β2 ) from intersecting the internal interface. This rarefaction fan is akin to the scenario without temperature where we introduce it to smooth the characteristics of the two initial conditions. The fan will have characteristics of the form, x = f2 (s)t,  s∗ < s < s+  where, like before, s∗ , the shock saturation value, is the point where a shock forms between the rarefaction and the vertical characteristics. The value s+ is the value where the rarefaction fan begins, i.e. it is such that, f2 (s+ ) = 1. Like in the isothermal case, the shock that forms (x = ξ) must be the fan characteristic associated with s∗ and so we know that, dξ = f2 (s∗ ) dt however the shock condition must also satisfy a Rankine-Hugoniot condition, dξ f2 (s∗ ) = f2 (s∗ ) = dt s∗ which determines s∗ . While there may be multiple solutions for s∗ we will be able to determine it uniquely, discarding any solutions that do not satisfy 0 ≤ s ≤ 1 and s∗ < s+ . Finally, we can determine the constant C via the Stefan condition (3.44) to get that, f1 (C) − f2 (s+ ) = 1. C − s+ 109  3.3. 1D Temperature Model While, like with s∗ , there may be multiple solutions for C, it can be determined uniquely since s+ < C < 1. If there are no valid solutions for C then the characteristics from the internal interface do not emerge into the hot region. Further more, any valid C must also satisfy f1 (C) < 1 as that was an original assumption. Having determined all the constants uniquely we would get that the saturation solution for this case is,    1,       −1   [f1 ]   s(x, t) = C,       [f2 ]−1       0,  x<0 x t  , 0 < x < f1 (C)t f1 (C)t < x < t .  x t  (3.45)  , t < x < f2 (s∗ )t x > f2 (s∗ )t  Figure 3.13 shows the characteristic diagram for this scenario.  110  3.3. 1D Temperature Model s=C s=1  x=0  x=t  x=ξ s=0 t=0  Figure 3.13: Characteristic diagram for saturation when internal interface characteristics emerge into the hot region. Here, solid thin lines represent the constant solution characteristics from the stationary boundaries, dotted lines represent rarefaction characteristics, dashed lines represent the constant solution characteristics from the internal interface, and the dot-dashed line represents the shock. Notice that the internal interface at x = t corresponds to the jump in temperature while the shock solution at x = ξ corresponds to the traditional Buckley-Leverett behaviour. As we will see shortly, this is not the case that our problem is categorized by. As such, we wish to present an artificial case that shows this type of behaviour occurring. If we consider a scenario where viscosities actually decrease with a decrease in temperature then we could imagine that for the same thermal problem, we would solve the same saturation problem as stated but with the viscous term in the cold region being of the form, exp (T − 1)  111  3.3. 1D Temperature Model instead of our current, exp (1 − T ) which has the viscosity increasing. This just places a term exp(1) into the cold water flux, f2 instead of exp(−1) so that, s3 s3 + δ(1 − s)3 s3 . f2 = 3 s + δ exp(1)(1 − s)3  f1 =  Again, this is not a physical case but an easy slight change in the model formulation to show how the hot emergent characteristic case holds. Recall then for this case we need to find a point where f2 (s+ ) = 1 to determine where the rarefaction fan begins. We compute this point numerically to get that, s+ = 0.287. We then compute that the shock saturation value satisfies, f2 (s∗ ) =  f2 (s∗ ) s∗  to get that s∗ = 0.194. Finally we use the Stefan condition to get that C = 0.348. Notice we indeed have that f1 (C) < 1, a requirement of this case and that C > s+ > s∗ > 0. Figure 3.14 plots the numerical solution to the saturation equations with the above flux functions along with an analytic solution in the form of (3.45) using the constants determined.  112  3.3. 1D Temperature Model Temperature influence on water saturation at t=0.0625  1  Numerical Temperature Analytical  0.8  s,T  0.6  0.4  0.2  0  0.45  0.5  0.55  0.6  0.65  0.7 x  0.75  0.8  0.85  0.9  0.95  Figure 3.14: This is an artificial example showing how the hot emergent characteristic scenario holds. The solid line is the analytical solution while the dashed line is the numerical one. The dot-dashed line is the temperature front. To compute the solution numerically, we used a standard shock capturing scheme which handles the solution of the system differently depending on the distance to the discontinuous shock solution. It can be seen from the Figure 3.14 that the scheme still does some smoothing near the points where the derivative is discontinuous.  113  3.3. 1D Temperature Model Cold Region Emergent Characteristics We now consider the situation where characteristics emerge from the internal interface x = t into the cold region. They will be of the form, x = f2 (C) + β with f2 (C) > 1. These characteristics which carry the value s = C will intersect the vertical characteristics with s = 0 forming a classic shock, ξ where the shock speed is given by the Rankine-Hugoniot condition, dξ f2 (C) = . dt C Here the shock saturation value C, unlike s∗ , is determined elsewhere. Now in the hot region, we will need a rarefaction fan to connect the vertical characteristics of s = 1 to the internal interface at x = t. These characteristics will be of the form, x = f1 (s)t,  s− < s < 1  where s− is the value of s where the expansion fan intersects the internal interface and thus it satisfies, f1 (s− ) = 1. We then use our Stefan condition (3.44) to determine C, 1=  f1 (s− ) − f2 (C) . s− − C  114  3.3. 1D Temperature Model Like before, invalid answers for C, as well as those that do not satisfy f2 (C) > 1 indicate that the characteristics do not emerge into the cold region. As with the hot region, with the constants determined, we now have the saturation solution,  s(x, t) =     1,       [f1 ]−1    C,      0,  x<0 x t  , 0<x<t .  (3.46)  t<x<ξ x>ξ  Figure 3.15 shows the characteristic diagram for this case. x=t s=1  x=0  x=ξ s=C  s=0  t=0  Figure 3.15: Characteristic diagram for saturation when internal interface characteristics emerge into the cold region. The thin solid lines represent the constant solution characteristics from the stationary boundaries. The dotted lines represent the characteristics from the rarefaction wave, the dashed lines represent the constant solution characteristics from the internal interface and the dot-dashed line represents the shock. Like in the hot region, the internal interface at x = t corresponds to the jump in temperature while the shock solution at x = ξ corresponds to the traditional Buckley-Leverett behaviour.  115  3.3. 1D Temperature Model  3.3.3  Other Cases  Before continuing forward it is worth noting that we have only discussed the consequences of f1 (C) < 0 or f2 (C) > 0 and have ignored any other possibilities. Firstly we have ignored the unphysical case where f1 (C) > 0(< 0) but f2 (C) < 0(> 0). In this case there are no characteristics emerging from the boundary and there is no solution. This case is also ignored on mathematical grounds. Assume that f1 (C) < 1, in which case, f1 (C) =  δ αo (C)αs (C) − αo (C)αs (C) <1 2 µ01 αs (C) + µδ01 αo (C)  which can be rearranged as, δ (αo (C)αs (C) − αo (C)αs (C)) < µ01  δ αo (C) αs (C) + µ01  2  .  (3.47)  Then, we can write f2 (C) as, f2 (C) =  δ αo (C)αs (C) − αo (C)αs (C) 2 µ02 αs (C) + µδ02 αo (C)  µ01 αs (C) + < µ02 αs (C) +  δ α (C) µ01 o δ α (C) µ02 o  2 2  using inequality (3.47). Now assume T1 > T2 so that µ01 < µ02 and then, αs (C) +  δ α (C) µ01 o  αs (C) +  δ α (C) µ02 o  f2 (C) <  2 2.  116  3.3. 1D Temperature Model Now if C has multiple solutions, the unique value chosen will be the one closest to s = 1 since it must satisfy s+ < s < 1 via the Stefan condition. Therefore it is reasonable that C is near 1 and thus, f2 (C)  αs (C) =1 αs (C)  and thus if f1 (C) < 1 then so too is f2 (C). As δ  1 this approximation  holds stronger even for C not near 1. We can extend these ideas to the other cases easily when δ  1 to show that if fi < 1(> 1) then fj < 1(> 1)  as well so that characteristics will either emerge into the hot or cold region exclusively. Admittedly the analysis doesn’t hold for δ not small and thus we are unable to prove mathematically that the two fluxes share an inequality, however we still physically expect that only one of f1 < 1 and f2 > 1 should hold. Lastly, there is the case when f1 (C) = f2 (C) = 1. In this case the boundary does indeed become the physical shock which is admissible because no characteristics emerge from the boundary. The shock speed of the saturation is at the same velocity as the temperature front; both effects are noted across the one shock. The solution for this problem will be akin to the isothermal problem, it will be a rarefaction wave that turns into a shock.  3.3.4  Solving The Saturation Problem  Now that we’ve analyzed the two interesting physical scenarios that could occur with our problem we must determine which case actually occurs. When the flux functions are not simple as is our case then we must run a small numeric code to either support or validate a case. If we assume that the  117  3.3. 1D Temperature Model characteristics emerge into the hot boundary then we assume that C is such that f1 (C) < 1. If this is the case then we solve for C using the Stefan condition and we get that the only admissible choice for C is, C = 0.045 based on our parameter choices11 . However we notice that f1 (C) = 2.18 > 1 and thus we contradict our assumption. Therefore we can conclude that the characteristics emerge into the cold region. While it is unnecessary to verify this since it is the only other physical option, it is easy to do so since the parameters need to be computed. We solve f1 = 1 to get that, s− = 0.233. Now we use this in the Stefan condition to get that, C = 0.133 which does indeed satisfy f2 (C) > 1. We can now use this to get that the velocity of the post temperature change saturation shock created by the intersection of s = C and s = 0 characteristics is, dξ f2 (C) = = 6.30. dt C 11  This may make our claims above that if f1 (C) < 1 then C is relatively close to 1 seem absurd, however, since δ 1 then C can be sufficiently far from 1.  118  3.3. 1D Temperature Model Figure 3.16 plots the numeric and analytic solution for the saturation problem as posed at t = 1/16. Temperature influence on water saturation at t=0.0625  1  Numerical Temperature Analytical  0.8  s,T  0.6  0.4  0.2  0  0.45  0.5  0.55  0.6  0.65  0.7 x  0.75  0.8  0.85  0.9  0.95  Figure 3.16: Analytic and numerical results for saturation as it depends on temperature. The solid curve is the analytic solution while the dashed curve is the numeric solution. The temperature is plotted as dot-dash in black. The analytic solution is of the form (3.46) with the appropriate constants (solved for numerically) substituted in. The numerical solution is the full solution of the one dimensional equations. To compute it, we used the same shock capturing scheme as with the artificial hot characteristic problem. We 119  3.3. 1D Temperature Model now have the full solution to our problem in the form (3.46). The problem has two shock fronts, one in the classic sense where the jump in saturation is due to the inherent viscosity difference of the two constituent fluids and the other associated with a change in viscosity due to the change in temperature. If we scale our problem in a travelling wave coordinate frame, η = x − vt for some wave speed v then we have a steady state solution if, d (f − vs) = 0 dη where it is assumed here that,  f=    f 1 , x < t  .   f 2 , x > t If, v=  dξ dt  then the coordinate frame is centered around the classical shock front and if we integrate, d (f − vs) dη = f − vs = D dη  120  3.3. 1D Temperature Model for D, a constant, then a steady state exists. We see that for D = 0, f2 (C) = v, C satisfies the condition for a steady state. This is precisely the condition that describes the classical shock front and thus the front is a steady state solution in this travelling wave reference. If instead we let v = 1 then the coordinate frame is centered around the temperature shock front and then by the same justification as the classical front, a steady state solution exists if f1 (s− ) − f2 (C) = 1, s− − C which once again precisely defines the temperature shock and thus this front is a steady state as well. However, each front is only steady with respect to a coordinate transform around its own wave speed and thus they are not mutually steady.  3.3.5  Linear Stability Analysis  Classical Shock We will begin our analysis with the classical shock problem since the temperature is constant across this interface. In this scenario we are centered around the travelling wave coordinate, η = x − vt  121  3.3. 1D Temperature Model with, v=  dξ . dt  Here we have a steady state solution that satisfies, d (f2 − vs0 ) = 0 dη which has an associated pressure steady state, 1 dP0 =− . dη Λ(s0 ) We will do the standard perturbations for s, P (using steady states as above), and the interface ηˆ such that, s ∼ s0 + εs1 (η) exp (iny + σt) P ∼ P0 + εP1 (η) exp (iny + σt) ηˆ ∼ 0 + εA exp (iny + σt) s1 ,  ds1 dP1 , P1 , → 0, dη dη  |η| → ∞  where ε is a small parameter to manage order of magnitude. Notice we do not need to perturb temperature here because the field is isothermal in this region and thus there is no mechanism for changes in temperature. We apply  122  3.3. 1D Temperature Model this to the two dimensional variant of the 1D equations, ∂s + ∇ · (qs − vsˆ x) = 0 ∂t qs = −αs ∇P αo qo = −δ ∇P µo ∇ · (qs + qo ) = 0 where we have once again redefined the gradient operator in the travelling wave as, ∇=  ∂ ∂ , ∂η ∂y  .  When we apply our perturbations to this system we get for our eigenvalue problem, d dη  α¯s  dP1 dP0 + α¯s s1 + vs1 − n2 α¯s P1 = σs1 dη dη d ¯ dP1 ¯ dPo ¯ 1=0 Λ +Λ s1 − n2 ΛP dη dη dη  (3.48a) (3.48b)  which is exactly (3.11) with Λ in place of λ. Once again we let over-bars represent functions evaluated at the base-state. Since the steady-state solutions are constant we can apply a similar methodology to the temperature  123  3.3. 1D Temperature Model independent case to get that, s1 =0, η = 0   B+ exp (−nη) , η > 0 . P1 =  B− exp (nη) , η<0 We then conserve total flux and pressure across the interface to generate the solvability conditions, ¯ dP1 = 0, Λ dη dP1 α¯s = CAσ, dη dP0 [P1 ] = −A , dη which once again are similar to the jump conditions for the temperature independent problem (3.14) except with the modified total permeability and the altered shock saturation value. This can be written as a matrix problem with a determinant that must vanish in order to admit non-trivial solutions for the eigenfunctions. This leads to the solution for the eigenvalue, ¯ ¯ σ Λ(C) − Λ(0) =v¯ ¯ . n Λ(C) + Λ(0) When δ  (3.49)  ¯ the leading term in the 1 then based on the expression for Λ,  numerator looks like C 3 > 0 and so we predict that fingering will occur. Once again we see the similarity to the temperature independent eigenvalue. If we were to use the same flux function here but in the isothermal case, we 124  3.3. 1D Temperature Model would see that the only difference would be the upstream shock saturation value used in the equation. Figure 3.17 shows the linear stability eigenvalue in the absence of capillary pressure with and without the effects of temperature included based on the parameters in Tables 2.1 and 2.2. Linear Stability Eigenvalue-No Capillary Pressure 20  18  16  14    12  10  8  6 Without Temperature With Temperature  4  2  0  0  0.5  1  1.5  2  2.5 n  3  3.5  4  4.5  5  Figure 3.17: Linear stability eigenvalue in the absence of capillary pressure with and without temperature effects. The solid line is the eigenvalue without temperature (3.15) and the dashed line is the eigenvalue with temperature (3.49). Here the upstream temperature change has the effect of promoting fingering across the classical shock.  125  3.3. 1D Temperature Model Temperature Shock We now perform a linear stability analysis around the shock associated with the change in temperature. In this case our travelling coordinate expansion, η = x − vt has v = 1. Our steady state satisfies, d (f − s0 ) = 0 dη with,  f=    f 1 , η < 0  f 2 , η > 0  like before. there is also an associated pressure steady state, 1 dP0 =− dη Λ(s0 , T0 ) but now there is also an associated temperature steady state,  T0 =    1, η < 0  .   0, η > 0 Now as we perform our perturbation analysis, we must consider that the temperature changes across the shock and so we must perturb that variable as well. Also we no longer have a sharp front solely for the saturation variable 126  3.3. 1D Temperature Model but for temperature as well. While, to highest order, these are both centered around η = 0 we have to allow for the possibility that they will degenerate into two local fronts and thus we will have a front variable describing the saturation front, ηˆs and one describing the temperature front, ηˆT . With this in mind we will expand as follows12 , s ∼ s0 + εs1 (η) exp (iny + σt) , P ∼ P0 + εP1 (η) exp (iny + σt) , T ∼ T0 + εT1 (η) exp (iny + σt) , ηˆs ∼ 0 + εAs exp (iny + σt) , ηˆT ∼ 0 + εAT exp (iny + σt) , s1 ,  ds1 dP1 dT1 , P1 , , T1 , → 0, dη dη dη  |η| → ∞  where once again ε is introduced to mange the order of terms. We apply this to the two dimensional variant of the one dimensional model, ∂s + ∇ · (qs − sˆ x) = 0 ∂t qs = −αs ∇P αo qo = −δ ∇P µo ∇ · (qs + qo ) = 0 ∂T ˆ ) · ∇T = 0 + (qs + qo − x ∂t 12  Like with velocity, at first it may seem odd that the temperature dependent viscosity is not being perturbed but indeed it is by incorporating the perturbation of temperature with the explicit form of viscosity  127  3.3. 1D Temperature Model where, like before, we take the gradient to be the derivatives in the travelling wave coordinate. Upon applying the perturbation variables, the eigenvalue problem to linear order is, dP0 dP1 + α¯s s1 + s1 − n2 α¯s P1 = σs1 (3.50a) dη dη   dµ ¯ d  ¯ dP1 ¯ dPo dT dP0 ¯ 1=0 +Λ s1 − δ α¯o 20 T1  − n2 ΛP (3.50b) Λ dη dη dη µ ¯ dη   dµ ¯ ¯ dP0 + 1 dT1 + Λ ¯ dT0 dP1 + Λ ¯ − δ α¯o dT0  dP0 dT0 = σT1 . (3.50c) Λ dη dη dη dη µ ¯2 dη dη d dη  α¯s  with over-bars once again denoting functions evaluated at their base state variables. While we specifically indicate temperature derivatives, primes still denote differentiation with respect to s in order to aesthetically recognize terms that are similar to the temperature independent case. Recall that these equations are valid away from η = 0 where steady-states are well defined. In this case then we can use the fact that, dT0 =0 dη dP0 1 = −¯ dη Λ and substitute it into (3.50c) to get, σT1 = 0  128  3.3. 1D Temperature Model which for non-zero eigenvalues means that, T1 ≡ 0,  η = 0.  This then removes the temperature dependence from the eigenvalue problem and so we can solve (3.50a) and (3.50b) to get the standard solutions, s1 =0, η = 0   B+ exp (−nη) , η > 0 . P1 =  B− exp (nη) , η<0 We must furnish solvability conditions to determine a unique solution for the eigenvalue. In what follows, we use the methodology from the isothermal case by forcing the continuity of physical variables directly. Instead, we could apply a methodology similar to [47] in which case the temperature eigenfunction would be of the form, T1 = AT δ(η). Because we now have two fronts that could form, we must be cautious with how we proceed. We expect that fluid fluxes will be continuous across the front that separates saturation values and that thermal fluxes will be continuous across the front that separates temperature values. Physically we expect that pressure must be continuous along both interfaces. If we impose continuity of fluid flux across the saturation interface as well as continuity of pressure then we generate the standard jump conditions that we’d expect 129  3.3. 1D Temperature Model from the previous problems, ¯ dP1 = 0 Λ dη dP1 α¯s = (s− − C)As σ dη dP0 [P1 ] = −As , dη taking into account the correct upstream and downstream saturation values. Now for the temperature front, pressure at η = ηˆTk where k = + or k = − depending on the sided limit of concern can be expanded as, P (ˆ ηTk ) ∼ P0 (0k ) + ε AT  dP0 dη  + P1 (0k ) exp (iny + σt) η=0k  as was the case for the saturation front. The jump in pressure across the temperature front must be zero and so we get a solvability condition for P1 , [P1 ] = −AT  dP0 . dη  Now the thermal flux in the travelling wave coordinate frame is given by, ˆ ) T. (qs + qo − x  130  3.3. 1D Temperature Model Physically the jump in normal thermal flux must be given by, [ˆ x · (qs + qo ) − T ] = total flux across η = ηˆT+ − total flux across η = ηˆT− ∼ total flux across η = 0+ − total flux across η = 0− = (s(0+ ) + o(0+ ))T (0+ )  d ηˆT dt  − (s(0− ) + o(0− ))T (0− ) η=0+  with, d ηˆT = εσAT exp (iny + σt) dt being the speed of the temperature interface. However, we also have that, ¯ [ˆ x · (qs + qo ) − T ] = −Λ  dP0 ¯ dP1 T0 exp (iny + σt) T0 − [T0 ] + ε −Λ dη dη  where we have explicitly used that, s1 = T1 =  dT0 = 0. dη  Combining the two expressions we get that, ¯ −Λ  dP0 T0 = [T0 ] dη  which is indeed satisfied and more importantly that, ¯ dP1 T0 = AT σ. Λ dη  131  d ηˆT dt  η=0−  3.3. 1D Temperature Model We then have a complete set of solvability conditions in order to have nontrivial solutions, ¯ dP1 = 0, Λ dη dP1 α¯s = (s− − C)As σ, dη ¯ dP1 T0 = AT σ, Λ dη  (3.51a) (3.51b) (3.51c) dP0 , dη  [P1 ] = −(As + AT )  (3.51d)  which are the jump conditions for this thermal problem. We can write (3.51) in a matrix problem as,          −  (s − 1 ¯ Λ(C)  −  C) nσ 1 ¯ −) Λ(s  −  α¯s (s ) α¯s (C)  0 1 ¯ Λ(C)  1  −1  −  0  ¯ −) Λ(s  ¯ Λ(C)  0  0  ¯ −) Λ(s  0  σ n    0 A  s        C−  0   =        C+  0     0 AT   1 ¯ −) Λ(s    where, like before, the determinant of the matrix must vanish to avoid trivial solutions. Setting the determinant to zero leads to finding the roots of a problem quadratic in the variable σ/n, 1 1 σ ¯ − ¯ −) − (s − C) Λ(C) − Λ(s − ¯ ¯ n Λ(C) Λ(s ) σ σ ¯ 1 1 ¯ −) − − (s− − C) Λ(C) + Λ(s − ¯ − ¯ n n Λ(C) Λ(s ) σ ¯ ¯ −) = 0 − α¯s (s− )Λ(C) − α¯s (C)Λ(s n  132  3.3. 1D Temperature Model which we can solve to get the solutions, σ=0 and ¯ − ) − Λ(C) ¯ σ Λ(s = 2¯ − . ¯ n Λ(s ) + Λ(C) When δ  (3.52)  1 then the leading order term in the numerator is (s− )3 − C 3  which is positive since s− > C and thus we expect that fingering will occur. If the speed of the temperature front is some general v then the eigenvalue generalizes to,  ¯ − ) − Λ(C) ¯ Λ(s σ . = 2v ¯ − ¯ n Λ(s ) + Λ(C)  Notice that this problem qualitatively behaves differently from the classic saturation shock. For all values of n there is a zero eigenvalue. Notice if we choose this σ = 0 eigenvalue then we get from the jump conditions (3.51) that there will only be a solution if, As = −AT i.e. if the jump in saturation is balanced by an equal but opposite jump in the temperature. This corresponds to the the two fronts separating away from each other and thus their individual instabilities will cause destructive interference into a neutrally stable front. Conversely if As = −AT then σ = 0 is the only possible eigenvalue. We do not anticipate that this solution is a physical one because the jump in temperature should immediately cause 133  3.3. 1D Temperature Model a jump in saturation which is only explained by a non-degenerate front. Therefore, we expect that the other eigenvalue is associated with a nondegenerate front. Indeed, if we substitute the expression for σ into the matrix we get that, As = AT . The non-zero eigenvalue then is associated with constructive interference; each front goes unstable in the classical shock sense and then collaborates to create an even larger fingering instability with twice the growth rate. The evidence of thermo-viscous effects are noted experimentally in [19] where they inject hot glycerin into cold glycerin. They choose glycerin due to the very thermal sensitive viscosity of the substance. They find that fingering occurs right away at the thermal interface which provides evidence that viscous fingering can indeed be due to sharp changes in temperature. Figure 3.18 is from [19] and it shows the thermal fingering.  134  3.3. 1D Temperature Model  Figure 3.18: Figure from [19] showing thermo-viscous fingering in glycerin. The snapshots are of different instances in time. Note that this is a single fluid and the media is not porous so we do not expect the same physics in that regard. In [37], they perform experiments in porous media by injecting water into glycerin, thus it is the two-phase problem that we have analyzed. They perform experiments in an isothermal setting and then by injecting water at a temperature 10K above that of the glycerin. They arrive at two interesting conclusions at the end of their study, firstly they say that the inclusion of the temperature gradient has no visible qualitative effects on finger pattern or finger growth. This result does not hinder our findings, in fact for such a low temperature gradient, secondary fingering patterns should be nearly negligible. In our model, as the temperature difference goes to zero we retrieve the isothermal case which 135  3.3. 1D Temperature Model has just a single fingering pattern. Indeed, Figure 3.19 shows the eigenvalue at the temperature front for decreasing thermal gradients. Linear Eigenvalue Growth Rate for Various Thermal Gradients 7  T=1 T=1/3 T=1/30 T=1/300  6  5    4  3  2  1  0  0  0.5  1  1.5  2  2.5 n  3  3.5  4  4.5  5  Figure 3.19: Fingering growth rate for various values of non-dimensional temperature gradient, ∆T . The solid line is ∆T = 1, the dashed line is ∆T = 1/3, the dot-dashed line is ∆T = 1/30 and the dotted line is ∆T = 1/300. The second conclusion that [37] draws is that despite the lack of qualitative difference in fingers, they notice that the single fingering pattern (the classical shock) travels at a slower transverse speed in the presence of a thermal gradient versus the isothermal case. Recall from our parameters that in the presence of a thermal gradient, our downstream classical shock had the velocity, v = 6.30 136  3.3. 1D Temperature Model when we used the flux function in the cold region as, f2 =  s3 . s3 + δ exp (−1) (1 − s)3  If we were to compute an isothermal calculation using the same flux function, i.e. finding the shock saturation value s∗ that satisfies, f2 (s∗ ) f2 (s ) = s∗ ∗  then we would get that the isothermal velocity vo is, vo = f2 (s∗ ) = 6.63. Thus we see that in our model, the inclusion of the upstream temperature shock does indeed decrease the travelling speed of the downstream classic shock. By comparing Figure 3.17 to Figure 3.19 we see that for the temperature gradient we have in our problem (∆T = 1), the classical shock admits larger eigenvalues than the temperature shock. This is not unsurprising since the classic shock moves faster than the temperature shock, and shock velocity is in the numerator of the eigenvalue expression. In general, for classic shock speeds vc and temperature shock speeds vT , the classic shock will admit larger eigenvalues provided that, 2 ¯ ¯ − ) − Λ(0)) ¯ ¯ − )Λ(0) ¯ ¯ (2vT − vc )Λ(C)( Λ(s + (2vT + vc )(Λ(s − Λ(C) ) < 0.  137  Chapter 4 Conclusions and Future Work We have presented a mathematical analysis of a problem as it relates to the remediation of contaminants from soil. We began by developing a model using local averaging techniques that explained the conservation of mass of water and oil if only convective mass transport was considered. To describe the fluid flux, we used Darcy’s law which states that in a porous media, mass transport velocities are due to pressure differences. In the end we had a set of partial differential equations describing water saturation in the soil. We then coupled this with a thermal model that took into account, the transfer of energy through conduction and convection based on a fixed heat injection into the soil. We concluded the model by using parameters that apply directly to the remediation problem which led to the advent of two small numbers. Firstly we saw that the viscosity ratio between water and oil (or benzene in our case) was quite small and that the inverse capillary number, an effective measurement of surface tension over the domain length scale was small as well. In fact, for larger contamination sites, this becomes even smaller. However, as our analysis demonstrated, the capillarity is still quite significant over short length scales in the model.  From the model we proceeded by considering a one-dimensional system of  138  Chapter 4. Conclusions and Future Work oil and water separation. We neglected the small capillary effects entirely and were led to the famous Buckley-Leverett formulation. This solution admits a shock emerging from a rarefaction wave which connects an upstream shock saturation value, s∗ to a downstream water absent (s = 0) case. We were able to asymptotically determine, based on the small viscosity ratio, an analytic form for the shock saturation value. We then computed a linear stability analysis for this case, noticing that any positive eigenvalues would be associated with unstable modes and produce viscous fingers. The result we obtained was the analogue to the Taylor-Saffman instability ([32]) but with a non-uniform upstream permeability caused by the rarefaction wave. The big difference between our result and the Taylor-Saffman instability was that the Taylor-Saffman instability admits viscous fingering as long as the displacing fluid is less viscous than the displaced one (which is indeed the case for water displacing oil). However, based on our asymptotic result for the upstream saturation value, we showed that instabilities persisted with sufficiently small viscosity ratio. Although, we demonstrated an example where the displacing fluid was less viscous than the displaced one (satisfying the Taylor-Saffman instability) but that only stable eigenvalues formed. This means that as long as the viscosity ratio is reasonably close to one but still less than one, then the Buckley-Leverett formulation allows for a less viscous fluid to displace one of higher viscosity stably.  With our completion of the Buckley-Leverett analysis we then incorporated the surface tension effects via a small capillary pressure between the oil and water. This formulation led to an outer problem where the capillary pres-  139  Chapter 4. Conclusions and Future Work sure was negligible and an inner problem where it was a dominating factor. The outer problem corresponded to the steady-state of the Buckley-Leverett formulation and thus appeared as a hyperbolic shock. However, the inner problem was a travelling wave solution that connected the upstream and downstream shock values. We computed a linear stability analysis on the wave front which also produced an inner and outer problem. The outer problem, like with the base state, provided no new information that was not already provided by the Buckley-Leverett analysis. However, the same methodology could not be used to determine fingering because the outer problem was only a limiting case of the region where capillary pressure dominated. In this inner region we numerically solved the eigenvalue problem and saw that the resulting eigenvalue differed from the Buckley-Leverett linear eigenvalue by having a maximal wavenumber and a finite range of instability. This means that viscous fingering will have a favoured wave mode and that there will be a wave number, nc , beyond which no fingers will form. Asymptotically we showed that to first order, the linear eigenvalue still dominated as expected so that in the limit of vanishing capillary effect, the BuckleyLeverett analysis holds. The reason for the critical wave number at which fingers stop forming is due to the competing nature of finger instability and capillary effects. The capillary number measures the surface tension effects on the porous level between water and oil. As the wave number increases, the wavelength of the fingers decreases and thus at some n, this capillary effect will dominate and hinder the growth of the fingers. With a large n scaling we saw the formation of an -invariant system which still maintained the maximal growth and cut-off wave number behaviour.  140  Chapter 4. Conclusions and Future Work  We next considered the inclusion of thermal effects in the model. To do this, we maintained a one-dimensional setting but released the requirement of capillary effects for both simplicity of the model and due to the fact that the inverse capillary number is small. We also neglected thermal diffusion on the basis that parameter choices showed it to be an irrelevant mechanism for heat transfer. A more controversial decision was to simplify the convection terms by equating the convection coefficients to one. In essence, this has the same effect as making the volumetric heat capacities, ρcp equal to one another. This was a reasonable assumption since each differed from unity by an order one function, precisely the type that multiplied each convection number, but the more pragmatic reason was that it simplified analysis of the model. With our choice of making the convection numbers on par and equal to unity, we decoupled saturation from the energy equation and thus were able to produce a trivial shock wave solution for the temperature front. We then incorporated this result into the temperature dependent viscosities and set to find solutions for saturation. We developed theory to show the different possibilities that could occur with the sharp temperature front and argued that the saturation went from an initial value through a rarefaction wave to the shock front associated with the drop in temperature at which point what followed was a constant saturation value that eventually formed another shock solution to s = 0. We showed that while each of these fronts of course were steady states, they were only steady-states in a frame of reference centered around their travelling velocities. Since the downstream, classical shock had a larger velocity then the two fronts would separate in time. We  141  Chapter 4. Conclusions and Future Work performed a linear stability analysis around each front solution. For the downstream classical shock we concluded that the eigenvalue was linear just as in the complete isothermal case but corrected for the appropriate viscosity of the temperature in that region. This is to be expected since the analysis of this front undergoes no temperature changes and this is quasi-isothermal. The reason it is not purely isothermal is because the upstream effects from temperature change the downstream shock saturation value for the classical shock from what it would have been in a true isothermal setting. Upstream, along the temperature shock we concluded that the eigenvalue was also linear but was double the value that would have been generated in an isothermal setting. This is due to the fact that when an instability forms, the saturation front fingers but because the temperature varies across this front, the thermal front fingers as well. As each of these fronts grow, they collaborate into constructive wave interference, thus doubling the growth of the visible front, hence accounting for this factor of 2. We postulated in the introduction that the inclusion of thermal effects would promote fingering and indeed, we see that we not only have induced thermal fingering across the temperature shock but as evidenced in Figure 3.17, the downstream classic shock produces larger eigenvalues in the presence of an upstream temperature change. We also postulated that the promotion of fingering would support soil remediation. While we cannot make any astute conclusions on that hypothesis without understanding the dynamics of the fingers post formation, it is worth noting that there is a competitive behaviour between the two shocks. As noted, the downstream shock travels faster than the temperature shock, giving rise to the idea that the fingers would spread out over time. However,  142  Chapter 4. Conclusions and Future Work along the temperature shock, the instability has doubled the growth rate and thus smaller wavelengths can produce larger instabilities, promoting higher amplitude fingers if mass conservation is considered. Therefore we have this competing nature between narrow longer fingers forming upstream trying to interact with fingers that are downstream moving at a faster speed. If there is a point where these fingers intersect then one could imagine that it would lead to an overall stabilizing front, pushing the oil towards the collection site.  Overall, we have attempted to analyze a thermo-saturation model with specific application to the parameters associated with remediation. We have demonstrated that the inclusion of surface tension between the oil and water, even if nearly negligible, has an immense effect on the stability properties of the associated linear eigenvalue. This statement holds true for any problem where the inverse capillary number can be assumed small. For instance, if the two substances are only weakly immiscible then it is reasonable to assume that there will be very little surface tension effect and hence there will be a small surface tension coefficient leading to a small inverse capillary number. This analysis shows that even for vanishing immiscibility, there is a drastic difference in the behaviour of fingers produced. We have also shown that the inclusion of thermal energy, even in a simple way, has altered the stability analysis by creating a dual finger profile, one upstream where the temperature changes and one down-stream that is associated in the classical isothermal sense. We see that even without considering the thermal fingering, the change in temperature upstream, effects the downstream shock speed and linear eigenvalue. This means that even for small thermal gradients where  143  4.1. Future Work upstream temperature shocks may be nearly negligible, there will still be an effect on the classical downstream fingering.  4.1  Future Work  In consideration of what we have presented, there are many mathematical avenues to explore both for completeness of the problem and interest to the application. Firstly, with the capillary pressure flow, for simplification we considered J(s), the Leverett function, to be chosen such that the saturation problem was non-degenerate. While we are not the first people to use such an assumption ([34]), by removing the restriction we would allow for a complete solution with the function fitted to real soils. We demonstrated that with degeneracy, the solution, u, had the behaviour, u ∼ u1/3 ,  u→0  leading to saturation profiles that decay to the base-state in finite space while having a divergent derivative at u = 0. The major interests with this problem are related to the numerical computation of the base-state. Since the derivative diverges, most numerical codes will not be able to resolve the sharp front that occurs. We suspect that there are two methods to resolving the sharp degenerate front. Firstly, knowing the form of the solution as it degenerates, allows us to perform an analytic-numeric matching technique. We compute the base-state solution numerically as before up to a certain tolerance where the numeric scheme has trouble resolving the near divergent derivative. We then insert the analytic solution using the terminus point 144  4.1. Future Work of the numerical simulation as an initial condition which transitions the numeric solution through its region of degeneracy. The second solution is to use a modified numeric scheme, similar to the way in which we handled the hyperbolic problem, a method that relaxes computation near smooth points and intensifies it near the sharp pseudo-hyperbolic front associated with the degeneracy. Once the base-state is computed we predict that the linear stability analysis will remain mostly unaffected, except that the inner solution domain will have a finite terminus point, the point where the degenerate front satisfies u = 0. However, because hyperbolicity is conserved near the degenerate point u = 0, then there will be a shock-like front there. We expect that we will require the appropriate conservation arguments across this front to determine how it is effected by linear perturbations  Secondly, while we have introduced thermal injection into the soil via heated water, we wish to be more specific about the heating mechanism. Most remediation companies heat their water via electrodes in the ground [27]. We therefore wish to introduce into the thermal equation, a source term that represents the heat generated inside the soil via the electric field generated by the electrodes. Of course, like with the viscosities, both water and oil will have different electrical conductivities and thus we must consider the ratio associated with that. We predict that with the inclusion of electric effects, fingering will once again promote remediation. It is anticipated that the fluid displacement will eventually lead to the most optimal path for the current from the electrode thus promoting heating and reducing the viscosity ratio between the fluids ultimately driving the oil towards the collection site. Also  145  4.1. Future Work in consideration of thermal effects, it may be worth including the boundary layer that arises from the small diffusion term in the temperature equation. While we demonstrated that diffusion was infinitesimally small, the boundary layer allows the temperature to be continuous. We postulate that the inclusion of such effects will mimic the analysis of small capillary pressure in that there will be a limit on the maximal fingering growth and a cut-off wave number where fingers no longer form. If this is true then it will be interesting to determine how these maximal and cut-off wave numbers relate to the downstream finger-growth, specifically if there will be wave numbers such that one type of finger development stabilizes while the other is unstable.  Finally we would like to perform numerical simulations for viscous fingering to confirm or refute many of the postulations we have regarding efficiency of soil remediation. With proper simulations we could also determine finger width which relates to the transverse wave number n and thus we could compare the simulated fingers with the predicted wave number of maximum growth. Recently in [16], they simulate liquid snowflake interfacial fingering via numerical level set methods which are interface tracking. Due to the fingering dependency on capillarity, there is no strict interface or shock to capture and thus these methods are not satisfactory to our problem. Instead, because there will be a localized smooth front, we expect that we will require methods that are adaptive in space and implicit in time. We will then use these methods with an initial front assumption composed of various wave modes and observe how the numerical computations resolve the front’s longtime behaviour.  146  4.1. Future Work  While there are many other mathematical considerations to develop, this set encompasses the most logical future goals from a pragmatic framework, with the end aim of contributing to a stronger understanding of the underlying physics involved in the oil extraction process, something that is of paramount concern as evidenced by [1].  147  Bibliography [1] In Situ Oil Sands Summit 2011. http://www.in-situ-oil-sands. com/, May 2011. [2] M. Assael, M. Papadaki, and W. Wakeham. Measurements of the viscosity of benzene, toluene, and m-xylene at pressure up to 80 mpa. International Journal of Thermophysics, 12:449–457, 1991. [3] J. Baba and P. D. Komar. Measurements and analysis of setting velocities of natural quartz sand grains. Journal of Sedimentary Research, 51(2):631–640, 1981. [4] S.A. Bowers and R.J. Hanks. Specific heat capacity of soils and minerals as determined with a radiation calorimeter. Soil Science, 94(6):392–396, 1962. [5] L. Bridge, R. Bradean, M.J. Ward, and Wetton B.R. The analysis of a two-phase zone with condensation in a porous medium. Journal of Engineering Mathematics, 45:247–268, 2003. [6] S.E. Buckley and M.C. Leverett. Mechanism of fluid displacement in sands. Transactions of the AIME, 146:107–116, 1941.  148  Bibliography [7] R.G. Carbonell and S. Whitaker. Heat and mass transfer in porous media. In Fundamentals of Transport Phenomena in Porous Media. 1984. [8] M.Y. Corapcioglu, K. Tuncay, and B.K. Ceylan. Oil mound spreading and migration with ambient groundwater flow in coarse porous media. Water Resources Research, 32(5):1299–1308. [9] A.H. Craven and Peletier L.A. Similarity solutions for degenerate quasilinear parabolic equations. Journal of Mathematical Analysis and Applications, 38:73–81, 1972. [10] M.B. Dusseault. Comparing Venezuelan and Canadian Heavy Oil and Tar Sands. Technical report, Canadian International Petroleum Conference, Calgary, Canada. [11] G. Elert. Viscosity. http://http://physics.info/viscosity/. Accessed August 17, 2011. [12] R. Haberman. Applied Partial Differential Equations. Pearson, 4th edition, 2004. [13] T. Harter. Basic Concepts of Groundwater Hydrology. University of California, 2003. http://groundwater.ucdavis.edu/Publications/ Harter_FWQFS_8083.pdf. [14] W.M Haynes, editor. Standard Density of Water. CRC Handbook of Chemistry and Physics, 91st Edition (Internet Version 2011). CRC Press/Taylor and Francis. 149  Bibliography [15] W.M Haynes, editor. Thermophysical Properties of Water and Steam. CRC Handbook of Chemistry and Physics, 91st Edition (Internet Version 2011). CRC Press/Taylor and Francis. [16] M.G. Hennessy. Liquid snowflake formation in superheated ice. Master’s thesis, University of Oxford, 2010. [17] R.N. Hills and P.H. Roberts. A note on the kinetic conditions at a supercooled interface. Int. Comm. Heat Mass Transfer, 20:407–416, 1993. [18] E.J. Hinch. Perturbation Methods. Cambridge Texts in Applied Mathematics. 1991. [19] Kristi E Holloway and John R de Bruyn. Viscous fingering with a single fluid. Canadian Journal of Physics, 83(5):551–564, 2005. [20] G.M. Homsy. Viscous fingering in porous media. Annual Review Fluid Mechanics, 19:271–311, 1987. [21] M. Kaviany. Principles of Heat Transfer in Porous Media. Mechanical Engineering Series. Springer, 2nd edition, 1995. [22] H. Kim and D.J. Burgess. Prediction of interfacial tension between oil mixtures and water.  Journal of Colloid and Interface Science,  241(2):509–513, September 2001. [23] L. Korson, W. Drost-Hansen, and F.J. Millero. Viscosity of water at various temperatures. J. Phys. Chem., 73(1):34–39, January 1969.  150  Bibliography [24] A.A. Lacey, J.R. Ockendon, and A.B. Tayler. “waiting-time” solutions of a nonlinear differential equation. SIAM Journal on Applied Mathematics, 42(6), December. [25] M.C. Leverett. Capillary behaviour in porous solids. Transactions of the AIME, 142:159–172, 1941. [26] S. F. Y. Li, G. C. Maitland, and W. A. Wakeham. Thermal conductivity of benzene and cyclohexane in the temperature range 3690c at pressures up to 0.33 gpa. International Journal of Thermophysics, 5:351–365, 1984. [27] C. Ma and J. Kingscott. Recent Developments for In Situ Treatment of Metal Contaminated Soils. Technical Report 68-W5-0055, U.S. Environment Protection Agency, Washington D.C., March 1997. [28] G. I. Makhatadze and P. L. Privalov. Partial specific heat capacity of benzene and of toluene in aqueous solution determined calorimetrically for a broad temperature range. The Journal of Chemical Thermodynamics, 20(4):405 – 412, 1988. [29] T. Nagai and M. Mimura. Asymptotic behaviour for a nonlinear diffusion equation in population dynamics. SIAM J. Appl. Math., 43(3):449–464, 1983. [30] I. Nazad, R.G. Carbonell, and S. Whitaker. Heat condition in multiphase systems i: Theory and experiments for two-phase systems. Chem. Engng. Sci., 40:843–855, 1985.  151  Bibliography [31] J. Ockendon, S.D. Howison, and A.B. Movchan. Applied Partial Differential Equations. Oxford University Press, 2nd edition, 2003. [32] Saffman P.G. and Taylor G.I. Viscous Fingering in Hele-Shaw Cells. Proceedings of the Royal Society of London, Series A, 245, 1958. [33] M.L.V. Ramires et al. Standard reference data for the thermal conductivity of water. Journal of Physical and Chemical Reference Data, 24(3):1377–1381, 1995. [34] A. Riaz and H.A. Tchelepi. Linear stability analysis of immiscible twophase flow in porous media with capillary dispersion and density variation. Phys. Fluids, 16(12):4727–4737, December 2004. [35] Amir Riaz and Hamdi A. Tchelepi. Numerical simulation of immiscible two-phase flow in porous media. 18(1):014104, 2006. [36] M.E. Ritter. The Physical Environment: an Introduction to Physical Geography.  http://www.uwsp.edu/geo/faculty/ritter/geog101/  textbook/title_page.html. Accessed May 10, 2011. [37] M. Z. Saghir, O. Chaalal, and M. R. Islam. Numerical and experimental modeling of viscous fingering during liquid-liquid miscible displacement. Journal of Petroleum Science and Engineering, 26(1-4):253 – 262, 2000. [38] J.S. Selker, T.S. Steenhuis, and J.Y. Parlange. Wetting front instability in homogeneous sandy soils under continuous infiltration. Soil Science Society of America Journal, 56(5):1346–1350, SEP-OCT 1992. 56th Annual Meeting of the Soil Science Soc of America, Minneapolis, MN, Nov 01-06, 1992. 152  Bibliography [39] Shanghai Niumag Corporation. Porous media. http://www.micromr. com/niumag/micromr\%20features.htm. Accessed July 28, 2011. [40] J.C. Slattery. Momentum, Energy, and Mass Transfer in Continua. R.F. Krieger, 2nd edition, 1981. [41] J.C. Slattery. Single-phase flow through porous media. AIChE J., 15:866–872, 1989. [42] J. Smoller. Shock Wave and Reaction-Diffusion Equations. A Series of Comprehensive Studies in Mathematics. Springer-Verlag, second edition, 1994. [43] M. Van Dyke. Perturbation Methods in Fluid Mechanics. The Parabolic Press, 1975. [44] W.F. Waite et al. Thermal conductivity measurements in porous mixtures of methane hydrate and quartz sand. Geophysical Research Letters, 29(24):82(4), 2002. [45] K. Woodruff and D. Miller. Newtown Creek/Greenpoint oil spill study. Technical report, Lockheed Martin/REAC, September 2007. [46] Y.C. Yortsos and F.J. Hickernell. Linear-stability of immiscible displacement in porous media. SIAM J. Appl. Math., 49(3):730–748, June 1989. [47] Y.C. Yortsos and A.B. Huang. Linear-stability analysis of immiscible displacement:part 1-simple basic flows. SPE Reservoir Engineering, 1:378–390, 1986. 153  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items