"Science, Faculty of"@en . "Mathematics, Department of"@en . "DSpace"@en . "UBCV"@en . "Moyles, Iain"@en . "2011-08-26T18:29:56Z"@en . "2011"@en . "Master of Science - MSc"@en . "University of British Columbia"@en . "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 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."@en . "https://circle.library.ubc.ca/rest/handle/2429/36926?expand=metadata"@en . "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\u00C2\u00A9 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 analy- sis 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 in- terface 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 im- portantly, 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 ther- mal gradients across the interface where temperature drops induce unstable viscous patterns with a higher wave number than would occur for an equiv- alent 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 . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 Volume Averaging Method . . . . . . . . . . . . . . . . . . . 7 2.2 Water Saturation . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2.1 Capillary Pressure . . . . . . . . . . . . . . . . . . . . 15 2.2.2 Non-Dimensionalization . . . . . . . . . . . . . . . . . 19 2.2.3 Cory-Type Permeabilities . . . . . . . . . . . . . . . . 21 2.3 Temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.3.1 Volume Averaged Temperature Equation . . . . . . . 23 2.3.2 Temperature Dependent Viscosity . . . . . . . . . . . 36 2.3.3 Non-Dimensionalization . . . . . . . . . . . . . . . . . 37 iv Table of Contents 3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.1 1D Model in the Absence of Capillary Pressure . . . . . . . . 42 3.1.1 Linear Stability Analysis . . . . . . . . . . . . . . . . 50 3.2 1D Model With Capillary Pressure . . . . . . . . . . . . . . . 61 3.2.1 Linear Stability Analysis . . . . . . . . . . . . . . . . 70 3.2.2 Numerical Solution . . . . . . . . . . . . . . . . . . . . 75 3.2.3 Asymptotic Analysis . . . . . . . . . . . . . . . . . . . 83 3.3 1D Temperature Model . . . . . . . . . . . . . . . . . . . . . 99 3.3.1 Stationary Boundary Characteristics . . . . . . . . . . 106 3.3.2 Internal Interface Characteristics . . . . . . . . . . . . 107 3.3.3 Other Cases . . . . . . . . . . . . . . . . . . . . . . . 116 3.3.4 Solving The Saturation Problem . . . . . . . . . . . . 117 3.3.5 Linear Stability Analysis . . . . . . . . . . . . . . . . 121 4 Conclusions and Future Work . . . . . . . . . . . . . . . . . . 138 4.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 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]. 3 2.1 Here we consider a control volume V which is broken up into its volume fractions. The porosity \u00CF\u0086 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, \u00CE\u0093, 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\u00E2\u0088\u0097 = 0.1443 and t = 1/8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.3 Saturation steady-state solution for the inner boundary layer centered around \u00CE\u00B7 = 0. It is a smooth heteroclinic orbit con- necting the steady-state values of the outer hyperbolic problem. 68 3.4 Pressure derivative steady-state solution for the inner bound- ary layer centered around \u00CE\u00B7 = 0. The derivative has the ap- propriate 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 discretiza- tion 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 in- formation flows to describe the i+ 1/2 cell. . . . . . . . . . . . 76 3.6 Numerical solution of u\u00CC\u00821 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 \u00CF\u0083 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 or- der (dash-line) for n = 1. 201 mesh points were used in the solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 3.10 Pressure eigenfunction computed numerically (solid line), asymp- totically 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\u00CC\u0082 = \u000Fn. . . . . . . . . . . 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 char- acteristics, dashed lines represent the constant solution char- acteristics 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 = \u00CE\u00BE corresponds to the traditional Buckley- Leverett 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 dot- dashed 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 charac- teristics from the rarefaction wave, the dashed lines represent the constant solution characteristics from the internal inter- face 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 = \u00CE\u00BE 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, \u00E2\u0088\u0086T . The solid line is \u00E2\u0088\u0086T = 1, the dashed line is \u00E2\u0088\u0086T = 1/3, the dot-dashed line is \u00E2\u0088\u0086T = 1/30 and the dotted line is \u00E2\u0088\u0086T = 1/300. . . . . . . . . . . . . . . . 136 x Acknowledgements Firstly, I\u00E2\u0080\u0099d 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 mathemat- ics, and as a trusted representative for academic-industrial relations, promot- ing 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\u00E2\u0080\u0099d 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 coun- try 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 associ- ated with every aspect of this research and my life. I thank her profusely for her love and belief in me. Finally I\u00E2\u0080\u0099d like to thank many of the friends who have provided welcome distractions through various points of the research process. I\u00E2\u0080\u0099d 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\u00E2\u0080\u0099d like to thank Matt Hennessy for providing both a strong pride and admiration in his aca- demic endeavours that inspires me to be successful, for reading and editing this thesis, and also for the much needed vacation that allowed me to maxi- mize 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 pro- vided to me is greatly appreciated. xii Chapter 1 Introduction The process of soil remediation is designed, through various means, to re- move 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 man- dated 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 be- coming 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 compa- nies 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 mo- tivated 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 or- der 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 varia- tions 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 experimen- tally ([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 to- wards 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. There- fore, 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 satura- tion 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 remedia- tion. 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 pa- rameter 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 under- stand 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 previ- ous 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 eigen- value 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 in- jection 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 fin- gering 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 improve- ments 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 multi- component 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, \u00E3\u0080\u0088w\u00E3\u0080\u0089 = 1 V \u00E2\u0088\u00AB Vi w dV. 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, \u00E3\u0080\u0088w\u00E3\u0080\u0089\u00E2\u0088\u0097 = 1 Vi \u00E2\u0088\u00AB Vi w dV. The volume averaging method is essentially a formulation of homogeniza- tion 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 pe- riodic 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 informa- tion regarding the global problem. Before we derive the effective temperature equation we will require an appropriate divergence theorem for averaged vari- ables. It is not difficult to show ([21],[41],[40]) that the following identities hold, \u00E3\u0080\u0088\u00E2\u0088\u0087w\u00E3\u0080\u0089 = \u00E2\u0088\u0087\u00E3\u0080\u0088w\u00E3\u0080\u0089+ 1 V \u00E2\u0088\u00AB Aint wn dA, (2.1a) \u00E3\u0080\u0088\u00E2\u0088\u0087 \u00C2\u00B7w\u00E3\u0080\u0089 = \u00E2\u0088\u0087 \u00C2\u00B7 \u00E3\u0080\u0088w\u00E3\u0080\u0089+ 1 V \u00E2\u0088\u00AB Aint w \u00C2\u00B7 n dA (2.1b) 8 2.2. Water Saturation where n is the standard outward normal vector and Aint is the area of interfa- cial contact surfaces between adjacent subdomains. It will also be important that [21], \u00E3\u0080\u0088\u00E3\u0080\u0088w\u00E3\u0080\u0089\u00E3\u0080\u0089 = \u00E3\u0080\u0088w\u00E3\u0080\u0089 \u00E3\u0080\u0088\u00E2\u0088\u0087\u00E3\u0080\u0088w\u00E3\u0080\u0089\u00E3\u0080\u0089 = \u00E2\u0088\u0087\u00E3\u0080\u0088w\u00E3\u0080\u0089 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 \u00CF\u0086 { s o } (1\u00E2\u0088\u0092 \u00CF\u0086) Figure 2.1: Here we consider a control volume V which is broken up into its volume fractions. The porosity \u00CF\u0086 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) = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B31, x is in fluid0, 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 \u00CF\u0086 = 1 V \u00E2\u0088\u00AB V a(x) dV = Vf V where Vf is the fluid volume made up of oil and water. This volume fraction \u00CF\u0086, is volume available for flow while the remaining portion (1\u00E2\u0088\u0092\u00CF\u0086) is the solid matrix or non-porous portion. We must further divide \u00CF\u0086 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 \u00E2\u0089\u00A4 s \u00E2\u0089\u00A4 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, \u00E2\u0088\u0082\u00CF\u0081s \u00E2\u0088\u0082t +\u00E2\u0088\u0087 \u00C2\u00B7 (\u00CF\u0081sqs) = 0 where \u00CF\u0081 is the density, q is the volumetric flux (velocity) and the subscript 1More 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\u00E2\u0080\u0099s association with the water phase. We now take a local volume average of this equation, \u00E2\u0088\u0082 \u00E2\u0088\u0082t \u00EF\u00A3\u00AB\u00EF\u00A3\u00AD 1 V \u00E2\u0088\u00AB Vs \u00CF\u0081s dV \u00EF\u00A3\u00B6\u00EF\u00A3\u00B8+ 1 V \u00E2\u0088\u00AB Vs \u00E2\u0088\u0087 \u00C2\u00B7 (\u00CF\u0081sqs) dV = 0. 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 \u00E2\u0088\u00AB V \u00E2\u0088\u0082\u00CF\u0081s \u00E2\u0088\u0082t dV and by the previous statement that V is constant in time we can remove the derivative from the integrand, \u00E2\u0088\u0082 \u00E2\u0088\u0082t \u00EF\u00A3\u00AB\u00EF\u00A3\u00AD 1 V \u00E2\u0088\u00AB V \u00CF\u0081s dV \u00EF\u00A3\u00B6\u00EF\u00A3\u00B8 . We complete the argument by noting that, \u00CF\u0081s = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3\u00CF\u0081s x \u00E2\u0088\u0088 Vs0 x /\u00E2\u0088\u0088 Vs . We may now use the divergence theorem for averages (2.1b) to simplify the 12 2.2. Water Saturation continuity equation further, \u00CF\u0086 \u00E2\u0088\u0082s \u00E2\u0088\u0082t +\u00E2\u0088\u0087 \u00C2\u00B7 \u00E3\u0080\u0088qs\u00E3\u0080\u0089+ 1 V \u00E2\u0088\u00AB Aint qs \u00C2\u00B7 n dA = 0 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, \u00CF\u0086s \u00E2\u0089\u00A1 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 im- posing the boundary condition that the normal component of the velocity along the interfacial boundaries vanishes, nij \u00C2\u00B7 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, \u00CF\u0086 \u00E2\u0088\u0082s \u00E2\u0088\u0082t +\u00E2\u0088\u0087 \u00C2\u00B7 qs = 0 (2.3) 13 2.2. Water Saturation where we have removed the average brackets for simplicity and assume here- after 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, \u00E2\u0088\u0087 \u00C2\u00B7 (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\u00E2\u0080\u0099s Law [21], qs = \u00E2\u0088\u0092\u00CE\u00BAs \u00C2\u00B5s \u00E2\u0088\u0087Ps, (2.5) qo = \u00E2\u0088\u0092\u00CE\u00BAo \u00C2\u00B5o \u00E2\u0088\u0087Po (2.6) where q is the water flux like before, \u00CE\u00BA is the permeability, \u00C2\u00B5 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, Pc = Po \u00E2\u0088\u0092 Ps. (2.7) We will also allow the permeability, \u00CE\u00BA 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 \u00CE\u0093 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 \u00CE\u0093 x = 0 x = L h(x) Figure 2.2: A 1D model interface contact, \u00CE\u0093, 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 \u00CE\u00B3ss be the energy per unit length of this surface. Secondly, there is the contact surface between oil and the substrate and \u00CE\u00B3os 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 \u00CE\u00B3 such 15 2.2. Water Saturation that the energy of the entire surface is2, E = \u00CE\u00B3 \u00E2\u0088\u00AB | d\u00CE\u0093|+ \u00CE\u00B3osL+ \u00CE\u00B3ssLss 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, | d\u00CE\u0093| = \u00E2\u0088\u009A 1 + [h\u00E2\u0080\u00B2(x)]2 dx 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, \u00E2\u0088\u00AB L 0 h(x) dx = A, 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, Edroplet = \u00CE\u00B3 \u00E2\u0088\u00AB L 0 \u00E2\u0088\u009A 1 + [h\u00E2\u0080\u00B2(x)]2 dx subject to the constraint. Let u = h(x) and v = h\u00E2\u0080\u00B2(x) and use the standard Euler-Lagrange equations to minimize this integral, i.e. find a solution to, d dx Fv = Fu 2We 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 = \u00CE\u00B3 \u00E2\u0088\u009A 1 + v2 + \u00CE\u00BBu with \u00CE\u00BB a Lagrange multiplier. Upon carrying out these operations we get that, \u00CE\u00B3 d dx v\u00E2\u0088\u009A 1 + v2 = \u00CE\u00BB. Now, we notice that the curve y = h(x) = u has normal vector, n\u00CC\u0082 = 1\u00E2\u0088\u009A 1 + v2 \u00E3\u0080\u0088v,\u00E2\u0088\u00921\u00E3\u0080\u0089 and so we can write our Euler-Lagrange expression as \u00CE\u00B3\u00E2\u0088\u0087 \u00C2\u00B7 n\u00CC\u0082 = \u00CE\u00BB, \u00CE\u00B3\u00CE\u00BD = \u00CE\u00BB where we define the curvature \u00CE\u00BD as the standard divergence of the normal vector, \u00CE\u00BD = \u00E2\u0088\u0087 \u00C2\u00B7 n\u00CC\u0082. This derivation can easily be extended to higher dimensions. In our case we have a 2-dimensional surface energy to minimize and the parameter \u00CE\u00B3 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, \u00CE\u00B3\u00CE\u00BD = Pc. Now Leverett in [25] determined an expression for the curvature in terms of the saturation s, \u00CE\u00BD = \u00E2\u0088\u009A \u00CF\u0086 \u00CE\u00BA 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 \u00CE\u00BA is due to the fact that this is the total permeability. Leverett\u00E2\u0080\u0099s motivation for creating such an expression was sim- ply 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 Lev- erett\u00E2\u0080\u0099s 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 = \u00CE\u00B3 \u00E2\u0088\u009A \u00CF\u0086 \u00CE\u00BA J(s). (2.8) 18 2.2. Water Saturation 2.2.2 Non-Dimensionalization Consider the rescaling, x\u00CC\u0084 = x L , t\u00CC\u0084 = t \u00CF\u0084 , q\u00CC\u0084i = qi Q , \u00CE\u00B1i = \u00CE\u00BAi \u00CE\u00BA , P\u00CC\u0084i = Pi \u00CE\u00BA \u00C2\u00B5sQL where L is a characteristic length, \u00CF\u0084 is a characteristic time, \u00CE\u00BA is the total permeability, \u00CE\u00B1i is the relative permeability of species i, and Q is a character- istic velocity. Using these non-dimensional variables, we get for our equations (after removing the overbars), \u00E2\u0088\u0087 \u00C2\u00B7 (qs + qo) = 0, (2.9a) \u00CF\u0086Fl \u00E2\u0088\u0082s \u00E2\u0088\u0082t +\u00E2\u0088\u0087 \u00C2\u00B7 qs = 0, (2.9b) qs = \u00E2\u0088\u0092\u00CE\u00B1s(s)\u00E2\u0088\u0087Ps, (2.9c) qo = \u00E2\u0088\u0092\u00CE\u00B4\u00CE\u00B1o(s)\u00E2\u0088\u0087Po, (2.9d) Po \u00E2\u0088\u0092 Ps = Ca\u00E2\u0088\u00921J(s), (2.9e) where here we introduce the non-dimensional numbers, \u00CE\u00B4 = \u00C2\u00B5s \u00C2\u00B5o , Viscosity ratio, (2.10a) Ca = \u00C2\u00B5sQL \u00CE\u00B3 \u00E2\u0088\u009A \u00CE\u00BA\u00CF\u0086 , Capillary Number, (2.10b) Fl = L Q\u00CF\u0084 , Flow Number. (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 \u00CF\u0086Fl = 1, i.e., Q = L\u00CF\u0086 \u00CF\u0084 . (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. Table 2.1: Parameter Values for Model. Symbol Value Description Reference \u00CF\u0086 0.49 Porosity [21] (average value) L 20m Characteristic Chosen length scale \u00CF\u0084 1 day Characteristic Chosen time scale Q 9.8m/day Characteristic velocity Eqn (2.11) \u00CE\u00BA 7.1\u00C3\u0097 10\u00E2\u0088\u009212 m2 Total permeability [21] (average value) \u00C2\u00B5s 1.00\u00C3\u0097 10\u00E2\u0088\u00923Pa\u00C2\u00B7s Water viscosity (300K) [23] \u00C2\u00B5o 0.523Pa\u00C2\u00B7s Oil viscosity(300K) Chosen \u00CE\u00B3 32.7mN/m Interfacial tension [22] of water & benzene \u00CE\u00B4 1/523 Viscosity Ratio Eqn (2.10a) Ca 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\u00E2\u0080\u0099s 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\u00C2\u00B7s ([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, \u000F = Ca\u00E2\u0088\u00921 = 0.03. We also note that by our choice of parameters, \u00CE\u00B4 \u001C 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, \u00CE\u00B1i(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, \u00CE\u00B1s(s) = s 3, (2.12a) \u00CE\u00B1o(s) = (1\u00E2\u0088\u0092 s)3, (2.12b) 21 2.3. Temperature which is always positive since 0 \u00E2\u0089\u00A4 s \u00E2\u0089\u00A4 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 \u00E2\u0080\u009Cpoint\u00E2\u0080\u009D heat equation which is a continuous model valid in each subdomain volume. We have, (\u00CF\u0081cp)r \u00E2\u0088\u0082Tr \u00E2\u0088\u0082t = \u00E2\u0088\u0092\u00E2\u0088\u0087 \u00C2\u00B7 JTr (\u00CF\u0081cp)s \u00E2\u0088\u0082Ts \u00E2\u0088\u0082t = \u00E2\u0088\u0092\u00E2\u0088\u0087 \u00C2\u00B7 JTs (\u00CF\u0081cp)o \u00E2\u0088\u0082To \u00E2\u0088\u0082t = \u00E2\u0088\u0092\u00E2\u0088\u0087 \u00C2\u00B7 JTo for rock, water, and oil respectively. We denote JT for thermal fluxes in the problem, cp as specific heat capacities and \u00CF\u0081 as the density. We assume that \u00CF\u0081cp is constant, hence its exclusion from the time derivative. We take the volume average of these equations to get, \u00E2\u008C\u00A9 (\u00CF\u0081cp)r \u00E2\u0088\u0082Tr \u00E2\u0088\u0082t \u00E2\u008C\u00AA = \u00E2\u0088\u0092\u00E3\u0080\u0088\u00E2\u0088\u0087 \u00C2\u00B7 JTr\u00E3\u0080\u0089 (2.13a)\u00E2\u008C\u00A9 (\u00CF\u0081cp)s \u00E2\u0088\u0082Ts \u00E2\u0088\u0082t \u00E2\u008C\u00AA = \u00E2\u0088\u0092\u00E3\u0080\u0088\u00E2\u0088\u0087 \u00C2\u00B7 JTs\u00E3\u0080\u0089 (2.13b)\u00E2\u008C\u00A9 (\u00CF\u0081cp)o \u00E2\u0088\u0082To \u00E2\u0088\u0082t \u00E2\u008C\u00AA = \u00E2\u0088\u0092\u00E3\u0080\u0088\u00E2\u0088\u0087 \u00C2\u00B7 JTo\u00E3\u0080\u0089 . (2.13c) 23 2.3. Temperature We will deal with the various fluxes individually. First consider the left-hand side of (2.13a), \u00E2\u008C\u00A9 (\u00CF\u0081cp)r \u00E2\u0088\u0082Tr \u00E2\u0088\u0082t \u00E2\u008C\u00AA = (\u00CF\u0081cp)r Vr V \u00E2\u0088\u0082 \u00E2\u0088\u0082t \u00EF\u00A3\u00AB\u00EF\u00A3\u00AD 1 Vr \u00E2\u0088\u00AB Vr Tr dV \u00EF\u00A3\u00B6\u00EF\u00A3\u00B8 = (1\u00E2\u0088\u0092 \u00CF\u0086)(\u00CF\u0081cp)r \u00E2\u0088\u0082 \u00E2\u0088\u0082t \u00E3\u0080\u0088Tr\u00E3\u0080\u0089\u00E2\u0088\u0097 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 sim- plification 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, \u00E2\u008C\u00A9 (\u00CF\u0081cp)s \u00E2\u0088\u0082Ts \u00E2\u0088\u0082t \u00E2\u008C\u00AA = (\u00CF\u0081cp)s \u00E2\u0088\u0082 \u00E2\u0088\u0082t \u00EF\u00A3\u00AB\u00EF\u00A3\u00AD 1 V \u00E2\u0088\u00AB Vs Ts dV \u00EF\u00A3\u00B6\u00EF\u00A3\u00B8 = \u00CF\u0086(\u00CF\u0081cp)s \u00E2\u0088\u0082 \u00E2\u0088\u0082t (s \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097) 24 2.3. Temperature where we have removed \u00CF\u0086 from the derivative under the assumption the porosity is constant. We can use the same methodology for (2.13c), \u00E2\u008C\u00A9 (\u00CF\u0081cp)o \u00E2\u0088\u0082To \u00E2\u0088\u0082t \u00E2\u008C\u00AA = \u00CF\u0086(\u00CF\u0081cp)o \u00E2\u0088\u0082 \u00E2\u0088\u0082t ((1\u00E2\u0088\u0092 s) \u00E3\u0080\u0088To\u00E3\u0080\u0089\u00E2\u0088\u0097) . 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) = (\u00CF\u0081cp)sqsTs Jconv(o) = (\u00CF\u0081cp)oqoTo where qi are the velocities of each fluid. We note from (2.4) that, \u00E2\u0088\u0087 \u00C2\u00B7 (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, \u00E3\u0080\u0088\u00E2\u0088\u0087 \u00C2\u00B7 ((\u00CF\u0081cp)sqsTs)\u00E3\u0080\u0089 = (\u00CF\u0081cp)s\u00E2\u0088\u0087 \u00C2\u00B7 \u00E3\u0080\u0088qsTs\u00E3\u0080\u0089+ (\u00CF\u0081cp)s 1 V \u00E2\u0088\u00AB Aint qsTs \u00C2\u00B7 n dA (2.14) 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 = \u00E3\u0080\u0088wi\u00E3\u0080\u0089\u00E2\u0088\u0097 + w\u00CC\u0083i which says that the true microscopic quantity is the intrinsic volume average of that quantity plus any spatial variations from the average denoted w\u00CC\u0083i. Notice if we take the intrinsic volume average of that definition, \u00E3\u0080\u0088wi\u00E3\u0080\u0089\u00E2\u0088\u0097 = \u00E3\u0080\u0088wi\u00E3\u0080\u0089\u00E2\u0088\u0097 + \u00E3\u0080\u0088w\u00CC\u0083i\u00E3\u0080\u0089\u00E2\u0088\u0097 \u00E3\u0080\u0088w\u00CC\u0083i\u00E3\u0080\u0089\u00E2\u0088\u0097 = 0. where we have used the fact that over Vi, \u00E3\u0080\u0088wi\u00E3\u0080\u0089\u00E2\u0088\u0097 is constant. Using these two facts we can write, qsTs = \u00E3\u0080\u0088qs\u00E3\u0080\u0089\u00E2\u0088\u0097 \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097 + \u00E3\u0080\u0088qs\u00E3\u0080\u0089\u00E2\u0088\u0097 T\u00CC\u0083s + \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097 q\u00CC\u0083s + q\u00CC\u0083sT\u00CC\u0083s. If we then take the intrinsic volume average of this, use the identities for w\u00CC\u0083i and the fact that a volume averaged quantity is constant over its intrinsic 26 2.3. Temperature volume we get, \u00E3\u0080\u0088qsTs\u00E3\u0080\u0089\u00E2\u0088\u0097 = \u00E3\u0080\u0088qs\u00E3\u0080\u0089\u00E2\u0088\u0097 \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097 + \u00E2\u008C\u00A9 q\u00CC\u0083sT\u00CC\u0083s \u00E2\u008C\u00AA\u00E2\u0088\u0097 \u00E3\u0080\u0088qsTs\u00E3\u0080\u0089 = \u00CF\u0086s ( \u00E3\u0080\u0088qs\u00E3\u0080\u0089\u00E2\u0088\u0097 \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097 + \u00E2\u008C\u00A9 q\u00CC\u0083sT\u00CC\u0083s \u00E2\u008C\u00AA\u00E2\u0088\u0097) . We may substitute this expression into (2.14) to get Jconv(s) = \u00E3\u0080\u0088\u00E2\u0088\u0087 \u00C2\u00B7 ((\u00CF\u0081cp)sqsTs)\u00E3\u0080\u0089 = (\u00CF\u0081cp)s\u00E2\u0088\u0087 \u00C2\u00B7 (\u00CF\u0086s \u00E3\u0080\u0088qs\u00E3\u0080\u0089\u00E2\u0088\u0097 \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097) + (\u00CF\u0081cp)s\u00E2\u0088\u0087 \u00C2\u00B7 ( \u00CF\u0086s \u00E2\u008C\u00A9 q\u00CC\u0083sT\u00CC\u0083s \u00E2\u008C\u00AA\u00E2\u0088\u0097) + (\u00CF\u0081cp)s 1 V \u00E2\u0088\u00AB Aint qsTs \u00C2\u00B7 n dA. We can define a similar convective flux for oil, Jconv(o) = \u00E3\u0080\u0088\u00E2\u0088\u0087 \u00C2\u00B7 ((\u00CF\u0081cp)oqoTo)\u00E3\u0080\u0089 = (\u00CF\u0081cp)o\u00E2\u0088\u0087 \u00C2\u00B7 (\u00CF\u0086(1\u00E2\u0088\u0092 s) \u00E3\u0080\u0088qo\u00E3\u0080\u0089\u00E2\u0088\u0097 \u00E3\u0080\u0088To\u00E3\u0080\u0089\u00E2\u0088\u0097) + (\u00CF\u0081cp)o\u00E2\u0088\u0087 \u00C2\u00B7 ( \u00CF\u0086(1\u00E2\u0088\u0092 s) \u00E2\u008C\u00A9 q\u00CC\u0083oT\u00CC\u0083o \u00E2\u008C\u00AA\u00E2\u0088\u0097) + (\u00CF\u0081cp)o 1 V \u00E2\u0088\u00AB Aint qoTo \u00C2\u00B7 n dA. 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 \u00C2\u00B7 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, \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097 = \u00E3\u0080\u0088To\u00E3\u0080\u0089\u00E2\u0088\u0097 = \u00E3\u0080\u0088T \u00E3\u0080\u0089 . In order to produce closure arguments, it is typical ([21],[30]) to rewrite the quantity T\u00CC\u0083i as, T\u00CC\u0083i = bi \u00C2\u00B7 \u00E2\u0088\u0087 \u00E3\u0080\u0088T \u00E3\u0080\u0089 with bi = bi(x) a vector in order to allow for anisotropy. This allows us to write the variations, T\u00CC\u0083i 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 \u00E2\u0088\u00AB Vi q\u00CC\u0083iT\u00CC\u0083i dV = 1 Vi \u00E2\u0088\u00AB Vi q\u00CC\u0083ibi \u00C2\u00B7 \u00E2\u0088\u0087 \u00E3\u0080\u0088T \u00E3\u0080\u0089 dV = \u00E3\u0080\u0088q\u00CC\u0083ibi\u00E3\u0080\u0089\u00E2\u0088\u0097 \u00C2\u00B7 \u00E2\u0088\u0087 \u00E3\u0080\u0088T \u00E3\u0080\u0089 28 2.3. Temperature where the product of two vectors u and v is to be understood as a dyadic product, uv = \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD u1v1 u1v2 u1v3 u2v1 u2v2 u2v3 u3v1 u3v2 u3v3 \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8 . 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, D = (\u00CF\u0081cp)s\u00CF\u0086s \u00E3\u0080\u0088q\u00CC\u0083sbs\u00E3\u0080\u0089\u00E2\u0088\u0097 + (\u00CF\u0081cp)o\u00CF\u0086(1\u00E2\u0088\u0092 s) \u00E3\u0080\u0088q\u00CC\u0083obo\u00E3\u0080\u0089\u00E2\u0088\u0097 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 = \u00E2\u0088\u0092nji, Jconv = \u00E2\u0088\u0087 \u00C2\u00B7 ([(\u00CF\u0081cp)s\u00CF\u0086s \u00E3\u0080\u0088qs\u00E3\u0080\u0089\u00E2\u0088\u0097 + (\u00CF\u0081cp)o\u00CF\u0086(1\u00E2\u0088\u0092 s) \u00E3\u0080\u0088qo\u00E3\u0080\u0089\u00E2\u0088\u0097] \u00E3\u0080\u0088T \u00E3\u0080\u0089) +\u00E2\u0088\u0087 \u00C2\u00B7 (D \u00C2\u00B7 \u00E2\u0088\u0087 \u00E3\u0080\u0088T \u00E3\u0080\u0089) . (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\u00E2\u0080\u0099s Law, Jcond(s) = \u00E2\u0088\u0092ks\u00E2\u0088\u0087Ts, Jcond(o) = \u00E2\u0088\u0092ko\u00E2\u0088\u0087To, Jcond(r) = \u00E2\u0088\u0092kr\u00E2\u0088\u0087Tr, 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 conduc- tivity 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, \u00E3\u0080\u0088\u00E2\u0088\u0087 \u00C2\u00B7 (ks\u00E2\u0088\u0087Ts)\u00E3\u0080\u0089 = \u00E2\u0088\u0087 \u00C2\u00B7 \u00E3\u0080\u0088ks\u00E2\u0088\u0087Ts\u00E3\u0080\u0089+ ks V \u00E2\u0088\u00AB Aint \u00E2\u0088\u0087Ts \u00C2\u00B7 n dA (2.16) 30 2.3. Temperature using identity (2.1b). Since we assume that ks is constant we can use identity (2.1a) to write, \u00E3\u0080\u0088ks\u00E2\u0088\u0087Ts\u00E3\u0080\u0089 = ks\u00E2\u0088\u0087\u00E3\u0080\u0088Ts\u00E3\u0080\u0089+ ks 1 V \u00E2\u0088\u00AB Aint Tsn dA = ks\u00E2\u0088\u0087(\u00CF\u0086s \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097) + ks 1 V \u00E2\u0088\u00AB Aint Tsn dA = ks\u00CF\u0086s\u00E2\u0088\u0087\u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097 + ks \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097\u00E2\u0088\u0087(\u00CF\u0086s) + ks 1 V \u00E2\u0088\u00AB Aint Tsn dA (2.17) 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 devi- ation via, Ts = \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097 + T\u00CC\u0083s and do a surface integral over the interfacial area, 1 V \u00E2\u0088\u00AB Aint Tsn dA = 1 V \u00E2\u0088\u00AB Aint \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097 n dA+ 1 V \u00E2\u0088\u00AB Aint T\u00CC\u0083sn dA. (2.18) Now for the second term in this expression we can reverse identity (2.1a) to get that, 1 V \u00E2\u0088\u00AB Aint \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097 n dA = \u00E3\u0080\u0088\u00E2\u0088\u0087 \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097\u00E3\u0080\u0089 \u00E2\u0088\u0092 \u00E2\u0088\u0087 \u00E3\u0080\u0088\u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097\u00E3\u0080\u0089 . We can write, \u00E2\u0088\u0087\u00E3\u0080\u0088\u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097\u00E3\u0080\u0089 = \u00E2\u0088\u0087 (\u00CF\u0086s \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097) = \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097\u00E2\u0088\u0087(\u00CF\u0086s) + \u00CF\u0086s\u00E2\u0088\u0087\u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097 31 2.3. Temperature by utilizing the fact that the intrinsic volume average is constant over its subdomain volume. Now we may write, \u00E3\u0080\u0088\u00E2\u0088\u0087 \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097\u00E3\u0080\u0089 = \u00CF\u0086s\u00E2\u0088\u0087\u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097 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 \u00E2\u0088\u00AB Aint Tsn dA = 1 V \u00E2\u0088\u00AB Aint T\u00CC\u0083sn dA\u00E2\u0088\u0092 \u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097\u00E2\u0088\u0087(\u00CF\u0086s). 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, \u00E3\u0080\u0088ks\u00E2\u0088\u0087Ts\u00E3\u0080\u0089 = ks\u00CF\u0086s\u00E2\u0088\u0087\u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097 + ks V \u00E2\u0088\u00AB Aint T\u00CC\u0083sn dA so that the conductive flux for water is, Jcond(s) = \u00E2\u0088\u0092\u00E2\u0088\u0087 \u00C2\u00B7 (ks\u00CF\u0086s\u00E2\u0088\u0087\u00E3\u0080\u0088Ts\u00E3\u0080\u0089\u00E2\u0088\u0097)\u00E2\u0088\u0092\u00E2\u0088\u0087 \u00C2\u00B7 \u00EF\u00A3\u00AB\u00EF\u00A3\u00ADks V \u00E2\u0088\u00AB Aint T\u00CC\u0083sn dA \u00EF\u00A3\u00B6\u00EF\u00A3\u00B8\u00E2\u0088\u0092 ks V \u00E2\u0088\u00AB Aint \u00E2\u0088\u0087Ts \u00C2\u00B7 n dA. 32 2.3. Temperature In a similar fashion we can derive conductive fluxes for oil and the rock, Jcond(o) = \u00E2\u0088\u0092\u00E2\u0088\u0087 \u00C2\u00B7 (ko\u00CF\u0086(1\u00E2\u0088\u0092 s)\u00E2\u0088\u0087\u00E3\u0080\u0088To\u00E3\u0080\u0089\u00E2\u0088\u0097)\u00E2\u0088\u0092\u00E2\u0088\u0087 \u00C2\u00B7 \u00EF\u00A3\u00AB\u00EF\u00A3\u00ADko V \u00E2\u0088\u00AB Aint T\u00CC\u0083on dA \u00EF\u00A3\u00B6\u00EF\u00A3\u00B8 \u00E2\u0088\u0092 ko V \u00E2\u0088\u00AB Aint \u00E2\u0088\u0087To \u00C2\u00B7 n dA, Jcond(r) = \u00E2\u0088\u0092\u00E2\u0088\u0087 \u00C2\u00B7 (kr(1\u00E2\u0088\u0092 \u00CF\u0086)\u00E2\u0088\u0087\u00E3\u0080\u0088Tr\u00E3\u0080\u0089\u00E2\u0088\u0097)\u00E2\u0088\u0092\u00E2\u0088\u0087 \u00C2\u00B7 \u00EF\u00A3\u00AB\u00EF\u00A3\u00ADkr V \u00E2\u0088\u00AB Aint T\u00CC\u0083rn dA \u00EF\u00A3\u00B6\u00EF\u00A3\u00B8 \u00E2\u0088\u0092 kr V \u00E2\u0088\u00AB Aint \u00E2\u0088\u0087Tr \u00C2\u00B7 n dA. Now as for the convective fluxes we will assume that we have thermal equi- librium 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 \u00C2\u00B7 (ki\u00E2\u0088\u0087Ti) = nij \u00C2\u00B7 (kj\u00E2\u0088\u0087Tj). Notice that the first boundary condition along with thermal equilibrium im- plies that T\u00CC\u0083i = T\u00CC\u0083j. Now like with convection, we will rewrite the temperature variations by, T\u00CC\u0083i = bi \u00C2\u00B7 \u00E2\u0088\u0087 \u00E3\u0080\u0088T \u00E3\u0080\u0089 33 2.3. Temperature so that, \u00E2\u0088\u00AB Aint T\u00CC\u0083in dA = \u00E2\u0088\u00AB Aint bin dA\u00E2\u0088\u0087\u00E3\u0080\u0088T \u00E3\u0080\u0089 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 = \u00E2\u0088\u0092\u00E2\u0088\u0087 \u00C2\u00B7 ((ks\u00CF\u0086s+ ko\u00CF\u0086(1\u00E2\u0088\u0092 s) + kr(1\u00E2\u0088\u0092 \u00CF\u0086))\u00E2\u0088\u0087\u00E3\u0080\u0088T \u00E3\u0080\u0089+ K \u00C2\u00B7 \u00E2\u0088\u0087 \u00E3\u0080\u0088T \u00E3\u0080\u0089) (2.19) where, K = ks \u00E2\u0088\u0092 ko V \u00E2\u0088\u00AB Aso bsnso dA+ ks \u00E2\u0088\u0092 kr V \u00E2\u0088\u00AB Asr bsnsr dA+ ko \u00E2\u0088\u0092 kr V \u00E2\u0088\u00AB Aor bonor dA is the conduction tensor. Notice we could incorporate all of the terms mul- tiplying 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 re- formulation using \u00E2\u0088\u0087\u00E3\u0080\u0088T \u00E3\u0080\u0089 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 con- vective (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, \u00CF\u0086(\u00CF\u0081cp)s \u00E2\u0088\u0082 \u00E2\u0088\u0082t (s \u00E3\u0080\u0088T \u00E3\u0080\u0089) + \u00CF\u0086(\u00CF\u0081cp)o \u00E2\u0088\u0082 \u00E2\u0088\u0082t ((1\u00E2\u0088\u0092 s) \u00E3\u0080\u0088T \u00E3\u0080\u0089) + (1\u00E2\u0088\u0092 \u00CF\u0086)(\u00CF\u0081cp)r \u00E2\u0088\u0082 \u00E3\u0080\u0088T \u00E3\u0080\u0089 \u00E2\u0088\u0082t +\u00E2\u0088\u0087 \u00C2\u00B7 ([(\u00CF\u0081cp)s\u00CF\u0086s \u00E3\u0080\u0088qs\u00E3\u0080\u0089\u00E2\u0088\u0097 + (\u00CF\u0081cp)o\u00CF\u0086(1\u00E2\u0088\u0092 s) \u00E3\u0080\u0088qo\u00E3\u0080\u0089\u00E2\u0088\u0097] \u00E3\u0080\u0088T \u00E3\u0080\u0089) = \u00E2\u0088\u0087 \u00C2\u00B7 ((ks\u00CF\u0086s+ ko\u00CF\u0086(1\u00E2\u0088\u0092 s) + kr(1\u00E2\u0088\u0092 \u00CF\u0086))\u00E2\u0088\u0087\u00E3\u0080\u0088T \u00E3\u0080\u0089) +\u00E2\u0088\u0087 \u00C2\u00B7 ((K\u00E2\u0088\u0092D) \u00C2\u00B7 \u00E2\u0088\u0087 \u00E3\u0080\u0088T \u00E3\u0080\u0089) . Now if we expand the time derivative out and use (2.3) for the following, \u00CF\u0086(\u00CF\u0081cp)s \u00E2\u0088\u0082s \u00E2\u0088\u0082t = \u00E2\u0088\u0092(\u00CF\u0081cp)s\u00E2\u0088\u0087 \u00C2\u00B7 \u00E3\u0080\u0088qs\u00E3\u0080\u0089 \u00CF\u0086(\u00CF\u0081cp)o \u00E2\u0088\u0082s \u00E2\u0088\u0082t = (\u00CF\u0081cp)o\u00E2\u0088\u0087 \u00C2\u00B7 \u00E3\u0080\u0088qo\u00E3\u0080\u0089 and take advantage of the fact that, \u00CF\u0086s \u00E3\u0080\u0088qs\u00E3\u0080\u0089\u00E2\u0088\u0097 = \u00E3\u0080\u0088qs\u00E3\u0080\u0089 \u00CF\u0086(1\u00E2\u0088\u0092 s) \u00E3\u0080\u0088qo\u00E3\u0080\u0089\u00E2\u0088\u0097 = \u00E3\u0080\u0088qo\u00E3\u0080\u0089 then we can drastically simplify the problem, (\u00CF\u0086(\u00CF\u0081cp)ss+ \u00CF\u0086(\u00CF\u0081cp)o(1\u00E2\u0088\u0092 s) + (1\u00E2\u0088\u0092 \u00CF\u0086)(\u00CF\u0081cp)r) \u00E2\u0088\u0082T \u00E2\u0088\u0082t + (\u00CF\u0086(\u00CF\u0081cp)sqs + \u00CF\u0086(\u00CF\u0081cp)oqo) \u00C2\u00B7 \u00E2\u0088\u0087T = \u00E2\u0088\u0087 \u00C2\u00B7 ((ks\u00CF\u0086s+ ko\u00CF\u0086(1\u00E2\u0088\u0092 s) + kr(1\u00E2\u0088\u0092 \u00CF\u0086))\u00E2\u0088\u0087T ) +\u00E2\u0088\u0087 \u00C2\u00B7 ((K\u00E2\u0088\u0092D) \u00C2\u00B7 \u00E2\u0088\u0087T ) . (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 en- compassing 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 expres- sion 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, \u00C2\u00B5i = \u00C2\u00B5i(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, \u00C2\u00B5s(T ) = \u00C2\u00B5s0 \u00C2\u00B5o(T ) = (\u00C2\u00B5o0 \u00E2\u0088\u0092 \u00C2\u00B5s0) exp ( \u00E2\u0088\u0092(T \u00E2\u0088\u0092 Ta) Ta ) + \u00C2\u00B5s0 36 2.3. Temperature where \u00C2\u00B5i0 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 vis- cosities [2]. 2.3.3 Non-Dimensionalization First we consider the same rescaling as before, x\u00CC\u0084 = x L , t\u00CC\u0084 = t \u00CF\u0084 , q\u00CC\u0084i = qi Q , \u00CE\u00B1i = \u00CE\u00BAi \u00CE\u00BA where all parameters are as defined in Table 2.1. We will rescale temperature by an ambient value, Ta, T\u00CC\u0084 = T Ta consistent with the introduction of a temperature dependent viscosity. We will also scale the viscosities by their respective reference values, \u00C2\u00B5\u00CC\u0084i = \u00C2\u00B5i \u00C2\u00B5i0 . With the use of this rescaling we may utilize the same scaling for pressure as in the temperature independent model, P\u00CC\u0084i = Pi \u00CE\u00BA \u00C2\u00B5s0QL 37 2.3. Temperature where we now consider the reference values since viscosity is no longer con- stant. Our non-dimensional equations can then be written as (removing overbars), \u00E2\u0088\u0087 \u00C2\u00B7 (qs + qo) = 0, (2.21a) \u00CF\u0086Fl \u00E2\u0088\u0082s \u00E2\u0088\u0082t +\u00E2\u0088\u0087 \u00C2\u00B7 qs = 0, (2.21b) \u00C2\u00B5s(T ) = 1, (2.21c) \u00C2\u00B5o(T ) = (1\u00E2\u0088\u0092 \u00CE\u00B4) exp (1\u00E2\u0088\u0092 T ) + \u00CE\u00B4, (2.21d) qs = \u00E2\u0088\u0092\u00CE\u00B1s(s)\u00E2\u0088\u0087Ps, (2.21e) qo = \u00E2\u0088\u0092\u00CE\u00B4 \u00CE\u00B1o(s) \u00C2\u00B5o(T ) \u00E2\u0088\u0087Po, (2.21f) Po \u00E2\u0088\u0092 Ps = Ca\u00E2\u0088\u00921J(s), (2.21g) \u00E2\u0088\u0082T \u00E2\u0088\u0082t + (CsC\u00CC\u0082s(s)qs + CoC\u00CC\u0082o(s)qo) \u00C2\u00B7 \u00E2\u0088\u0087T = D\u00E2\u0088\u0087 \u00C2\u00B7 (D\u00CC\u0082(s)\u00E2\u0088\u0087T ) (2.21h) where the nondimensional numbers \u00CE\u00B4, 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 = \u00CF\u0086(\u00CF\u0081cp)s Fl \u00E2\u0088\u00AB 1 0 (\u00CF\u0086(\u00CF\u0081cp)ss+ \u00CF\u0086(\u00CF\u0081cp)o(1\u00E2\u0088\u0092 s) + (1\u00E2\u0088\u0092 \u00CF\u0086)(\u00CF\u0081cp)r) ds (2.22) Co = \u00CF\u0086(\u00CF\u0081cp)o Fl \u00E2\u0088\u00AB 1 0 (\u00CF\u0086(\u00CF\u0081cp)ss+ \u00CF\u0086(\u00CF\u0081cp)o(1\u00E2\u0088\u0092 s) + (1\u00E2\u0088\u0092 \u00CF\u0086)(\u00CF\u0081cp)r) ds (2.23) 38 2.3. Temperature and the diffusion number D as the average ratio of diffusive energy transport to total heat energy, D = \u00E2\u0088\u00AB 1 0 (ks\u00CF\u0086s+ ko\u00CF\u0086(1\u00E2\u0088\u0092 s) + kr(1\u00E2\u0088\u0092 \u00CF\u0086)) ds QLFl \u00E2\u0088\u00AB 1 0 (\u00CF\u0086(\u00CF\u0081cp)ss+ \u00CF\u0086(\u00CF\u0081cp)o(1\u00E2\u0088\u0092 s) + (1\u00E2\u0088\u0092 \u00CF\u0086)(\u00CF\u0081cp)r) ds . (2.24) Notice that the diffusion number acts like an averaged inverse Peclet number. The functions C\u00CC\u0082i and D\u00CC\u0082 are linear order one functions that contain the s dependence. They are order one because 0 \u00E2\u0089\u00A4 s \u00E2\u0089\u00A4 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. Symbol Value Description Reference s o r s o r cp 4.18\u00C3\u0097103 1.76\u00C3\u0097103 794 Heat Capacity [15] [28] [4](Jkg\u00E2\u0088\u00921K\u00E2\u0088\u00921) \u00CF\u0081 1000 877 2.65\u00C3\u0097103 Density [14] [26] [3] (kgm\u00E2\u0088\u00923) k 0.610 0.151 1.04 Thermal [33] [26] [44]Conductivity (Wm\u00E2\u0088\u00921K\u00E2\u0088\u00921) Ta 300K Reference Chosen Temperature \u00C2\u00B5s0 1.00\u00C3\u009710\u00E2\u0088\u00923Pa\u00C2\u00B7s Reference [23]Water Viscosity \u00C2\u00B5o0 0.523Pa\u00C2\u00B7s Reference ChosenOil Viscosity Cs 0.405 Water Eqn (2.22)Convection Number Co 0.150 Oil Eqn (2.23)Convection Number D 7.24\u00C3\u009710\u00E2\u0088\u009210 Diffusion 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 \u000F\u001C 1. In this case, the water and oil pressures are the same and we denote them as P to get a 1D model, \u00E2\u0088\u0082 \u00E2\u0088\u0082x (qs + qo) = 0, (3.1a) \u00E2\u0088\u0082s \u00E2\u0088\u0082t + \u00E2\u0088\u0082qs \u00E2\u0088\u0082x = 0, (3.1b) qs = \u00E2\u0088\u0092\u00CE\u00B1s(s)dP dx , (3.1c) qo = \u00E2\u0088\u0092\u00CE\u00B4\u00CE\u00B1o(s)dP dx . (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, \u00E2\u0088\u0082s \u00E2\u0088\u0082t + \u00E2\u0088\u0082 \u00E2\u0088\u0082x f(s) = 0, \u00E2\u0088\u0092\u00E2\u0088\u009E < x <\u00E2\u0088\u009E, t > 0 (3.3) where, for the relative permeability functions assumed here, f(s) = s3 s3 + \u00CE\u00B4(1\u00E2\u0088\u0092 s)3 = \u00CE\u00B1s \u00CE\u00BB (3.4) with, \u00CE\u00BB = \u00CE\u00B1s + \u00CE\u00B4\u00CE\u00B1o = s 3 + \u00CE\u00B4(1\u00E2\u0088\u0092 s)3 being denoted as the total permeability. Note that this is the famous Buckley- Leverett 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) = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B31 x < 00 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 x\u00E2\u0086\u0092\u00E2\u0088\u009E s(x, t) = 0 lim x\u00E2\u0086\u0092\u00E2\u0088\u0092\u00E2\u0088\u009E s(x, t) = 1. 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 \u00E2\u0080\u00B2\u00E2\u0080\u00B2(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 dt = 0 along characteristics that satisfy, dx dt = f \u00E2\u0080\u00B2(s). 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 = \u00CE\u00B2 44 3.1. 1D Model in the Absence of Capillary Pressure where x(0) = \u00CE\u00B2 and the solution is then, s(\u00CE\u00B2) = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B31 \u00CE\u00B2 < 00 \u00CE\u00B2 > 0 . All the characteristics that satisfy \u00CE\u00B2 < 0 have s = 1 while \u00CE\u00B2 > 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 \u00CE\u00B2 = 0 and so, knowing that the characteristic lines that form the rarefaction must emerge from the origin, we seek a solution, s = s (x t ) = s(\u00CF\u0087) there. Substituting this into (3.3) we get that, ds d\u00CF\u0087 (f \u00E2\u0080\u00B2(s)\u00E2\u0088\u0092 \u00CF\u0087) = 0. Since we expect that s = s(\u00CF\u0087) we assume that the derivative is non-zero and thus we get for our rarefaction fan, f \u00E2\u0080\u00B2(s) = x t , 0 \u00E2\u0089\u00A4 s \u00E2\u0089\u00A4 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 t = 0 x = 0 x = \u00CE\u00BE s = 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 = \u00CE\u00BE(t) forms (Figure 3.1). We know that the conditions across the shock must satisfy a Rankine-Hugoniot condition ([31]), v = d\u00CE\u00BE dt = f(s+)\u00E2\u0088\u0092 f(s\u00E2\u0088\u0092) s+ \u00E2\u0088\u0092 s\u00E2\u0088\u0092 where v is the shock velocity, s\u00E2\u0088\u0092 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 \u00E2\u0080\u00B2(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\u00E2\u0088\u0092 = s\u00E2\u0088\u0097 < 1 that satisfies, f \u00E2\u0080\u00B2(s\u00E2\u0088\u0097) = f(s\u00E2\u0088\u0097) s\u00E2\u0088\u0097 . Now using the expression for f(s) we get that a non-trivial s\u00E2\u0088\u0097 must satisfy (\u00CE\u00B4 \u00E2\u0088\u0092 1)s\u00E2\u0088\u00973 \u00E2\u0088\u0092 3\u00CE\u00B4s\u00E2\u0088\u0097 + 2\u00CE\u00B4 = 0. (3.5) In the case where \u00CE\u00B4 \u001C 1 then (\u00CE\u00B4 \u00E2\u0088\u0092 1) \u00E2\u0089\u0088 \u00E2\u0088\u0092 1 and so, s\u00E2\u0088\u00973 + 3\u00CE\u00B4s\u00E2\u0088\u0097 \u00E2\u0088\u0092 2\u00CE\u00B4 \u00E2\u0089\u0088 0. To leading order in \u00CE\u00B4, this only produces trivial roots but if we scale s\u00E2\u0088\u0097 = \u00CE\u00B41/3k\u00E2\u0088\u0097 and expand, k\u00E2\u0088\u0097 = k0 + \u00CE\u00B41/3k1 then we get that the asymptotic shock saturation value s\u00E2\u0088\u0097a is, s\u00E2\u0088\u0097a \u00E2\u0088\u00BC 21/3\u00CE\u00B41/3 \u00E2\u0088\u0092 2\u00E2\u0088\u00921/3\u00CE\u00B42/3. (3.6) For the specific choice of \u00CE\u00B4 from Table 2.1, we get that s\u00E2\u0088\u0097a \u00E2\u0089\u0088 0.1442 and by numerically solving (3.5) we get the actual shock saturation is s\u00E2\u0088\u0097 \u00E2\u0089\u0088 0.1443, 47 3.1. 1D Model in the Absence of Capillary Pressure a good agreement. Now for the shock velocity, v = f(s\u00E2\u0088\u0097) s\u00E2\u0088\u0097 = s\u00E2\u0088\u00972 s\u00E2\u0088\u00973 + \u00CE\u00B4 (1\u00E2\u0088\u0092 s\u00E2\u0088\u0097)3 (3.7) which we can approximate using the asymptotic shock saturation value. For simplification, we will take only the first term in s\u00E2\u0088\u0097a which will provide accu- racy up to o(\u00CE\u00B41/3). By substituting this into the expression for shock velocity (3.7) and expanding for \u00CE\u00B4 \u001C 1 we get that the asymptotic shock velocity va is, va \u00E2\u0088\u00BC 2 2/3 3\u00CE\u00B41/3 + 2 3 . For our choice of parameters, we get that va \u00E2\u0089\u0088 4.93 while numerically we calculate from (3.7) that v \u00E2\u0089\u0088 4.95, another good agreement. By finding s\u00E2\u0088\u0097 and v, the problem is complete and has solution, s(x, t) = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3 1 x \u00E2\u0089\u00A4 0 [f \u00E2\u0080\u00B2]\u00E2\u0088\u00921 ( x t ) 0 < x < f \u00E2\u0080\u00B2(s\u00E2\u0088\u0097)t 0 x > f \u00E2\u0080\u00B2(s\u00E2\u0088\u0097)t . (3.8) For the \u00CE\u00B4 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 -2 -1 0 1 2 3 4 5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Saturation Profile at t=0.125 x s Figure 3.2: Saturation profile from solution (3.8) with s\u00E2\u0088\u0097 = 0.1443 and t = 1/8. Now, notice that if we consider a moving coordinate frame, \u00CE\u00B7 = x\u00E2\u0088\u0092 vt then (3.3) becomes, \u00E2\u0088\u0082s \u00E2\u0088\u0082t + \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00CE\u00B7 (f(s)\u00E2\u0088\u0092 vs) = 0 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\u00E2\u0088\u00AB d d\u00CE\u00B7 (f(s0)\u00E2\u0088\u0092 vs0) d\u00CE\u00B7 = f(s0)\u00E2\u0088\u0092 vs0 = C, 49 3.1. 1D Model in the Absence of Capillary Pressure a constant. When C = 0 then f(s\u00E2\u0088\u0097) s\u00E2\u0088\u0097 = v satisfies this condition. These are precisely the conditions that describe the shock and thus the shock front, s0 = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3s \u00E2\u0088\u0097, \u00CE\u00B7 < 0 s\u00E2\u0088\u0092, \u00CE\u00B7 > 0 (3.9a) 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\u00E2\u0088\u0097 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 d\u00CE\u00B7 = \u00E2\u0088\u0092 1 s30 + \u00CE\u00B4 (1\u00E2\u0088\u0092 s0)3 = \u00E2\u0088\u00921 \u00CE\u00BB , \u00CE\u00B7 6= 0. (3.9b) We will now analyze the stability of this travelling wave solution under per- turbations 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 \u00E2\u0088\u00BC s0 + \u00CE\u00B5s1(\u00CE\u00B7) exp (iny + \u00CF\u0083t) , \u00CE\u00B7 6= 0 (3.10a) P \u00E2\u0088\u00BC P0 + \u00CE\u00B5P1(\u00CE\u00B7) exp (iny + \u00CF\u0083t) , \u00CE\u00B7 6= 0 (3.10b) s1, ds1 d\u00CE\u00B7 , P1, dP1 d\u00CE\u00B7 \u00E2\u0086\u0092 0, |\u00CE\u00B7| \u00E2\u0086\u0092 \u00E2\u0088\u009E (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 \u00CF\u0083 then determines the stability of the shock to disturbances of wave number n. The parameter \u00CE\u00B5 is an artificial parameter introduced just to control the balance of terms, it is not to be confused with \u000F, the in- verse 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 di- rection 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 3It 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\u00E2\u0080\u0099s Law. 51 3.1. 1D Model in the Absence of Capillary Pressure excluding capillary pressure)4, \u00E2\u0088\u0082s \u00E2\u0088\u0082t +\u00E2\u0088\u0087 \u00C2\u00B7 (qs \u00E2\u0088\u0092 vsx\u00CC\u0082) = 0 qs = \u00E2\u0088\u0092\u00CE\u00B1s\u00E2\u0088\u0087P qo = \u00E2\u0088\u0092\u00CE\u00B4\u00CE\u00B1o\u00E2\u0088\u0087P \u00E2\u0088\u0087 \u00C2\u00B7 (qs + qo) = 0 where x\u00CC\u0082 is a unit vector in the \u00CE\u00B7 direction, then, keeping only linear terms (i.e. terms of O(\u00CE\u00B5)), we get for our eigenvalue problem, d d\u00CE\u00B7 [ \u00CE\u00B1\u00CC\u0084s dP1 d\u00CE\u00B7 + \u00CE\u00B1\u00CC\u0084s \u00E2\u0080\u00B2dP0 d\u00CE\u00B7 s1 + vs1 ] \u00E2\u0088\u0092 n2\u00CE\u00B1\u00CC\u0084sP1 = \u00CF\u0083s1 (3.11a) d d\u00CE\u00B7 [ \u00CE\u00BB\u00CC\u0084 dP1 d\u00CE\u00B7 + \u00CE\u00BB\u00CC\u0084\u00E2\u0080\u00B2 dP0 d\u00CE\u00B7 s1 ] \u00E2\u0088\u0092 n2\u00CE\u00BB\u00CC\u0084P1 = 0. (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 \u00CE\u00B7 = 0 is simply constant (3.9) so, \u00CE\u00B1\u00CC\u0084s, \u00CE\u00BB\u00CC\u0084, and its derivatives are all constants and can be removed from the derivatives. We can rearrange (3.11) solely for s1,( 1\u00E2\u0088\u0092 f \u00E2\u0080\u00B2(s0) v ) ds1 d\u00CE\u00B7 = \u00CF\u0083 v s1 4We have redefined the gradient operator as, \u00E2\u0088\u0087 = \u00E2\u008C\u00A9 \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00CE\u00B7 , \u00E2\u0088\u0082 \u00E2\u0088\u0082y \u00E2\u008C\u00AA . 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 \u00CE\u00B7 > 0 then s0 = 0 and we get that, s1 = D+ exp (\u00CF\u0083 v \u00CE\u00B7 ) but for any unstable modes5 (\u00CF\u0083 > 0) this cannot vanish as \u00CE\u00B7 \u00E2\u0086\u0092\u00E2\u0088\u009E and so we conclude that D+ = 0. If \u00CE\u00B7 < 0 then s0 = s \u00E2\u0088\u0097 < 1 and f \u00E2\u0080\u00B2(s\u00E2\u0088\u0097) = v so therefore, 0 = \u00CF\u0083 v s1. Thus for general \u00CF\u0083, s1 = 0. Therefore, s1 = 0, \u00CE\u00B7 6= 0. (3.12a) Now, knowing that s1 = 0 away from \u00CE\u00B7 = 0 means that we get a simplified expression for P1, d2P1 d\u00CE\u00B72 = n2P1 5Fingering 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 = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3C+ exp (\u00E2\u0088\u0092n\u00CE\u00B7) , \u00CE\u00B7 > 0C\u00E2\u0088\u0092 exp (n\u00CE\u00B7) , \u00CE\u00B7 < 0 (3.12b) 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 \u00CE\u00B7 = 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\u00CE\u00B4(\u00CE\u00B7) (3.13) in order to satisfy the appropriate conditions and conclusion about s1 away from \u00CE\u00B7 = 0. Upon doing so, they generate a series of jump conditions across the shock front that are solvability conditions for the eigenvalue, [ \u00CE\u00BB\u00CC\u0084 dP1 d\u00CE\u00B7 ] = 0, (3.14a)[ \u00CE\u00B1\u00CC\u0084s dP1 d\u00CE\u00B7 ] = s\u00E2\u0088\u0097A\u00CF\u0083, (3.14b) [P1] = \u00E2\u0088\u0092A [ dP0 d\u00CE\u00B7 ] , (3.14c) (3.14d) which we have generalized to suit our problem. Note that we let [\u00C2\u00B7] 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 \u00CE\u00B7 = 0 but also by perturbing the interface to get, s \u00E2\u0088\u00BC s0 + \u00CE\u00B5s1 exp (iny + \u00CF\u0083t) P \u00E2\u0088\u00BC P0 + \u00CE\u00B5P1 exp (iny + \u00CF\u0083t) \u00CE\u00B7\u00CC\u0082 \u00E2\u0088\u00BC 0 + \u00CE\u00B5A\u00CC\u0083 exp (iny + \u00CF\u0083t) where once again \u00CE\u00B5 is a small parameter introduced to control term balance, then we will still conclude that s1 = 0 and that P1 is exponential away from \u00CE\u00B7 = \u00CE\u00B7\u00CC\u0082 \u00E2\u0089\u0088 0 and thus the eigenfunctions are still given by (3.12). However at the boundary, \u00CE\u00B7\u00CC\u0082, 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 \u00E2\u0088\u0092 vsx\u00CC\u0082 qo \u00E2\u0088\u0092 v(1\u00E2\u0088\u0092 s)x\u00CC\u0082 for water and oil respectively. We find the normal fluxes by taking an in- ner product with the normal direction of the interface, n\u00CC\u0082. However at the interface \u00CE\u00B7 = \u00CE\u00B7\u00CC\u0082 an we have that the normal vector is, n\u00CC\u0082 = \u00E3\u0080\u00881, \u00CE\u00B5Ainy exp (iny + \u00CF\u0083t)\u00E3\u0080\u0089\u00E2\u0088\u009A 1 + (\u00CE\u00B5Ainy exp (iny + \u00CF\u0083t))2 55 3.1. 1D Model in the Absence of Capillary Pressure which we can expand as, n\u00CC\u0082 \u00E2\u0089\u0088 x\u00CC\u0082 + \u00CE\u00B5Ainy exp (iny + \u00CF\u0083t) y\u00CC\u0082 +O(\u00CE\u00B52) where y\u00CC\u0082 is a unit vector in the y-direction. However, to first order, the flux terms only have components in the x-direction, qi = \u00E3\u0080\u0088qi0, 0\u00E3\u0080\u0089+ \u00CE\u00B5\u00E3\u0080\u0088qi1x, qi1y\u00E3\u0080\u0089 then the only contribution for normal fluxes up to O(\u00CE\u00B5) is to have n\u00CC\u0082 = x\u00CC\u0082. Therefore by taking the inner product in the x\u00CC\u0082 direction, the total jump in water across the interface \u00CE\u00B7\u00CC\u0082 = 0 is, [x\u00CC\u0082 \u00C2\u00B7 qs \u00E2\u0088\u0092 vs] = ( total flux across \u00CE\u00B7 = \u00CE\u00B7\u00CC\u0082+ )\u00E2\u0088\u0092 (total flux across \u00CE\u00B7 = \u00CE\u00B7\u00CC\u0082\u00E2\u0088\u0092) \u00E2\u0088\u00BC (total flux across \u00CE\u00B7 = 0+)\u00E2\u0088\u0092 (total flux across \u00CE\u00B7 = 0\u00E2\u0088\u0092) = s(0+) d\u00CE\u00B7\u00CC\u0082 dt \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 \u00CE\u00B7=0+ \u00E2\u0088\u0092 s(0\u00E2\u0088\u0092) d\u00CE\u00B7\u00CC\u0082 dt \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 \u00CE\u00B7=0\u00E2\u0088\u0092 with, d\u00CE\u00B7\u00CC\u0082 dt = \u00CE\u00B5\u00CF\u0083A exp (iny + \u00CF\u0083t) being the speed of the interface. However we also know that, [x\u00CC\u0082 \u00C2\u00B7 qs \u00E2\u0088\u0092 vs] = [ \u00E2\u0088\u0092\u00CE\u00B1\u00CC\u0084sdP0 d\u00CE\u00B7 ] \u00E2\u0088\u0092 v [s0] + \u00CE\u00B5 [ \u00E2\u0088\u0092\u00CE\u00B1\u00CC\u0084sdP1 d\u00CE\u00B7 ] exp (iny + \u00CF\u0083t) +O(\u00CE\u00B52) where we have explicitly used the fact that, s1 = ds1 d\u00CE\u00B7 = d2P0 d\u00CE\u00B72 = 0. 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\u00CC\u0082 \u00C2\u00B7 qs \u00E2\u0088\u0092 vs] = [ \u00E2\u0088\u0092\u00CE\u00B1\u00CC\u0084sdP0 d\u00CE\u00B7 ] + \u00CE\u00B5 [ \u00E2\u0088\u0092\u00CE\u00B1\u00CC\u0084sdP1 d\u00CE\u00B7 ] exp (iny + \u00CF\u0083t) = \u00E2\u0088\u0092vs\u00E2\u0088\u0097 \u00E2\u0088\u0092 \u00CE\u00B5s\u00E2\u0088\u0097\u00CF\u0083A\u00CC\u0083 exp (iny + \u00CF\u0083t) which states that [ \u00CE\u00B1\u00CC\u0084s dP0 d\u00CE\u00B7 ] = vs\u00E2\u0088\u0097[ \u00CE\u00B1\u00CC\u0084s dP1 d\u00CE\u00B7 ] = s\u00E2\u0088\u0097\u00CF\u0083A\u00CC\u0083. The first of these expressions is satisfied by applying the known expression for the steady-state pressure and the second equation provides the first solvabil- ity condition for \u00CF\u0083. 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 \u00CE\u00B7 = \u00CE\u00B7\u00CC\u0082 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\u00CC\u0082 \u00C2\u00B7 qo \u00E2\u0088\u0092 v(1\u00E2\u0088\u0092 s)] = [ \u00E2\u0088\u0092\u00CE\u00B4\u00CE\u00B1\u00CC\u0084odP0 d\u00CE\u00B7 ] + \u00CE\u00B5 [ \u00E2\u0088\u0092\u00CE\u00B4\u00CE\u00B1\u00CC\u0084odP1 d\u00CE\u00B7 ] exp (iny + \u00CF\u0083t) = vs\u00E2\u0088\u0097 + \u00CE\u00B5s\u00E2\u0088\u0097\u00CF\u0083A\u00CC\u0083 exp (iny + \u00CF\u0083t) . 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, [x\u00CC\u0082 \u00C2\u00B7 (qs + qo)\u00E2\u0088\u0092 v] = [ \u00CE\u00BB\u00CC\u0084 dP0 d\u00CE\u00B7 ] + \u00CE\u00B5 [ \u00CE\u00BB\u00CC\u0084 dP1 d\u00CE\u00B7 ] exp (iny + \u00CF\u0083t) = 0. Here we see that the jump in total fluid flux across the boundary is con- served as we expected. From this condition we generate another solvability condition, [ \u00CE\u00BB\u00CC\u0084 dP1 d\u00CE\u00B7 ] = 0. Finally we look at the continuity of pressure. On either side of the interface we can expand P (\u00CE\u00B7\u00CC\u0082k) with \u00CE\u00B7\u00CC\u0082 \u00E2\u0089\u0088 0 where k = + or \u00E2\u0088\u0092 depending on the left or right near \u00CE\u00B7\u00CC\u0082 to get, P (\u00CE\u00B7\u00CC\u0082k) \u00E2\u0088\u00BC P k0 (0) + \u00CE\u00B5 dPo d\u00CE\u00B7 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 \u00CE\u00B7=0k A\u00CC\u0083 exp (iny + \u00CF\u0083t) + \u00CE\u00B5P k1 (0) exp (iny + \u00CF\u0083t) so that we can write the jump in pressure as, [P ] = [P0] + \u00CE\u00B5 ( [P1] + A\u00CC\u0083 [ dP0 d\u00CE\u00B7 ]) exp (iny + \u00CF\u0083t) = 0. This generates the conditions that, [P0] = 0 [P1] = \u00E2\u0088\u0092A\u00CC\u0083 [ dP0 d\u00CE\u00B7 ] 58 3.1. 1D Model in the Absence of Capillary Pressure and therefore we get for our three solvability conditions, [ \u00CE\u00BB\u00CC\u0084 dP1 d\u00CE\u00B7 ] = 0[ \u00CE\u00B1\u00CC\u0084s dP1 d\u00CE\u00B7 ] = s\u00E2\u0088\u0097\u00CF\u0083A [P1] = \u00E2\u0088\u0092A [ dP0 d\u00CE\u00B7 ] 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, \u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 s\u00E2\u0088\u0097 \u00CF\u0083 n \u00CE\u00B1\u00CC\u0084s(s \u00E2\u0088\u0097) \u00CE\u00B1\u00CC\u0084s(0)[ 1 \u00CE\u00BB\u00CC\u0084(0) \u00E2\u0088\u0092 1 \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097) ] 1 \u00E2\u0088\u00921 0 \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097) \u00CE\u00BB\u00CC\u0084(0) \u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB \u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 A C\u00E2\u0088\u0092 C+ \u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB \u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 0 0 0 \u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB and if we wish to have non-trivial eigenfunctions then the determinant of this matrix must be zero. This means that, \u00CF\u0083 n = ( \u00CE\u00BB\u00CC\u0084(0)\u00CE\u00B1\u00CC\u0084s(s \u00E2\u0088\u0097)\u00E2\u0088\u0092 \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097)\u00CE\u00B1\u00CC\u0084s(0) ) ( \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097)\u00E2\u0088\u0092 \u00CE\u00BB\u00CC\u0084(0)) s\u00E2\u0088\u0097\u00CE\u00BB\u00CC\u0084(0)\u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097)(\u00CE\u00BB\u00CC\u0084(0) + \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097)) which we can rearrange to \u00CF\u0083 n = v ( \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097)\u00E2\u0088\u0092 \u00CE\u00BB\u00CC\u0084(0))( \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097) + \u00CE\u00BB\u00CC\u0084(0) ) (3.15) by recalling that, f(s) = \u00CE\u00B1s \u00CE\u00BB 59 3.1. 1D Model in the Absence of Capillary Pressure and v = f(s\u00E2\u0088\u0097) s\u00E2\u0088\u0097 = f(s\u00E2\u0088\u0097)\u00E2\u0088\u0092 f(0) s\u00E2\u0088\u0097 \u00E2\u0088\u0092 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 \u00CE\u00B4 is sufficiently small then this will always be true. Consider the numerator for the eigenvalue, \u00CE\u00BB(s\u00E2\u0088\u0097)\u00E2\u0088\u0092 \u00CE\u00BB(0) = s\u00E2\u0088\u00973 + \u00CE\u00B4(1\u00E2\u0088\u0092 s\u00E2\u0088\u0097)3 \u00E2\u0088\u0092 \u00CE\u00B4. Recall from (3.6) that for \u00CE\u00B4 \u001C 1 the solution for the shock saturation value looks like, s\u00E2\u0088\u0097 \u00E2\u0088\u00BC 21/3\u00CE\u00B41/3 where we have taken only the first term which will provide us accuracy up to o(\u00CE\u00B4) but this is all that is required. With this in consideration we get that, s\u00E2\u0088\u00973 + \u00CE\u00B4(1\u00E2\u0088\u0092 s\u00E2\u0088\u0097)3 \u00E2\u0088\u0092 \u00CE\u00B4 \u00E2\u0089\u0088 2\u00CE\u00B4 > 0 and thus for large viscosity differences (small \u00CE\u00B4), fingering will always occur. As mentioned, this eigenvalue result is a generalized analysis of a front so- lution which transitions from s = 1 to s = 0 carried out in [47], however it yet generalizes even further when \u00CE\u00B1s(s \u00E2\u0088\u0097) = \u00CE\u00B1o(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 insta- bility occurs whenever \u00CE\u00B4 < 1; however in this problem we find values of \u00CE\u00B4 < 1 for which the system is stable. For example, consider \u00CE\u00B4 = 1/5, if we repeat 60 3.2. 1D Model With Capillary Pressure the analysis we get that s\u00E2\u0088\u0097 = 1/2. We then get that, \u00CE\u00BB(s\u00E2\u0088\u0097)\u00E2\u0088\u0092 \u00CE\u00BB(0) = \u00E2\u0088\u0092 1 20 < 0. Therefore we see in this quick example that there is a \u00CE\u00B4 < 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 permeabil- ities 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 \u00E2\u0088\u0092 Ps = \u000FJ(s) where J(s) is the Leverett function. As mentioned, it is chosen so that it is monotonically decreasing (J \u00E2\u0080\u00B2(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 \u00E2\u0080\u00B2(s) = \u00E2\u0088\u0092|J \u00E2\u0080\u00B2(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, qs = \u00CE\u00B1s \u00CE\u00BB \u00E2\u0088\u0092 \u000F\u00CE\u00B4\u00CE\u00B1s\u00CE\u00B1o \u00CE\u00BB |J \u00E2\u0080\u00B2(s)|\u00E2\u0088\u0082s \u00E2\u0088\u0082x , qo = \u00CE\u00B4\u00CE\u00B1o \u00CE\u00BB + \u000F \u00CE\u00B4\u00CE\u00B1o\u00CE\u00B1s \u00CE\u00BB |J \u00E2\u0080\u00B2(s)|\u00E2\u0088\u0082s \u00E2\u0088\u0082x . Now we expect that since \u000F \u001C 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, \u00CE\u00B7 = x\u00E2\u0088\u0092 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, \u00E2\u0088\u0082s \u00E2\u0088\u0082t + \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00CE\u00B7 (f(s)\u00E2\u0088\u0092 vs) = \u000F\u00CE\u00B4 \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00CE\u00B7 [(\u00CE\u00B1s\u00CE\u00B1o \u00CE\u00BB ) |J \u00E2\u0080\u00B2(s)|\u00E2\u0088\u0082s \u00E2\u0088\u0082\u00CE\u00B7 ] (3.16) of which we will seek steady-state solutions. Firstly we note that, for \u00CE\u00B7 of size O(1), if we were to expand s in powers of \u000F then to first order we just 62 3.2. 1D Model With Capillary Pressure get for the steady state, \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00CE\u00B7 (f(s)\u00E2\u0088\u0092 vs) = 0, (3.17a) \u00E2\u0088\u0082P \u00E2\u0088\u0082\u00CE\u00B7 = \u00E2\u0088\u00921 \u00CE\u00BB , (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 d\u00CE\u00B7 = 0 then if we continue the expansion we would see that the first order solution would persist at all orders of \u000F 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 \u00CE\u00B7 \u00E2\u0086\u0092 0 and so if we rescale our space variable, z = \u00CE\u00B7 \u000F , \u00E2\u0088\u0092\u00E2\u0088\u009E < z <\u00E2\u0088\u009E u(z) = s(\u000Fz) w(z) = P (\u000Fz) 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 dz (f(u)\u00E2\u0088\u0092 vu) = \u00CE\u00B4 d dz [(\u00CE\u00B1s\u00CE\u00B10 \u00CE\u00BB ) |J \u00E2\u0080\u00B2(u)|du 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 \u00E2\u0086\u0092\u00E2\u0088\u009E) to get a first order problem, du dz = ( \u00CE\u00BB \u00CE\u00B1s\u00CE\u00B1o )( f(u)\u00E2\u0088\u0092 vu \u00CE\u00B4|J \u00E2\u0080\u00B2(u)| ) . (3.18a) We may also define a steady-state pressure, dw dz = \u00E2\u0088\u0092\u000F ( 1 + \u00CE\u00B1s|J \u00E2\u0080\u00B2(u)|dudz ) \u00CE\u00BB 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 dz = \u000F vu\u00E2\u0088\u0092 1 \u00CE\u00B4\u00CE\u00B1o . (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 dz = \u00E2\u0088\u0092\u000Fvu \u00CE\u00B1s which diverges as s\u00E2\u0086\u0092 0. Since the oil pressure only diverges as s\u00E2\u0086\u0092 1, which is outside the domain of interest, the function is analytic where we are con- cerned 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 \u000F since \u00CE\u00B1s(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 \u00E2\u0080\u009Cwaiting-time\u00E2\u0080\u009D 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 cre- ates an infinite slope at the point where the travelling wave intersects s = 0. Physically this corresponds to the system attempting to preserve the hyper- bolicity 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 \u001C 1 (more specifically, we require u \u001C \u00CE\u00B41/3, a reasonable assumption since we are near u = 0) then, \u00CE\u00B1s(u)\u00CE\u00B1o(u) \u00CE\u00BB(u) \u00E2\u0088\u00BC u 3 \u00CE\u00B4 and f(u) \u00E2\u0088\u00BC u 3 \u00CE\u00B4 as well and therefore, du dz \u00E2\u0088\u00BC 1|J \u00E2\u0080\u00B2(0)| ( 1 \u00CE\u00B4 \u00E2\u0088\u0092 v u2 ) where we assume J(u) is smooth enough that J \u00E2\u0080\u00B2(0) is defined. Since u\u001C \u00CE\u00B41/3 then it is reasonable to assume that u2 \u001C \u00CE\u00B4 and so we can approximate 65 3.2. 1D Model With Capillary Pressure further, du dz \u00E2\u0089\u0088 \u00E2\u0088\u0092 v|J \u00E2\u0080\u00B2(0)| 1 u2 which has solution, u \u00E2\u0088\u00BC ( \u00E2\u0088\u0092 3v|J \u00E2\u0080\u00B2(0)|(z \u00E2\u0088\u0092 a) )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 \u00E2\u0080\u00B2(s)| = 1 \u00CE\u00B1s or some multiple there of then we still get for small u that, f(u) \u00E2\u0088\u00BC u 3 \u00CE\u00B4 but now, \u00CE\u00B1o(u) \u00CE\u00BB(u) \u00E2\u0088\u00BC 1 \u00CE\u00B4 so that when u\u001C \u00CE\u00B41/3 \u001C 1 that, du dz \u00E2\u0088\u00BC u3 \u00E2\u0088\u0092 \u00CE\u00B4vu. 6While 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\u001C \u00CE\u00B41/3 then it is reasonable to assume that u3 \u001C \u00CE\u00B4u and so we can approximate further to finally get that, du dz \u00E2\u0088\u00BC \u00E2\u0088\u0092\u00CE\u00B4vu so that without degeneracy we have the standard linear decay rate which corresponds to exponential decay, u(z) \u00E2\u0088\u00BC exp (\u00E2\u0088\u0092\u00CE\u00B4vz) near u = 0. For completeness we will look how the solution decays to the other constant state, u = s\u00E2\u0088\u0097. In this case, f(u)\u00E2\u0088\u0092 vu \u00E2\u0089\u0088 (f(s\u00E2\u0088\u0097)\u00E2\u0088\u0092 vs\u00E2\u0088\u0097)\u00EF\u00B8\u00B8 \u00EF\u00B8\u00B7\u00EF\u00B8\u00B7 \u00EF\u00B8\u00B8 =0 + (f \u00E2\u0080\u00B2(s\u00E2\u0088\u0097)\u00E2\u0088\u0092 v)\u00EF\u00B8\u00B8 \u00EF\u00B8\u00B7\u00EF\u00B8\u00B7 \u00EF\u00B8\u00B8 =0 (u\u00E2\u0088\u0092 s\u00E2\u0088\u0097) + 1 2 f \u00E2\u0080\u00B2\u00E2\u0080\u00B2(s\u00E2\u0088\u0097) (u\u00E2\u0088\u0092 s\u00E2\u0088\u0097)2 = 1 2 f \u00E2\u0080\u00B2\u00E2\u0080\u00B2(s\u00E2\u0088\u0097) (u\u00E2\u0088\u0092 s\u00E2\u0088\u0097)2 so therefore, \u00E2\u0088\u0082u \u00E2\u0088\u0082z \u00E2\u0088\u00BC \u00CE\u00BB(s \u00E2\u0088\u0097) \u00CE\u00B1s(s\u00E2\u0088\u0097)\u00CE\u00B1o(s\u00E2\u0088\u0097) f \u00E2\u0080\u00B2\u00E2\u0080\u00B2(s\u00E2\u0088\u0097) 2\u00CE\u00B4|J \u00E2\u0080\u00B2(s\u00E2\u0088\u0097)| (u\u00E2\u0088\u0092 s \u00E2\u0088\u0097)2 for u \u00E2\u0089\u0088 s\u00E2\u0088\u0097 where once again we require J(u) smooth enough so that J \u00E2\u0080\u00B2(s\u00E2\u0088\u0097) 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., lim \u00CE\u00B7\u00E2\u0086\u00920 ds d\u00CE\u00B7 = 0 = lim z\u00E2\u0086\u0092\u00C2\u00B1\u00E2\u0088\u009E du dz lim \u00CE\u00B7\u00E2\u0086\u00920 dP d\u00CE\u00B7 = \u000F ( \u00E2\u0088\u00921 \u00CE\u00BB ) = lim z\u00E2\u0086\u0092\u00C2\u00B1\u00E2\u0088\u009E dw dz . 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. -10 -8 -6 -4 -2 0 2 4 6 8 10 -0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Saturation Steady-State:Boundary Layer z u 0 Figure 3.3: Saturation steady-state solution for the inner boundary layer centered around \u00CE\u00B7 = 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 \u00E2\u0088\u009210 \u00E2\u0088\u00928 \u00E2\u0088\u00926 \u00E2\u0088\u00924 \u00E2\u0088\u00922 0 2 4 6 8 10 \u00E2\u0088\u00925.5 \u00E2\u0088\u00925 \u00E2\u0088\u00924.5 \u00E2\u0088\u00924 \u00E2\u0088\u00923.5 \u00E2\u0088\u00923 \u00E2\u0088\u00922.5 \u00E2\u0088\u00922 Pressure Derivative Steady\u00E2\u0088\u0092State:Boundary Layer z d w 0 d z Figure 3.4: Pressure derivative steady-state solution for the inner boundary layer centered around \u00CE\u00B7 = 0. The derivative has the appropriate limiting behaviour to the outer hyperbolic problem. The saturation is solved numerically using ODE45 in MATLAB and the pres- sure 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 \u00E2\u0088\u00BC s0 + \u00CE\u00B5s1(\u00CE\u00B7) exp (iny + \u00CF\u0083t) , \u00CE\u00B7 6= 0 P \u00E2\u0088\u00BC P0 + \u00CE\u00B5P1(\u00CE\u00B7) exp (iny + \u00CF\u0083t) , \u00CE\u00B7 6= 0 s1, ds1 d\u00CE\u00B7 , P1, dP1 d\u00CE\u00B7 \u00E2\u0086\u0092 0, |\u00CE\u00B7| \u00E2\u0086\u0092 \u00E2\u0088\u009E where as before, the conditions are chosen to have smooth decay to the base state for large |\u00CE\u00B7| and \u00CE\u00B5 is used to control term balance. We now wish to substitute this into the full PDE (2.9), \u00E2\u0088\u0082s \u00E2\u0088\u0082t \u00E2\u0088\u0092\u00E2\u0088\u0087 \u00C2\u00B7 (qo + vsx\u00CC\u0082) = 0 qs = \u00E2\u0088\u0092\u00CE\u00B1s\u00E2\u0088\u0087(P \u00E2\u0088\u0092 \u000FJ(s)) qo = \u00E2\u0088\u0092\u00CE\u00B4\u00CE\u00B1o\u00E2\u0088\u0087P \u00E2\u0088\u0087 \u00C2\u00B7 (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 \u00CE\u00B5, d d\u00CE\u00B7 ( \u00CE\u00BB\u00CC\u0084 dP1 d\u00CE\u00B7 + \u00CE\u00BB\u00CC\u0084\u00E2\u0080\u00B2 dP0 d\u00CE\u00B7 s1 ) + \u000F d2 d\u00CE\u00B72 (\u00CE\u00B1\u00CC\u0084s|J \u00E2\u0080\u00B2(s0)|s1) \u00E2\u0088\u0092n2 (\u00CE\u00BB\u00CC\u0084P1 + \u000F\u00CE\u00B1\u00CC\u0084s|J \u00E2\u0080\u00B2(s0)|s1) = 0 (3.19a) d d\u00CE\u00B7 ( \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o dP1 d\u00CE\u00B7 + \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o \u00E2\u0080\u00B2dP0 d\u00CE\u00B7 s1 \u00E2\u0088\u0092 vs1 ) \u00E2\u0088\u0092 n2\u00CE\u00B4\u00CE\u00B1\u00CC\u0084oP1 = \u00E2\u0088\u0092\u00CF\u0083s1 (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 \u000F = 0, with the use of the identity \u00CE\u00B4\u00CE\u00B1o = \u00CE\u00BB\u00E2\u0088\u0092 \u00CE\u00B1s, 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 \u00CE\u00B7 to be O(1) then in this region we take for the steady state (3.17), ds0 d\u00CE\u00B7 = 0 dP0 d\u00CE\u00B7 = \u00E2\u0088\u00921 \u00CE\u00BB\u00CC\u0084 and if we proceed in an expansion sense for s1, P1, and \u00CF\u0083 as, s1 \u00E2\u0088\u00BC s10 + \u000Fs11 + . . . P1 \u00E2\u0088\u00BC P10 + \u000FP11 + . . . \u00CF\u0083 \u00E2\u0088\u00BC \u00CF\u00830 + \u000F\u00CF\u00831 + . . . 71 3.2. 1D Model With Capillary Pressure then we get to first order the eigenvalue problem for \u000F = 0 and thus, s10 = 0 P10 = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3C0+ exp (\u00E2\u0088\u0092n\u00CE\u00B7) , \u00CE\u00B7 > 0C0\u00E2\u0088\u0092 exp (n\u00CE\u00B7) , \u00CE\u00B7 < 0 . Now in the \u000F = 0 case we proceeded to find the eigenvalue via jump conditions across the shock interface. However, with \u000F 6= 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(\u000F) 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 \u000F = 0 holds at all orders. Thus we get for the outer solution, s1 = 0 (3.20a) P1 = ( \u00E2\u0088\u009E\u00E2\u0088\u0091 p=0 \u000FpCpk ) exp (\u00E2\u0088\u0092kn\u00CE\u00B7) (3.20b) where k = + or k = \u00E2\u0088\u0092 to differentiate between the solutions for \u00CE\u00B7 > 0 and \u00CE\u00B7 < 0. Inner Solution, Non-Degenerate Case Like with the steady-state solution we are motivated to pose an inner vari- able 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 \u00CE\u00B7 = 0 and let, z = \u00CE\u00B7 \u000F , \u00E2\u0088\u0092\u00E2\u0088\u009E < z <\u00E2\u0088\u009E u(z) = s(\u000Fz) w(z) = P (\u000Fz). We consider the problem for the non-degenerate steady-state solution so that we have analyticity throughout the domain. We therefore will assume that, \u00CE\u00B1\u00CC\u0084s|J \u00E2\u0080\u00B2(u0)| = 1 and we take as our steady-state solution in this region, (3.18), du0 dz = \u00CE\u00BB\u00CC\u0084 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o (f(u0)\u00E2\u0088\u0092 vu0) dw0 dz = \u000F vu0 \u00E2\u0088\u0092 1 \u00CE\u00B4\u00CE\u00B1o = \u000F dw\u00CC\u00830 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 ( \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o dw1 dz + \u000F\u00CE\u00B4\u00CE\u00B1\u00CC\u0084o \u00E2\u0080\u00B2dw\u00CC\u00830 dz u1 \u00E2\u0088\u0092 \u000Fvu1 ) \u00E2\u0088\u0092 \u000F2n2\u00CE\u00B4\u00CE\u00B1\u00CC\u0084ow1 = \u00E2\u0088\u0092\u000F2u1\u00CF\u0083 d dz ( \u00CE\u00BB\u00CC\u0084 dw1 dz + \u000F\u00CE\u00BB\u00CC\u0084\u00E2\u0080\u00B2 dw\u00CC\u00830 dz u1 ) + d2 dz2 (\u000Fu1) \u00E2\u0088\u0092\u000F2n2 (\u00CE\u00BB\u00CC\u0084w1 + \u000Fu1) = 0 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\u00CC\u00821 \u000F so that we can write as our final inner variable eigenvalue problem, d dz ( \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o dw1 dz + \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o \u00E2\u0080\u00B2dw\u00CC\u00830 dz u\u00CC\u00821 \u00E2\u0088\u0092 vu\u00CC\u00821 ) \u00E2\u0088\u0092 \u000F2n2\u00CE\u00B4\u00CE\u00B1\u00CC\u0084ow1 = \u00E2\u0088\u0092\u000Fu\u00CC\u00821\u00CF\u0083 (3.21a) d dz ( \u00CE\u00BB\u00CC\u0084 dw1 dz + \u00CE\u00BB\u00CC\u0084\u00E2\u0080\u00B2 dw\u00CC\u00830 dz u\u00CC\u00821 ) + d2 dz2 (u\u00CC\u00821) \u00E2\u0088\u0092\u000F2n2 (\u00CE\u00BB\u00CC\u0084w1 + u\u00CC\u00821) = 0. (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, \u00CE\u00B4(x) = lim \u000F\u00E2\u0086\u00920 1 \u000F h (x \u000F ) 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 foun- dation for which to compare any analytic results as well as to have a qualita- tive 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 val- ues, we place all information relating to u\u00CC\u00821 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 ( \u00CE\u00BB\u00CC\u0084 dw1 dz ) we denote, J = \u00CE\u00BB\u00CC\u0084 dw1 dz as a flux term and place it on the nodes as in Figure 3.5. 75 3.2. 1D Model With Capillary Pressure \u00CE\u00B1\u00CC\u0084oi+3/2 w1(i+1/2) \u00CE\u00BB\u00CC\u0084i+1/2 \u00CE\u00B1\u00CC\u0084oi+1/2 u\u00CC\u00821(i\u00E2\u0088\u00921/2) w1(i\u00E2\u0088\u00921/2) \u00CE\u00BB\u00CC\u0084i\u00E2\u0088\u00921/2 \u00CE\u00B1\u00CC\u0084oi\u00E2\u0088\u00921/2 xi Ji xi+1 Ji+1 u\u00CC\u00821(i+3/2) w1(i+3/2) \u00CE\u00BB\u00CC\u0084i+3/2 u\u00CC\u00821(i+1/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 discretiza- tion 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, dJ dz = Ji+1 \u00E2\u0088\u0092 Ji 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 dif- ference operator on the nodes so that the two flux terms share a cell-center. By definition, Ji = \u00CE\u00BB\u00CC\u0084i dw1(i) dz but since \u00CE\u00BB and w1 are function of u\u00CC\u00821 and w1 then they are only defined on the cell-centers and thus we must extrapolate them to the nodes. For \u00CE\u00BB\u00CC\u0084 we 76 3.2. 1D Model With Capillary Pressure simply take an average, \u00CE\u00BB\u00CC\u0084i = \u00CE\u00BB\u00CC\u0084i\u00E2\u0088\u00921/2 + \u00CE\u00BB\u00CC\u0084i+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) dz = w1(i+1/2) \u00E2\u0088\u0092 w1(i\u00E2\u0088\u00921/2) 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 dz ( \u00CE\u00BB\u00CC\u0084 dw1 dz ) = (\u00CE\u00BB\u00CC\u0084i+3/2 + \u00CE\u00BB\u00CC\u0084i+1/2)(w1(i+3/2) \u00E2\u0088\u0092 w1(i+1/2))\u00E2\u0088\u0092 (\u00CE\u00BB\u00CC\u0084i+1/2 + \u00CE\u00BB\u00CC\u0084i\u00E2\u0088\u00921/2)(w1(i+1/2) \u00E2\u0088\u0092 w1(i\u00E2\u0088\u00921/2)) 2h2 . Notice if \u00CE\u00BB\u00CC\u0084 = 1 we recover the standard second-order centered operator, d2w1 dz2 = w1(i+3/2) \u00E2\u0088\u0092 2w1(i+1/2) + w1(1\u00E2\u0088\u00921/2) h2 . We could discretize the other differential operators in a similar fashion to generate a generalized numerical eigenvalue problem, Au = \u00CF\u0083Bu 77 3.2. 1D Model With Capillary Pressure where the vectors u are, u = [u\u00CC\u00821,w1] T . Boundary Conditions Now in order to solve the eigenvalue problem we require homogeneous bound- ary conditions so that we do not add any inhomogeneous terms to our right- hand 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 dz = \u00E2\u0088\u0092\u000Fnw1, z \u00E2\u0086\u0092\u00E2\u0088\u009E dw1 dz = \u000Fnw1, z \u00E2\u0086\u0092 \u00E2\u0088\u0092\u00E2\u0088\u009E. Notice that since the domain endpoints are nodal points, our far-field con- ditions 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\u00E2\u0088\u00921/2 to the left of the first node x1 then the far-field condition tells us that, w1(1\u00E2\u0088\u00921/2) = ( 1\u00E2\u0088\u0092 \u000Fnh 2 1 + \u000Fnh 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 Nth node, xN+1/2 and the far field condition 78 3.2. 1D Model With Capillary Pressure imposes that, w1(N+1/2) = ( 1\u00E2\u0088\u0092 \u000Fnh 2 1 + \u000Fnh 2 ) w1(N\u00E2\u0088\u00921/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\u00CC\u00821(1+1/2) = 0 u\u00CC\u00821(N\u00E2\u0088\u00921/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 \u000F = 0 case, when \u00CF\u0083 > 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. 7While technically it is necessary to also use the ghost point ideology for this eigenfunc- tion 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 \u00E2\u0088\u009210 \u00E2\u0088\u00928 \u00E2\u0088\u00926 \u00E2\u0088\u00924 \u00E2\u0088\u00922 0 2 4 6 8 10 \u00E2\u0088\u00920.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 z u\u00CC\u0082 1 Saturation Eigenfunction for n=1 Figure 3.6: Numerical solution of u\u00CC\u00821 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 -10 -8 -6 -4 -2 0 2 4 6 8 10 -600 -400 -200 0 200 400 600 800 1000 1200 1400 z w 1 Pressure Eigenfunction for n=1 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 \u00CF\u0083 = 1.81. Figure 3.8 plots the value of \u00CF\u0083 for various values of n. 81 3.2. 1D Model With Capillary Pressure 0 20 40 60 80 100 120 140 160 180 0 10 20 30 40 50 60 \u00CF\u0083 n Largest Eigenvalue vs. n Figure 3.8: This figure shows the largest value of \u00CF\u0083 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 nu- merical problem. While the properties of \u00CF\u0083 in the presence of capillary effects has been noted, often in the literature ([34]) a highly robust scheme is re- quired 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 \u000F is and thus is handled quite easily through standard techniques. We conducted a numerical convergence study to test robustness by first letting h \u00E2\u0086\u0092 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 \u00CF\u0083 and thus we proceed in the usual sense by expanding the eigenfunctions and eigenvalue of (3.21) as, u\u00CC\u00821 \u00E2\u0088\u00BC u\u00CC\u008210 + \u000Fu\u00CC\u008211 + . . . w1 \u00E2\u0088\u00BC w10 + \u000Fw11 + . . . \u00CF\u0083 \u00E2\u0088\u00BC \u00CF\u00830 + . . . . To first order we get that, d dz ( \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o dw10 dz + \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o \u00E2\u0080\u00B2dw\u00CC\u008310 dz u\u00CC\u008210 \u00E2\u0088\u0092 vu\u00CC\u008210 ) = 0 (3.22a) d dz ( \u00CE\u00BB\u00CC\u0084 dw10 dz + \u00CE\u00BB\u00CC\u0084\u00E2\u0080\u00B2 dw\u00CC\u00830 dz u\u00CC\u008210 ) + d2 dz2 (u\u00CC\u008210) = 0. (3.22b) This is valid as long as n is not sufficiently large (more specifically as long as n is not O(1/\u000F)). 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 P1 = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3 \u00E2\u0088\u009E\u00E2\u0088\u0091 p=0 \u000FpCp+ exp (\u00E2\u0088\u0092n\u00CE\u00B7) , \u00CE\u00B7 > 0 \u00E2\u0088\u009E\u00E2\u0088\u0091 p=0 \u000FpCp\u00E2\u0088\u0092 exp (n\u00CE\u00B7) , \u00CE\u00B7 < 0 and so if we Taylor expand around \u00CE\u00B7 = 0 we generate the following matching conditions, u\u00CC\u008210 \u00E2\u0086\u0092 0, |z| \u00E2\u0086\u0092 \u00E2\u0088\u009E (3.23a) w10 \u00E2\u0086\u0092 C0+, z \u00E2\u0086\u0092\u00E2\u0088\u009E (3.23b) w10 \u00E2\u0086\u0092 C0\u00E2\u0088\u0092, z \u00E2\u0086\u0092 \u00E2\u0088\u0092\u00E2\u0088\u009E. (3.23c) Since the outer solution for the saturation eigenfunction vanishes smoothly, we assume that the derivatives of u\u00CC\u008210 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 con- ditions we get, ( \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o dw10 dz )\u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E = 0. (3.24a) If we integrate (3.22b) and use the matching conditions we get that, ( \u00CE\u00BB\u00CC\u0084 dw10 dz )\u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E = 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, \u00CE\u00BB\u00CC\u0084 dw10 dz + \u00CE\u00BB\u00CC\u0084\u00E2\u0080\u00B2 dw\u00CC\u00830 dz u\u00CC\u008210 + du\u00CC\u00821 dz = B0 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o dw10 dz + \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o \u00E2\u0080\u00B2dw\u00CC\u00830 dz u\u00CC\u008210 \u00E2\u0088\u0092 vu\u00CC\u008210 = D0 where B0 and D0 are constants of integration. We can rearrange the second equation to get, dw10 dz = ( v \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o \u00E2\u0088\u0092 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o \u00E2\u0080\u00B2 (\u00CE\u00B4\u00CE\u00B1\u00CC\u0084o)2 (vu0 \u00E2\u0088\u0092 1) ) u\u00CC\u008210 + D0 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o = 1 du0 dz d dz ( vu0 \u00E2\u0088\u0092 1 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o ) u\u00CC\u008210 + D0 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o = 1 du0 dz d2w\u00CC\u00830 dz2 u\u00CC\u008210 + D0 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o (3.25) where we have used the knowledge that for any g(u(z)), g\u00E2\u0080\u00B2 = dg du0 = 1 du0 dz dg dz . Notice firstly we can integrate this expression to get another solvability con- dition, w10|\u00E2\u0088\u009E\u00E2\u0088\u0092\u00E2\u0088\u009E = \u00E2\u0088\u00AB \u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E ( 1 du0 dz d2w\u00CC\u00830 dz2 u\u00CC\u008210 + D0 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o ) dz (3.26) 85 3.2. 1D Model With Capillary Pressure but also that we can solve solely for the saturation eigenfunction, 1 du0 dz ( \u00CE\u00BB\u00CC\u0084 d2w\u00CC\u00830 dz2 + d\u00CE\u00BB\u00CC\u0084 dz dw\u00CC\u00830 dz ) u\u00CC\u008210 + du\u00CC\u008210 dz = B0 \u00E2\u0088\u0092 \u00CE\u00BB\u00CC\u0084 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o D0. (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, du\u00CC\u008210 dz = \u00E2\u0088\u0092 1 du0 dz d dz ( \u00CE\u00BB\u00CC\u0084 dw\u00CC\u00830 dz ) u\u00CC\u008210. In order to proceed further we notice that we can relate the base state deriva- tives via, du0 dz = \u00E2\u0088\u0092 ( \u00CE\u00BB\u00CC\u0084 dw\u00CC\u00830 dz + 1 ) so that we get, du\u00CC\u008210 dz = 1 du0 dz d2u0 dz2 = d dz ln ( du0 dz ) 8It may seem audacious to state this result since u\u00CC\u008210 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\u00CC\u008210 = A0 du0 dz . (3.28) We can substitute this back into (3.25) to solve for w10, w10 = A0 dw\u00CC\u00830 dz + E0 (3.29) where we can use the matching conditions to get that the constant is, E0 = ( 1\u00E2\u0088\u0092 \u00CE\u00BB\u00CC\u0084(0) \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097) )\u00E2\u0088\u00921( C\u00E2\u0088\u0092 \u00E2\u0088\u0092 \u00CE\u00BB\u00CC\u0084(0) \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097) C+ ) . We can now simplify solvability conditions (3.24) and (3.26), ( \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o dw10 dz )\u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E = 0 (3.30a)( \u00CE\u00BB\u00CC\u0084 dw10 dz )\u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E = 0 (3.30b) w10|\u00E2\u0088\u009E\u00E2\u0088\u0092\u00E2\u0088\u009E = A0 dw\u00CC\u00830 dz \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E (3.30c) which are indeed satisfied with our solution for any constants A0, C0+ and C0\u00E2\u0088\u0092. 87 3.2. 1D Model With Capillary Pressure Determining \u00CF\u00830 We are now able to look at the next order, O(\u000F) to get that our eigenvalue problem is, d dz ( \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o dw11 dz + \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o \u00E2\u0080\u00B2dw\u00CC\u00830 dz u\u00CC\u008211 \u00E2\u0088\u0092 vu\u00CC\u008211 ) = \u00E2\u0088\u0092\u00CF\u00830u\u00CC\u008210 (3.31a) d dz ( \u00CE\u00BB\u00CC\u0084 dw11 dz + \u00CE\u00BB\u00CC\u0084\u00E2\u0080\u00B2 dw\u00CC\u00830 dz u\u00CC\u008211 ) + d2 dz2 (u\u00CC\u008211) = 0. (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\u00CC\u008211 \u00E2\u0086\u0092 0, |z| \u00E2\u0086\u0092 \u00E2\u0088\u009E (3.32a) w11 \u00E2\u0086\u0092 C1+ \u00E2\u0088\u0092 nC0+z, z \u00E2\u0086\u0092\u00E2\u0088\u009E (3.32b) w11 \u00E2\u0086\u0092 C1\u00E2\u0088\u0092 + nC0\u00E2\u0088\u0092z, z \u00E2\u0086\u0092 \u00E2\u0088\u0092\u00E2\u0088\u009E. (3.32c) However, in what follows we will set C1+ = C1\u00E2\u0088\u0092 = 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), ( \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o dw11 dz )\u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E = \u00E2\u0088\u0092\u00CF\u00830 \u00E2\u0088\u00AB \u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E u\u00CC\u008210 dz = A0\u00CF\u00830s \u00E2\u0088\u0097 (3.33a)( \u00CE\u00BB\u00CC\u0084 dw11 dz )\u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E = 0. (3.33b) 88 3.2. 1D Model With Capillary Pressure To solve (3.31) we take a first integral to get, \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o dw11 dz + \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o \u00E2\u0080\u00B2dw\u00CC\u00830 dz u\u00CC\u008211 \u00E2\u0088\u0092 vu\u00CC\u008211 = \u00E2\u0088\u0092A0\u00CF\u00830u0 +D1 \u00CE\u00BB\u00CC\u0084 dw11 dz + \u00CE\u00BB\u00CC\u0084\u00E2\u0080\u00B2 dw\u00CC\u00830 dz u\u00CC\u008211 + du\u00CC\u008211 dz = B1 where B1 and D1 are constants of integration. We can rearrange the first equation to get, dw11 dz = 1 du0 dz d2w\u00CC\u00830 dz2 u\u00CC\u008211 \u00E2\u0088\u0092 A0\u00CF\u00830 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o u0 + D1 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o (3.34) which we can use to solve entirely for the saturation eigenfunction, du\u00CC\u008211 dz \u00E2\u0088\u0092 1 du0 dz ( d2u0 dz2 ) u\u00CC\u008211 = B1 + \u00CE\u00BB\u00CC\u0084A0\u00CF\u00830 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o u0 \u00E2\u0088\u0092 \u00CE\u00BB\u00CC\u0084 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o D1. (3.35) Notice in both of these simplifications we have recycled ideas from the pre- vious section regarding relationships between steady states. Now when we apply the matching conditions we get that the constants must satisfy, B1 = D1 = A0\u00CF\u00830 v where we have used the fact that, s\u00E2\u0088\u0097 f(s\u00E2\u0088\u0097) = 1 v . 89 3.2. 1D Model With Capillary Pressure We can then write the solution for (3.35) as, u\u00CC\u008211 = A1 du0 dz + A0\u00CF\u00830 v du0 dz \u00E2\u0088\u00AB ( 1 + \u00CE\u00BB\u00CC\u0084 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o (vu0 \u00E2\u0088\u0092 1) ) 1 du0 dz dz. However, 1 + \u00CE\u00BB\u00CC\u0084 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o (vu0 \u00E2\u0088\u0092 1) = \u00E2\u0088\u0092du0 dz and so finally we get that, u\u00CC\u008211 = A1 du0 dz \u00E2\u0088\u0092 A0\u00CF\u00830 v du0 dz z. (3.36) We can substitute this into (3.34) to get that, dw11 dz = A1 d2w\u00CC\u00830 dz2 \u00E2\u0088\u0092 A0\u00CF\u00830 v ( z d2w\u00CC\u00830 dz2 + vu0 \u00E2\u0088\u0092 1 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o ) = A1 d2w\u00CC\u00830 dz2 \u00E2\u0088\u0092 A0\u00CF\u00830 v d dz ( z dw\u00CC\u00830 dz ) where we have used that, vu0 \u00E2\u0088\u0092 1 \u00CE\u00B4\u00CE\u00B1\u00CC\u00840 = dw\u00CC\u00830 dz . Finally then we get that, w11 = A1 dw\u00CC\u00830 dz \u00E2\u0088\u0092 A0\u00CF\u00830 v dw\u00CC\u00830 dz z + E1 (3.37) 90 3.2. 1D Model With Capillary Pressure and after using the matching conditions (3.32), we conclude that9, E1 = 0 A0\u00CF\u00830 v = \u00E2\u0088\u0092\u00CE\u00BB\u00CC\u0084(0)nC+ = \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097)nC\u00E2\u0088\u0092. 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, ( \u00CE\u00BB\u00CC\u0084 dw11 dz )\u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E = 0 which after using the matching conditions (3.32) becomes, \u00E2\u0088\u0092\u00CE\u00BB\u00CC\u0084(0)nC+ = \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097)nC\u00E2\u0088\u0092, precisely the condition we need. Having the solutions u\u00CC\u008211 and w11 we can write the complete solvability conditions at this order as, ( \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o dw11 dz )\u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E = A0\u00CF\u00830s \u00E2\u0088\u0097 (3.38a)( \u00CE\u00BB\u00CC\u0084 dw11 dz )\u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E = 0 (3.38b) (w11) \u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E = 0. (3.38c) 9This 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(\u000F) as, u\u00CC\u00821 = A0 du0 dz + \u000FA1 du0 dz \u00E2\u0088\u0092 \u000FA0\u00CF\u00830 v du0 dz z (3.39a) w1 = A0 dw\u00CC\u00830 dz + E0 \u00E2\u0088\u0092 \u000FA0\u00CF\u00830 v dw\u00CC\u00830 dz z (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, ( \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o dw1 dz )\u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E = \u000FA0\u00CF\u00830s \u00E2\u0088\u0097 ( \u00CE\u00BB\u00CC\u0084 dw1 dz )\u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E = 0 (w1) \u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E = A0 ( dw\u00CC\u00830 dz )\u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E . If we then scale back to the outer variables, let A\u00CC\u00830 = \u00E2\u0088\u0092A0, and use the identity \u00CE\u00BB = \u00CE\u00B1s + \u00CE\u00B4\u00CE\u00B1o, we get,[ \u00CE\u00BB\u00CC\u0084 dP1 d\u00CE\u00B7 ] = 0 (3.40a)[ \u00CE\u00B1\u00CC\u0084s dP1 d\u00CE\u00B7 ] = A\u00CC\u00830\u00CF\u00830s \u00E2\u0088\u0097 (3.40b) [P1] = \u00E2\u0088\u0092A\u00CC\u00830 [ dP0 d\u00CE\u00B7 ] (3.40c) where [\u00C2\u00B7] indicates a jump across \u00CE\u00B7 = 0 like before. These are in fact the same jump conditions (3.14) for the \u000F = 0 problem and thus we conclude 92 3.2. 1D Model With Capillary Pressure that, \u00CF\u00830 = nv \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097)\u00E2\u0088\u0092 \u00CE\u00BB\u00CC\u0084(0) \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097) + \u00CE\u00BB\u00CC\u0084(0) (3.41) as we would have expected from the onset. Knowing \u00CF\u00830 we can use (3.40) to get that for any A\u00CC\u00830, we require, C+ = \u00E2\u0088\u0092 A\u00CC\u00830 \u00CE\u00BB\u00CC\u0084(0) ( 1 + \u00CE\u00BB\u00CC\u0084(0) \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097) )\u00E2\u0088\u00921( 1\u00E2\u0088\u0092 \u00CE\u00BB\u00CC\u0084(0) \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097) ) C\u00E2\u0088\u0092 = A\u00CC\u00830 \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097) ( 1 + \u00CE\u00BB\u00CC\u0084(0) \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097) )\u00E2\u0088\u00921( 1\u00E2\u0088\u0092 \u00CE\u00BB\u00CC\u0084(0) \u00CE\u00BB\u00CC\u0084(s\u00E2\u0088\u0097) ) . The fact that we are still left over with an arbitrary constant is not unsur- prising 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\u00CC\u00821(0), we instead impose a normalizing condition, \u00E2\u0088\u00AB \u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E u\u00CC\u00821 dz = 1. This supports the eigenfunction approaching a Dirac-delta function as \u000F\u00E2\u0086\u0092 0. If we impose such a normalization condition we get that, \u00E2\u0088\u00AB \u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E u\u00CC\u008210 dz = \u00E2\u0088\u0092A0s\u00E2\u0088\u0097 = 1 and so, A0 = \u00E2\u0088\u0092 1 s\u00E2\u0088\u0097 . 93 3.2. 1D Model With Capillary Pressure At the next order then we get that, \u00E2\u0088\u00AB \u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E u\u00CC\u008211 dz = 0 which requires that, \u00E2\u0088\u00AB \u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E ( A1 du0 dz \u00E2\u0088\u0092 zdu0 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 \u00E2\u0088\u00BC \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3s \u00E2\u0088\u0097z, 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\u00E2\u0088\u0092 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 satu- ration eigenfunction obeys our solvability condition, we compute it and then numerically integrate over the domain. We then divide our entire eigenvec- tor by this result and this ensures that the saturation eigenfunction has the required property that, \u00E2\u0088\u00AB \u00E2\u0088\u009E \u00E2\u0088\u0092\u00E2\u0088\u009E u\u00CC\u00821 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(\u000F). 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, auto- matically, 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 \u000F 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\u00CC\u008210 and u\u00CC\u008211 computed asymptotically. Figure 3.10 shows the same result but for the pressure saturation eigenfunction. \u00E2\u0088\u009210 \u00E2\u0088\u00928 \u00E2\u0088\u00926 \u00E2\u0088\u00924 \u00E2\u0088\u00922 0 2 4 6 8 10 \u00E2\u0088\u00920.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 z u\u00CC\u0082 1 Saturation Eigenfunction for n=1 Numeric u10 u10+\u00CE\u00B5u11 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 \u00E2\u0088\u009210 \u00E2\u0088\u00928 \u00E2\u0088\u00926 \u00E2\u0088\u00924 \u00E2\u0088\u00922 0 2 4 6 8 10 \u00E2\u0088\u00921000 \u00E2\u0088\u0092500 0 500 1000 1500 z w 1 Pressure Eigenfunction for n=1 Numeric w10 w10+\u00CE\u00B5w11 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 \u00CF\u00830 = 1.81 while asymptotically we estimate that \u00CF\u00830 = 1.85, a good agreement. Large n solution If n is large, i.e., O(1/\u000F) then terms previously neglected in (3.21) are no longer inadmissible. With this in mind we scale, n\u00CC\u0082 = \u000Fn \u00CF\u0083\u00CC\u0082 = \u000F\u00CF\u0083 96 3.2. 1D Model With Capillary Pressure and we get a completely O(1) problem, free of \u000F, d dz ( \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o dw1 dz + \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o \u00E2\u0080\u00B2dw\u00CC\u00830 dz u\u00CC\u00821 \u00E2\u0088\u0092 vu\u00CC\u00821 ) \u00E2\u0088\u0092 n\u00CC\u00822\u00CE\u00B4\u00CE\u00B1\u00CC\u0084ow1 = \u00E2\u0088\u0092u\u00CC\u00821\u00CF\u0083\u00CC\u0082 (3.42a) d dz ( \u00CE\u00BB\u00CC\u0084 dw1 dz + \u00CE\u00BB\u00CC\u0084\u00E2\u0080\u00B2 dw\u00CC\u00830 dz u\u00CC\u00821 ) + d2 dz2 (u\u00CC\u00821) \u00E2\u0088\u0092n\u00CC\u00822 (\u00CE\u00BB\u00CC\u0084w1 + u\u00CC\u00821) = 0. (3.42b) We can use our finite difference solver and compute the largest eigenvalue and eigenfunction for each n\u00CC\u0082. This is plotted in Figure 3.11. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 0 0.1 0.2 0.3 0.4 0.5 0.6 \u00CF\u0083\u00CC\u0082 n\u00CC\u0082 Largest Eigenvalue vs. n\u00CC\u0082 Figure 3.11: Eigenvalue plot for various values of n\u00CC\u0082 = \u000Fn. 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 \u000F. This formulation extends that idea fur- ther; by considering a strictly O(1) problem; \u000F has truly been removed and the numerical methods will always be solving a smooth problem. Therefore Figure 3.11 holds for all values of \u000F but of course, for fixed n, n\u00CC\u0082 \u00E2\u0086\u0092 0 as \u000F \u00E2\u0086\u0092 0. This means that as \u000F \u00E2\u0086\u0092 0, the linear eigenvalue persists for larger values of n as expected. The universality of this formulation being indepen- dent of \u000F 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\u00CC\u0082, 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 \u000F\u00E2\u0086\u0092\u00E2\u0088\u009E, then all wave-numbers produce stable fingering patterns, an interesting result although not entirely unsurprising since \u000F \u00E2\u0086\u0092 \u00E2\u0088\u009E is equivalent to large surface tension which is a strong stabilizing mechanism. There is a caveat associated with this formulation. If the remaining asymp- totic terms to (3.21) are algebraic with respect to n, i.e. they can be written in the form, \u00CF\u0083 \u00E2\u0088\u00BC \u00E2\u0088\u009E\u00E2\u0088\u0091 p=0 \u000Fp\u00CF\u0083pn p+1 and if we rescale using the large n formulation as, n\u00CC\u0082 = \u000Fn \u00CF\u0083\u00CC\u0082 = \u000F\u00CF\u0083 98 3.3. 1D Temperature Model then, \u00CF\u0083\u00CC\u0082 \u00E2\u0088\u00BC \u00E2\u0088\u009E\u00E2\u0088\u0091 p=0 \u00CF\u0083pn p+1. 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 \u00CF\u0083 in the large n case nor vice- versa 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 \u00CF\u0083 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 per- sisted. 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), \u00E2\u0088\u0082T \u00E2\u0088\u0082t + (CsC\u00CC\u0082s(s)qs + CoC\u00CC\u0082o(s)qo)\u00E2\u0088\u0082T \u00E2\u0088\u0082x = D \u00E2\u0088\u0082 \u00E2\u0088\u0082x (D\u00CC\u0082(s) \u00E2\u0088\u0082T \u00E2\u0088\u0082x ) 99 3.3. 1D Temperature Model where, recalling from Table 2.2 we have that, Cs = 0.405 Co = 0.105 D = 7.24\u00C3\u0097 10\u00E2\u0088\u009210. 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, CsC\u00CC\u0082s = CoC\u00CC\u0082o = 1, which essentially amounts to assuming the volumetric heat capacities, \u00CF\u0081cp are equal. Now since we are analyzing the 1D model, the conservation of volume (2.4) once again amounts to, qs + qo = 1 so that we are able to simplify the temperature dependence to simply being the hyperbolic partial differential equation, \u00E2\u0088\u0082T \u00E2\u0088\u0082t + \u00E2\u0088\u0082T \u00E2\u0088\u0082x = 0, \u00E2\u0088\u0092\u00E2\u0088\u009E < x <\u00E2\u0088\u009E. 10There will be a boundary layer of width O(\u00E2\u0088\u009ADt) which means that as time progresses the temperature front will become smooth but this is on a long time scale and we neglect it. 100 3.3. 1D Temperature Model We will take for initial conditions to this problem, Riemann data such that there is a \u00E2\u0080\u009Chot\u00E2\u0080\u009D front coming in from the left, x < 0 with T = 1 and then everything on the right x > 0 is \u00E2\u0080\u009Ccold\u00E2\u0080\u009D with T = 0, T (x, 0) = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B31, x < 00, x > 0 . With this consideration then we can easily solve the temperature profile to get that, T (x, t) = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B31, x < t0, x > t . We now have solved for the temperature and thus can easily compute the viscosity. We have from (2.21d) that, \u00C2\u00B5o(T ) = (1\u00E2\u0088\u0092 \u00CE\u00B4) exp(1\u00E2\u0088\u0092 T ) + \u00CE\u00B4 however upon the knowledge that for our particular problem \u00CE\u00B4 \u001C 1 we will simply ignore \u00CE\u00B4 in the viscosity. Therefore using the temperature we can get that the oil viscosity is approximately, \u00C2\u00B5o(T ) = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B31, x < texp(1), x > t 101 3.3. 1D Temperature Model and of course we have that the water viscosity is \u00C2\u00B5s(T ) \u00E2\u0089\u00A1 1. Having defined this we can now write our velocities as, qs = \u00E2\u0088\u0092\u00CE\u00B1sdP dx qo = \u00E2\u0088\u0092\u00CE\u00B4 \u00CE\u00B1o \u00C2\u00B5o(T ) dP dx 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 ) = \u00CE\u00B1s(s) \u00CE\u009B(s, T ) where we define \u00CE\u009B as the modified total mobility, \u00CE\u009B(s, T ) = \u00CE\u00B1s(s) + \u00CE\u00B4 \u00CE\u00B1o(s) \u00C2\u00B5o(T ) . For our currently specified temperature we have that, \u00CE\u009B = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3\u00CE\u00B1s + \u00CE\u00B4\u00CE\u00B1o = \u00CE\u00BB, x < t\u00CE\u00B1s + \u00CE\u00B4 exp(\u00E2\u0088\u00921)\u00CE\u00B1o, x > t and therefore we have that, qs = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3 s3 s3+\u00CE\u00B4(1\u00E2\u0088\u0092s)3 = f1, x < t s3 s3+\u00CE\u00B4 exp(\u00E2\u0088\u00921)(1\u00E2\u0088\u0092s)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, \u00E2\u0088\u0082s1 \u00E2\u0088\u0082t + \u00E2\u0088\u0082f1 \u00E2\u0088\u0082x = 0, \u00E2\u0088\u0092\u00E2\u0088\u009E < x < t (3.43a) \u00E2\u0088\u0082s2 \u00E2\u0088\u0082t + \u00E2\u0088\u0082f2 \u00E2\u0088\u0082x = 0, t < x <\u00E2\u0088\u009E (3.43b) where, for initial conditions, we will impose the Riemann data as in the scenario which did not include temperature, s(x, 0) = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B31, x < 00, 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 bound- ary 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 grow- ing 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+ \u00E2\u0088\u0086t as shown in Figure 3.12. 103 3.3. 1D Temperature Model cold t t+ \u00E2\u0088\u0086t vn hot H 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\u00E2\u0088\u0086t 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+\u00E2\u0088\u0086t and thus we can write this saturation differential as, Sd = (s2 \u00E2\u0088\u0092 s1)Hvn\u00E2\u0088\u0086t. Now as the front moves there is also a net change in the flux of water over the time period t+ \u00E2\u0088\u0086t and thus we get a flux differential, Sf = (f2 \u00E2\u0088\u0092 f1)H\u00E2\u0088\u0086t. 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 \u00E2\u0088\u0092 f1 s2 \u00E2\u0088\u0092 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, \u00E2\u0088\u0082s1 \u00E2\u0088\u0082t + f \u00E2\u0080\u00B21(s1) \u00E2\u0088\u0082s1 \u00E2\u0088\u0082x = 0, \u00E2\u0088\u0092\u00E2\u0088\u009E < 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 \u00E2\u0088\u0092\u00E2\u0088\u009E < x < 0 which satisfies the requirements. After applying the method we get that along vertical charac- teristics x = \u00CE\u00B21, 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 = \u00CE\u00B22. 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, \u00E2\u0088\u0082si \u00E2\u0088\u0082t + f \u00E2\u0080\u00B2i \u00E2\u0088\u0082si \u00E2\u0088\u0082x = 0 without any initial data. We see that the solution is, si = C, a constant, along any characteristics, x = f \u00E2\u0080\u00B2i(C)t+ \u00CE\u00B2. 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 dt = 1 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., f \u00E2\u0080\u00B21(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, f \u00E2\u0080\u00B22(C) > 1. Each scenario produces a unique characteristic plot and hence a unique so- lution 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 = f \u00E2\u0080\u00B21(C)t+ \u00CE\u00B2 with f \u00E2\u0080\u00B21(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 = \u00CE\u00B21) to these new linear ones. For the rarefaction fan the characteristics will emanate from the origin and thus be a function of \u00CF\u0087 = 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 = f \u00E2\u0080\u00B21(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 = \u00CE\u00B22) 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 = f \u00E2\u0080\u00B22(s)t, s \u00E2\u0088\u0097 < s < s+ where, like before, s\u00E2\u0088\u0097, 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, f \u00E2\u0080\u00B22(s +) = 1. Like in the isothermal case, the shock that forms (x = \u00CE\u00BE) must be the fan characteristic associated with s\u00E2\u0088\u0097 and so we know that, d\u00CE\u00BE dt = f \u00E2\u0080\u00B22(s \u00E2\u0088\u0097) however the shock condition must also satisfy a Rankine-Hugoniot condition, d\u00CE\u00BE dt = f \u00E2\u0080\u00B22(s \u00E2\u0088\u0097) = f2(s \u00E2\u0088\u0097) s\u00E2\u0088\u0097 which determines s\u00E2\u0088\u0097. While there may be multiple solutions for s\u00E2\u0088\u0097 we will be able to determine it uniquely, discarding any solutions that do not satisfy 0 \u00E2\u0089\u00A4 s \u00E2\u0089\u00A4 1 and s\u00E2\u0088\u0097 < s+. Finally, we can determine the constant C via the Stefan condition (3.44) to get that, f1(C)\u00E2\u0088\u0092 f2(s+) C \u00E2\u0088\u0092 s+ = 1. 109 3.3. 1D Temperature Model While, like with s\u00E2\u0088\u0097, there may be multiple solutions for C, it can be deter- mined 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 f \u00E2\u0080\u00B21(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, s(x, t) = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3 1, x < 0 [f \u00E2\u0080\u00B21] \u00E2\u0088\u00921 (x t ) , 0 < x < f \u00E2\u0080\u00B21(C)t C, f \u00E2\u0080\u00B21(C)t < x < t [f \u00E2\u0080\u00B22] \u00E2\u0088\u00921 (x t ) , t < x < f \u00E2\u0080\u00B22(s \u00E2\u0088\u0097)t 0, x > f \u00E2\u0080\u00B22(s \u00E2\u0088\u0097)t . (3.45) Figure 3.13 shows the characteristic diagram for this scenario. 110 3.3. 1D Temperature Model t = 0 x = t x = 0s = 1 s = C x = \u00CE\u00BE s = 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 = \u00CE\u00BE 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 \u00E2\u0088\u0092 1) 111 3.3. 1D Temperature Model instead of our current, exp (1\u00E2\u0088\u0092 T ) which has the viscosity increasing. This just places a term exp(1) into the cold water flux, f2 instead of exp(\u00E2\u0088\u00921) so that, f1 = s3 s3 + \u00CE\u00B4(1\u00E2\u0088\u0092 s)3 f2 = s3 s3 + \u00CE\u00B4 exp(1)(1\u00E2\u0088\u0092 s)3 . 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 f \u00E2\u0080\u00B22(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, f \u00E2\u0080\u00B22(s \u00E2\u0088\u0097) = f2(s \u00E2\u0088\u0097) s\u00E2\u0088\u0097 to get that s\u00E2\u0088\u0097 = 0.194. Finally we use the Stefan condition to get that C = 0.348. Notice we indeed have that f \u00E2\u0080\u00B21(C) < 1, a requirement of this case and that C > s+ > s\u00E2\u0088\u0097 > 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 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 0 0.2 0.4 0.6 0.8 1 x s, T Temperature influence on water saturation at t=0.0625 Numerical Temperature Analytical 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 = f \u00E2\u0080\u00B22(C) + \u00CE\u00B2 with f \u00E2\u0080\u00B22(C) > 1. These characteristics which carry the value s = C will intersect the vertical characteristics with s = 0 forming a classic shock, \u00CE\u00BE where the shock speed is given by the Rankine-Hugoniot condition, d\u00CE\u00BE dt = f2(C) C . Here the shock saturation value C, unlike s\u00E2\u0088\u0097, is determined elsewhere. Now in the hot region, we will need a rarefaction fan to connect the vertical char- acteristics of s = 1 to the internal interface at x = t. These characteristics will be of the form, x = f \u00E2\u0080\u00B21(s)t, s \u00E2\u0088\u0092 < s < 1 where s\u00E2\u0088\u0092 is the value of s where the expansion fan intersects the internal interface and thus it satisfies, f \u00E2\u0080\u00B21(s \u00E2\u0088\u0092) = 1. We then use our Stefan condition (3.44) to determine C, 1 = f1(s \u00E2\u0088\u0092)\u00E2\u0088\u0092 f2(C) s\u00E2\u0088\u0092 \u00E2\u0088\u0092 C . 114 3.3. 1D Temperature Model Like before, invalid answers for C, as well as those that do not satisfy f \u00E2\u0080\u00B22(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) = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3 1, x < 0 [f \u00E2\u0080\u00B21] \u00E2\u0088\u00921 (x t ) , 0 < x < t C, t < x < \u00CE\u00BE 0, x > \u00CE\u00BE . (3.46) Figure 3.15 shows the characteristic diagram for this case. s = 0 x = 0s = 1 x = t t = 0 s = C x = \u00CE\u00BE 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 = \u00CE\u00BE 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 f \u00E2\u0080\u00B21(C) < 0 or f \u00E2\u0080\u00B2 2(C) > 0 and have ignored any other possibilities. Firstly we have ignored the unphysical case where f \u00E2\u0080\u00B21(C) > 0(< 0) but f \u00E2\u0080\u00B22(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 f \u00E2\u0080\u00B21(C) < 1, in which case, f \u00E2\u0080\u00B21(C) = \u00CE\u00B4 \u00C2\u00B501 \u00CE\u00B1o(C)\u00CE\u00B1 \u00E2\u0080\u00B2 s(C)\u00E2\u0088\u0092 \u00CE\u00B1\u00E2\u0080\u00B2o(C)\u00CE\u00B1s(C)( \u00CE\u00B1s(C) + \u00CE\u00B4 \u00C2\u00B501 \u00CE\u00B1o(C) )2 < 1 which can be rearranged as, \u00CE\u00B4 (\u00CE\u00B1o(C)\u00CE\u00B1 \u00E2\u0080\u00B2 s(C)\u00E2\u0088\u0092 \u00CE\u00B1\u00E2\u0080\u00B2o(C)\u00CE\u00B1s(C)) < \u00C2\u00B501 ( \u00CE\u00B1s(C) + \u00CE\u00B4 \u00C2\u00B501 \u00CE\u00B1o(C) )2 . (3.47) Then, we can write f \u00E2\u0080\u00B22(C) as, f \u00E2\u0080\u00B22(C) = \u00CE\u00B4 \u00C2\u00B502 \u00CE\u00B1o(C)\u00CE\u00B1 \u00E2\u0080\u00B2 s(C)\u00E2\u0088\u0092 \u00CE\u00B1\u00E2\u0080\u00B2o(C)\u00CE\u00B1s(C)( \u00CE\u00B1s(C) + \u00CE\u00B4 \u00C2\u00B502 \u00CE\u00B1o(C) )2 < \u00C2\u00B501 \u00C2\u00B502 ( \u00CE\u00B1s(C) + \u00CE\u00B4 \u00C2\u00B501 \u00CE\u00B1o(C) )2 ( \u00CE\u00B1s(C) + \u00CE\u00B4 \u00C2\u00B502 \u00CE\u00B1o(C) )2 using inequality (3.47). Now assume T1 > T2 so that \u00C2\u00B501 < \u00C2\u00B502 and then, f \u00E2\u0080\u00B22(C) < ( \u00CE\u00B1s(C) + \u00CE\u00B4 \u00C2\u00B501 \u00CE\u00B1o(C) )2 ( \u00CE\u00B1s(C) + \u00CE\u00B4 \u00C2\u00B502 \u00CE\u00B1o(C) )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, f \u00E2\u0080\u00B22(C) / \u00CE\u00B1s(C) \u00CE\u00B1s(C) = 1 and thus if f \u00E2\u0080\u00B21(C) < 1 then so too is f \u00E2\u0080\u00B2 2(C). As \u00CE\u00B4 \u001C 1 this approximation holds stronger even for C not near 1. We can extend these ideas to the other cases easily when \u00CE\u00B4 \u001C 1 to show that if f \u00E2\u0080\u00B2i < 1(> 1) then f \u00E2\u0080\u00B2j < 1(> 1) as well so that characteristics will either emerge into the hot or cold region exclusively. Admittedly the analysis doesn\u00E2\u0080\u0099t hold for \u00CE\u00B4 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 f \u00E2\u0080\u00B21 < 1 and f \u00E2\u0080\u00B2 2 > 1 should hold. Lastly, there is the case when f \u00E2\u0080\u00B21(C) = f \u00E2\u0080\u00B2 2(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\u00E2\u0080\u0099ve 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 f \u00E2\u0080\u00B21(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 f \u00E2\u0080\u00B21(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 f \u00E2\u0080\u00B21 = 1 to get that, s\u00E2\u0088\u0092 = 0.233. Now we use this in the Stefan condition to get that, C = 0.133 which does indeed satisfy f \u00E2\u0080\u00B22(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\u00CE\u00BE dt = f2(C) C = 6.30. 11This may make our claims above that if f \u00E2\u0080\u00B21(C) < 1 then C is relatively close to 1 seem absurd, however, since \u00CE\u00B4 \u001C 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. 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 0 0.2 0.4 0.6 0.8 1 x s, T Temperature influence on water saturation at t=0.0625 Numerical Temperature Analytical 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, \u00CE\u00B7 = x\u00E2\u0088\u0092 vt for some wave speed v then we have a steady state solution if, d d\u00CE\u00B7 (f \u00E2\u0088\u0092 vs) = 0 where it is assumed here that, f = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3f1, x < tf2, x > t . If, v = d\u00CE\u00BE dt then the coordinate frame is centered around the classical shock front and if we integrate, \u00E2\u0088\u00AB d d\u00CE\u00B7 (f \u00E2\u0088\u0092 vs) d\u00CE\u00B7 = f \u00E2\u0088\u0092 vs = D 120 3.3. 1D Temperature Model for D, a constant, then a steady state exists. We see that for D = 0, f2(C) C = v, 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 \u00E2\u0088\u0092)\u00E2\u0088\u0092 f2(C) s\u00E2\u0088\u0092 \u00E2\u0088\u0092 C = 1, 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 tem- perature is constant across this interface. In this scenario we are centered around the travelling wave coordinate, \u00CE\u00B7 = x\u00E2\u0088\u0092 vt 121 3.3. 1D Temperature Model with, v = d\u00CE\u00BE dt . Here we have a steady state solution that satisfies, d d\u00CE\u00B7 (f2 \u00E2\u0088\u0092 vs0) = 0 which has an associated pressure steady state, dP0 d\u00CE\u00B7 = \u00E2\u0088\u0092 1 \u00CE\u009B(s0) . We will do the standard perturbations for s, P (using steady states as above), and the interface \u00CE\u00B7\u00CC\u0082 such that, s \u00E2\u0088\u00BC s0 + \u00CE\u00B5s1(\u00CE\u00B7) exp (iny + \u00CF\u0083t) P \u00E2\u0088\u00BC P0 + \u00CE\u00B5P1(\u00CE\u00B7) exp (iny + \u00CF\u0083t) \u00CE\u00B7\u00CC\u0082 \u00E2\u0088\u00BC 0 + \u00CE\u00B5A exp (iny + \u00CF\u0083t) s1, ds1 d\u00CE\u00B7 , P1, dP1 d\u00CE\u00B7 \u00E2\u0086\u0092 0, |\u00CE\u00B7| \u00E2\u0086\u0092 \u00E2\u0088\u009E where \u00CE\u00B5 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, \u00E2\u0088\u0082s \u00E2\u0088\u0082t +\u00E2\u0088\u0087 \u00C2\u00B7 (qs \u00E2\u0088\u0092 vsx\u00CC\u0082) = 0 qs = \u00E2\u0088\u0092\u00CE\u00B1s\u00E2\u0088\u0087P qo = \u00E2\u0088\u0092\u00CE\u00B4\u00CE\u00B1o \u00C2\u00B5o \u00E2\u0088\u0087P \u00E2\u0088\u0087 \u00C2\u00B7 (qs + qo) = 0 where we have once again redefined the gradient operator in the travelling wave as, \u00E2\u0088\u0087 = \u00E2\u008C\u00A9 \u00E2\u0088\u0082 \u00E2\u0088\u0082\u00CE\u00B7 , \u00E2\u0088\u0082 \u00E2\u0088\u0082y \u00E2\u008C\u00AA . When we apply our perturbations to this system we get for our eigenvalue problem, d d\u00CE\u00B7 ( \u00CE\u00B1\u00CC\u0084s dP1 d\u00CE\u00B7 + \u00CE\u00B1\u00CC\u0084s \u00E2\u0080\u00B2dP0 d\u00CE\u00B7 s1 + vs1 ) \u00E2\u0088\u0092 n2\u00CE\u00B1\u00CC\u0084sP1 = \u00CF\u0083s1 (3.48a) d d\u00CE\u00B7 ( \u00CE\u009B\u00CC\u0084 dP1 d\u00CE\u00B7 + \u00CE\u009B\u00CC\u0084\u00E2\u0080\u00B2 dPo d\u00CE\u00B7 s1 ) \u00E2\u0088\u0092 n2\u00CE\u009B\u00CC\u0084P1 = 0 (3.48b) which is exactly (3.11) with \u00CE\u009B in place of \u00CE\u00BB. Once again we let over-bars represent functions evaluated at the base-state. Since the steady-state solu- tions are constant we can apply a similar methodology to the temperature 123 3.3. 1D Temperature Model independent case to get that, s1 =0, \u00CE\u00B7 6= 0 P1 = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3B+ exp (\u00E2\u0088\u0092n\u00CE\u00B7) , \u00CE\u00B7 > 0B\u00E2\u0088\u0092 exp (n\u00CE\u00B7) , \u00CE\u00B7 < 0 . We then conserve total flux and pressure across the interface to generate the solvability conditions, [ \u00CE\u009B\u00CC\u0084 dP1 d\u00CE\u00B7 ] = 0,[ \u00CE\u00B1\u00CC\u0084s dP1 d\u00CE\u00B7 ] = CA\u00CF\u0083, [P1] = \u00E2\u0088\u0092A [ dP0 d\u00CE\u00B7 ] , 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, \u00CF\u0083 n = v \u00CE\u009B\u00CC\u0084(C)\u00E2\u0088\u0092 \u00CE\u009B\u00CC\u0084(0) \u00CE\u009B\u00CC\u0084(C) + \u00CE\u009B\u00CC\u0084(0) . (3.49) When \u00CE\u00B4 \u001C 1 then based on the expression for \u00CE\u009B\u00CC\u0084, the leading term in the numerator looks like C3 > 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. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 2 4 6 8 10 12 14 16 18 20 n \u00EF\u0081\u00B3 Linear Stability Eigenvalue-No Capillary Pressure Without Temperature With Temperature 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, \u00CE\u00B7 = x\u00E2\u0088\u0092 vt has v = 1. Our steady state satisfies, d d\u00CE\u00B7 (f \u00E2\u0088\u0092 s0) = 0 with, f = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3f1, \u00CE\u00B7 < 0f2, \u00CE\u00B7 > 0 like before. there is also an associated pressure steady state, dP0 d\u00CE\u00B7 = \u00E2\u0088\u0092 1 \u00CE\u009B(s0, T0) but now there is also an associated temperature steady state, T0 = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B31, \u00CE\u00B7 < 00, \u00CE\u00B7 > 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 \u00CE\u00B7 = 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, \u00CE\u00B7\u00CC\u0082s and one describing the temperature front, \u00CE\u00B7\u00CC\u0082T . With this in mind we will expand as follows12, s \u00E2\u0088\u00BC s0 + \u00CE\u00B5s1(\u00CE\u00B7) exp (iny + \u00CF\u0083t) , P \u00E2\u0088\u00BC P0 + \u00CE\u00B5P1(\u00CE\u00B7) exp (iny + \u00CF\u0083t) , T \u00E2\u0088\u00BC T0 + \u00CE\u00B5T1(\u00CE\u00B7) exp (iny + \u00CF\u0083t) , \u00CE\u00B7\u00CC\u0082s \u00E2\u0088\u00BC 0 + \u00CE\u00B5As exp (iny + \u00CF\u0083t) , \u00CE\u00B7\u00CC\u0082T \u00E2\u0088\u00BC 0 + \u00CE\u00B5AT exp (iny + \u00CF\u0083t) , s1, ds1 d\u00CE\u00B7 , P1, dP1 d\u00CE\u00B7 , T1, dT1 d\u00CE\u00B7 \u00E2\u0086\u0092 0, |\u00CE\u00B7| \u00E2\u0086\u0092 \u00E2\u0088\u009E where once again \u00CE\u00B5 is introduced to mange the order of terms. We apply this to the two dimensional variant of the one dimensional model, \u00E2\u0088\u0082s \u00E2\u0088\u0082t +\u00E2\u0088\u0087 \u00C2\u00B7 (qs \u00E2\u0088\u0092 sx\u00CC\u0082) = 0 qs = \u00E2\u0088\u0092\u00CE\u00B1s\u00E2\u0088\u0087P qo = \u00E2\u0088\u0092\u00CE\u00B4\u00CE\u00B1o \u00C2\u00B5o \u00E2\u0088\u0087P \u00E2\u0088\u0087 \u00C2\u00B7 (qs + qo) = 0 \u00E2\u0088\u0082T \u00E2\u0088\u0082t + (qs + qo \u00E2\u0088\u0092 x\u00CC\u0082) \u00C2\u00B7 \u00E2\u0088\u0087T = 0 12Like 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, d d\u00CE\u00B7 ( \u00CE\u00B1\u00CC\u0084s dP1 d\u00CE\u00B7 + \u00CE\u00B1\u00CC\u0084s \u00E2\u0080\u00B2dP0 d\u00CE\u00B7 s1 + s1 ) \u00E2\u0088\u0092 n2\u00CE\u00B1\u00CC\u0084sP1 = \u00CF\u0083s1 (3.50a) d d\u00CE\u00B7 \u00EF\u00A3\u00AB\u00EF\u00A3\u00AD\u00CE\u009B\u00CC\u0084dP1 d\u00CE\u00B7 + \u00CE\u009B\u00CC\u0084\u00E2\u0080\u00B2 dPo d\u00CE\u00B7 s1 \u00E2\u0088\u0092 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o d\u00C2\u00B5\u00CC\u0084 dT0 \u00C2\u00B5\u00CC\u00842 dP0 d\u00CE\u00B7 T1 \u00EF\u00A3\u00B6\u00EF\u00A3\u00B8\u00E2\u0088\u0092 n2\u00CE\u009B\u00CC\u0084P1 = 0 (3.50b) \u00CE\u009B\u00CC\u0084 dT0 d\u00CE\u00B7 dP1 d\u00CE\u00B7 + ( \u00CE\u009B\u00CC\u0084 dP0 d\u00CE\u00B7 + 1 ) dT1 d\u00CE\u00B7 + \u00EF\u00A3\u00AB\u00EF\u00A3\u00AD\u00CE\u009B\u00CC\u0084\u00E2\u0080\u00B2 \u00E2\u0088\u0092 \u00CE\u00B4\u00CE\u00B1\u00CC\u0084o d\u00C2\u00B5\u00CC\u0084dT0 \u00C2\u00B5\u00CC\u00842 \u00EF\u00A3\u00B6\u00EF\u00A3\u00B8 dP0 d\u00CE\u00B7 dT0 d\u00CE\u00B7 = \u00CF\u0083T1. (3.50c) 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 \u00CE\u00B7 = 0 where steady-states are well defined. In this case then we can use the fact that, dT0 d\u00CE\u00B7 = 0 dP0 d\u00CE\u00B7 = \u00E2\u0088\u0092 1 \u00CE\u009B\u00CC\u0084 and substitute it into (3.50c) to get, \u00CF\u0083T1 = 0 128 3.3. 1D Temperature Model which for non-zero eigenvalues means that, T1 \u00E2\u0089\u00A1 0, \u00CE\u00B7 6= 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, \u00CE\u00B7 6= 0 P1 = \u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3B+ exp (\u00E2\u0088\u0092n\u00CE\u00B7) , \u00CE\u00B7 > 0B\u00E2\u0088\u0092 exp (n\u00CE\u00B7) , \u00CE\u00B7 < 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 \u00CE\u00B4(\u00CE\u00B7). 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 con- tinuous 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\u00E2\u0080\u0099d expect 129 3.3. 1D Temperature Model from the previous problems, [ \u00CE\u009B\u00CC\u0084 dP1 d\u00CE\u00B7 ] = 0[ \u00CE\u00B1\u00CC\u0084s dP1 d\u00CE\u00B7 ] = (s\u00E2\u0088\u0092 \u00E2\u0088\u0092 C)As\u00CF\u0083 [P1] = \u00E2\u0088\u0092As [ dP0 d\u00CE\u00B7 ] , taking into account the correct upstream and downstream saturation values. Now for the temperature front, pressure at \u00CE\u00B7 = \u00CE\u00B7\u00CC\u0082kT where k = + or k = \u00E2\u0088\u0092 depending on the sided limit of concern can be expanded as, P (\u00CE\u00B7\u00CC\u0082kT ) \u00E2\u0088\u00BC P0(0k) + \u00CE\u00B5 ( AT dP0 d\u00CE\u00B7 \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 \u00CE\u00B7=0k + P1(0 k) ) exp (iny + \u00CF\u0083t) 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] = \u00E2\u0088\u0092AT [ dP0 d\u00CE\u00B7 ] . Now the thermal flux in the travelling wave coordinate frame is given by, (qs + qo \u00E2\u0088\u0092 x\u00CC\u0082)T. 130 3.3. 1D Temperature Model Physically the jump in normal thermal flux must be given by, [x\u00CC\u0082 \u00C2\u00B7 (qs + qo)\u00E2\u0088\u0092 T ] = ( total flux across \u00CE\u00B7 = \u00CE\u00B7\u00CC\u0082+T )\u00E2\u0088\u0092 (total flux across \u00CE\u00B7 = \u00CE\u00B7\u00CC\u0082\u00E2\u0088\u0092T ) \u00E2\u0088\u00BC (total flux across \u00CE\u00B7 = 0+)\u00E2\u0088\u0092 (total flux across \u00CE\u00B7 = 0\u00E2\u0088\u0092) = (s(0+) + o(0+))T (0+) d\u00CE\u00B7\u00CC\u0082T dt \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 \u00CE\u00B7=0+ \u00E2\u0088\u0092 (s(0\u00E2\u0088\u0092) + o(0\u00E2\u0088\u0092))T (0\u00E2\u0088\u0092) d\u00CE\u00B7\u00CC\u0082T dt \u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3\u00E2\u0088\u00A3 \u00CE\u00B7=0\u00E2\u0088\u0092 with, d\u00CE\u00B7\u00CC\u0082T dt = \u00CE\u00B5\u00CF\u0083AT exp (iny + \u00CF\u0083t) being the speed of the temperature interface. However, we also have that, [x\u00CC\u0082 \u00C2\u00B7 (qs + qo)\u00E2\u0088\u0092 T ] = [ \u00E2\u0088\u0092\u00CE\u009B\u00CC\u0084dP0 d\u00CE\u00B7 T0 ] \u00E2\u0088\u0092 [T0] + \u00CE\u00B5 ( \u00E2\u0088\u0092\u00CE\u009B\u00CC\u0084dP1 d\u00CE\u00B7 T0 ) exp (iny + \u00CF\u0083t) where we have explicitly used that, s1 = T1 = dT0 d\u00CE\u00B7 = 0. Combining the two expressions we get that, [ \u00E2\u0088\u0092\u00CE\u009B\u00CC\u0084dP0 d\u00CE\u00B7 T0 ] = [T0] which is indeed satisfied and more importantly that, [ \u00CE\u009B\u00CC\u0084 dP1 d\u00CE\u00B7 T0 ] = AT\u00CF\u0083. 131 3.3. 1D Temperature Model We then have a complete set of solvability conditions in order to have non- trivial solutions, [ \u00CE\u009B\u00CC\u0084 dP1 d\u00CE\u00B7 ] = 0, (3.51a)[ \u00CE\u00B1\u00CC\u0084s dP1 d\u00CE\u00B7 ] = (s\u00E2\u0088\u0092 \u00E2\u0088\u0092 C)As\u00CF\u0083, (3.51b)[ \u00CE\u009B\u00CC\u0084 dP1 d\u00CE\u00B7 T0 ] = AT\u00CF\u0083, (3.51c) [P1] = \u00E2\u0088\u0092(As + AT ) [ dP0 d\u00CE\u00B7 ] , (3.51d) which are the jump conditions for this thermal problem. We can write (3.51) in a matrix problem as, \u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 (s\u00E2\u0088\u0092 \u00E2\u0088\u0092 C)\u00CF\u0083 n \u00CE\u00B1\u00CC\u0084s(s \u00E2\u0088\u0092) \u00CE\u00B1\u00CC\u0084s(C) 0( 1 \u00CE\u009B\u00CC\u0084(C) \u00E2\u0088\u0092 1 \u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092) ) 1 \u00E2\u0088\u00921 ( 1 \u00CE\u009B\u00CC\u0084(C) \u00E2\u0088\u0092 1 \u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092) ) 0 \u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092) \u00CE\u009B\u00CC\u0084(C) 0 0 \u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092) 0 \u00CF\u0083 n \u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB \u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 As C\u00E2\u0088\u0092 C+ AT \u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB = \u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 0 0 0 0 \u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB 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 \u00CF\u0083/n, ( 1 \u00CE\u009B\u00CC\u0084(C) \u00E2\u0088\u0092 1 \u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092) ) (s\u00E2\u0088\u0092 \u00E2\u0088\u0092 C)\u00CF\u0083 n ( \u00CE\u009B\u00CC\u0084(C)\u00E2\u0088\u0092 \u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092)) \u00E2\u0088\u0092 \u00CF\u0083 n ( (s\u00E2\u0088\u0092 \u00E2\u0088\u0092 C)\u00CF\u0083 n ( \u00CE\u009B\u00CC\u0084(C) + \u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092) )\u00E2\u0088\u0092 ( 1 \u00CE\u009B\u00CC\u0084(C) \u00E2\u0088\u0092 1 \u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092) )) \u00E2\u0088\u0092 \u00CF\u0083 n (( \u00CE\u00B1\u00CC\u0084s(s \u00E2\u0088\u0092)\u00CE\u009B\u00CC\u0084(C)\u00E2\u0088\u0092 \u00CE\u00B1\u00CC\u0084s(C)\u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092) )) = 0 132 3.3. 1D Temperature Model which we can solve to get the solutions, \u00CF\u0083 = 0 and \u00CF\u0083 n = 2 \u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092)\u00E2\u0088\u0092 \u00CE\u009B\u00CC\u0084(C) \u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092) + \u00CE\u009B\u00CC\u0084(C) . (3.52) When \u00CE\u00B4 \u001C 1 then the leading order term in the numerator is (s\u00E2\u0088\u0092)3 \u00E2\u0088\u0092 C3 which is positive since s\u00E2\u0088\u0092 > 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, \u00CF\u0083 n = 2v \u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092)\u00E2\u0088\u0092 \u00CE\u009B\u00CC\u0084(C) \u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092) + \u00CE\u009B\u00CC\u0084(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 \u00CF\u0083 = 0 eigenvalue then we get from the jump conditions (3.51) that there will only be a solution if, As = \u00E2\u0088\u0092AT 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 = \u00E2\u0088\u0092AT then \u00CF\u0083 = 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 non- degenerate front. Indeed, if we substitute the expression for \u00CF\u0083 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. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 1 2 3 4 5 6 7 n \u00EF\u0081\u00B3 Linear Eigenvalue Growth Rate for Various Thermal Gradients \u00EF\u0081\u0084T=1 \u00EF\u0081\u0084T=1/3 \u00EF\u0081\u0084T=1/30 \u00EF\u0081\u0084T=1/300 Figure 3.19: Fingering growth rate for various values of non-dimensional temperature gradient, \u00E2\u0088\u0086T . The solid line is \u00E2\u0088\u0086T = 1, the dashed line is \u00E2\u0088\u0086T = 1/3, the dot-dashed line is \u00E2\u0088\u0086T = 1/30 and the dotted line is \u00E2\u0088\u0086T = 1/300. The second conclusion that [37] draws is that despite the lack of qualita- tive difference in fingers, they notice that the single fingering pattern (the classical shock) travels at a slower transverse speed in the presence of a ther- mal 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 + \u00CE\u00B4 exp (\u00E2\u0088\u00921) (1\u00E2\u0088\u0092 s)3 . If we were to compute an isothermal calculation using the same flux function, i.e. finding the shock saturation value s\u00E2\u0088\u0097 that satisfies, f \u00E2\u0080\u00B22(s \u00E2\u0088\u0097) = f2(s \u00E2\u0088\u0097) s\u00E2\u0088\u0097 then we would get that the isothermal velocity vo is, vo = f2(s \u00E2\u0088\u0097) = 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 temper- ature gradient we have in our problem (\u00E2\u0088\u0086T = 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 veloc- ity 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, (2vT \u00E2\u0088\u0092 vc)\u00CE\u009B\u00CC\u0084(C)(\u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092)\u00E2\u0088\u0092 \u00CE\u009B\u00CC\u0084(0)) + (2vT + vc)(\u00CE\u009B\u00CC\u0084(s\u00E2\u0088\u0092)\u00CE\u009B\u00CC\u0084(0)\u00E2\u0088\u0092 \u00CE\u009B\u00CC\u0084(C)2) < 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\u00E2\u0080\u0099s 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 en- ergy 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. How- ever, 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 ad- mits a shock emerging from a rarefaction wave which connects an upstream shock saturation value, s\u00E2\u0088\u0097 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 sta- bility 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 dom- inated. 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. Asymp- totically we showed that to first order, the linear eigenvalue still dominated as expected so that in the limit of vanishing capillary effect, the Buckley- Leverett 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 \u000F-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, \u00CF\u0081cp 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 satura- tion 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 intro- duction 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 hy- pothesis 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 spe- cific application to the parameters associated with remediation. We have demonstrated that the inclusion of surface tension between the oil and wa- ter, even if nearly negligible, has an immense effect on the stability properties of the associated linear eigenvalue. This statement holds true for any prob- lem 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 \u00E2\u0088\u00BC u1/3, u\u00E2\u0086\u0092 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 nu- meric 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 sta- bility 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 re- mediation 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 bound- ary 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 fin- ger 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\u00E2\u0080\u0099s long- time 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 un- derlying 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 vis- cosity of benzene, toluene, and m-xylene at pressure up to 80 mpa. International Journal of Thermophysics, 12:449\u00E2\u0080\u0093457, 1991. [3] J. Baba and P. D. Komar. Measurements and analysis of setting veloc- ities of natural quartz sand grains. Journal of Sedimentary Research, 51(2):631\u00E2\u0080\u0093640, 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\u00E2\u0080\u0093396, 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\u00E2\u0080\u0093268, 2003. [6] S.E. Buckley and M.C. Leverett. Mechanism of fluid displacement in sands. Transactions of the AIME, 146:107\u00E2\u0080\u0093116, 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\u00E2\u0080\u00931308. [9] A.H. Craven and Peletier L.A. Similarity solutions for degenerate quasi- linear parabolic equations. Journal of Mathematical Analysis and Ap- plications, 38:73\u00E2\u0080\u009381, 1972. [10] M.B. Dusseault. Comparing Venezuelan and Canadian Heavy Oil and Tar Sands. Technical report, Canadian International Petroleum Con- ference, Calgary, Canada. [11] G. Elert. Viscosity. http://http://physics.info/viscosity/. Ac- cessed 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 Ver- sion 2011). CRC Press/Taylor and Francis. [16] M.G. Hennessy. Liquid snowflake formation in superheated ice. Master\u00E2\u0080\u0099s thesis, University of Oxford, 2010. [17] R.N. Hills and P.H. Roberts. A note on the kinetic conditions at a su- percooled interface. Int. Comm. Heat Mass Transfer, 20:407\u00E2\u0080\u0093416, 1993. [18] E.J. Hinch. Perturbation Methods. Cambridge Texts in Applied Math- ematics. 1991. [19] Kristi E Holloway and John R de Bruyn. Viscous fingering with a single fluid. Canadian Journal of Physics, 83(5):551\u00E2\u0080\u0093564, 2005. [20] G.M. Homsy. Viscous fingering in porous media. Annual Review Fluid Mechanics, 19:271\u00E2\u0080\u0093311, 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\u00E2\u0080\u0093513, September 2001. [23] L. Korson, W. Drost-Hansen, and F.J. Millero. Viscosity of water at various temperatures. J. Phys. Chem., 73(1):34\u00E2\u0080\u009339, January 1969. 150 Bibliography [24] A.A. Lacey, J.R. Ockendon, and A.B. Tayler. \u00E2\u0080\u009Cwaiting-time\u00E2\u0080\u009D solutions of a nonlinear differential equation. SIAM Journal on Applied Mathe- matics, 42(6), December. [25] M.C. Leverett. Capillary behaviour in porous solids. Transactions of the AIME, 142:159\u00E2\u0080\u0093172, 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\u00E2\u0080\u0093365, 1984. [27] C. Ma and J. Kingscott. Recent Developments for In Situ Treatment of Metal Contaminated Soils. Technical Report 68-W5-0055, U.S. En- vironment 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 Thermody- namics, 20(4):405 \u00E2\u0080\u0093 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\u00E2\u0080\u0093464, 1983. [30] I. Nazad, R.G. Carbonell, and S. Whitaker. Heat condition in multi- phase systems i: Theory and experiments for two-phase systems. Chem. Engng. Sci., 40:843\u00E2\u0080\u0093855, 1985. 151 Bibliography [31] J. Ockendon, S.D. Howison, and A.B. Movchan. Applied Partial Differ- ential 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 con- ductivity of water. Journal of Physical and Chemical Reference Data, 24(3):1377\u00E2\u0080\u00931381, 1995. [34] A. Riaz and H.A. Tchelepi. Linear stability analysis of immiscible two- phase flow in porous media with capillary dispersion and density varia- tion. Phys. Fluids, 16(12):4727\u00E2\u0080\u00934737, 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 \u00E2\u0080\u0093 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\u00E2\u0080\u00931350, 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\u00E2\u0080\u0093872, 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 mix- tures 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 dis- placement in porous media. SIAM J. Appl. Math., 49(3):730\u00E2\u0080\u0093748, June 1989. [47] Y.C. Yortsos and A.B. Huang. Linear-stability analysis of immisci- ble displacement:part 1-simple basic flows. SPE Reservoir Engineering, 1:378\u00E2\u0080\u0093390, 1986. 153"@en . "Thesis/Dissertation"@en . "2011-11"@en . "10.14288/1.0080593"@en . "eng"@en . "Mathematics"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "Attribution-NonCommercial-NoDerivs 3.0 Unported"@en . "http://creativecommons.org/licenses/by-nc-nd/3.0/"@en . "Graduate"@en . "Thermo-viscous fingering in porous media and in-situ soil remediation"@en . "Text"@en . "http://hdl.handle.net/2429/36926"@en .