Fluid Displacement and Solid Formation in a Porous Medium Using

Article pubs.acs.org/EF

Fluid Displacement and Solid Formation in a Porous Medium Using Invasion Percolation in a Gradient with Pore Blocking Ioannis N. Tsimpanogiannis*,†,§ and Peter C. Lichtner‡ †

Environmental Research Laboratory, National Center for Scientific Research “Demokritos”, Aghia Paraskevi, Athens, Greece 15310 Earth and Environmental Sciences Division, LosAlamos National Laboratory, Los Alamos, New Mexico 87545

‡

ABSTRACT: We utilize a modified form of the theory of invasion percolation in a gradient (IPG) in pore networks to simulate numerically immiscible displacement and reaction (solid formation and precipitation) in porous media. We examine the twodimensional invasion patterns that result, as well as the change in prosity due to the solid precipitation and the subsequent pore blocking and how they depend on various parameters including the solid formation kinetics and the strength of the percolation gradient. This work clearly demonstrates that for two-dimensional systems it is not the overall change of porosity in the system that controls the ability of the porous medium to accommodate flow, but its transverse average. This is because localized pore blocking occurring in narrow zones perpendicular to the invasion front can completely eliminate flow.

■

INTRODUCTION Fluid displacement within confined geometries, such as it occurs in porous media, is encountered in a large number of processes that are of interest to the engineering and natural sciences. Typical examples include water movement in the subsurface,1 oil/gas production,2 and soil remediation.3 In addition to nonreactive displacement, of interest are also cases where chemical reactions occur within the porous medium (e.g., flow and reaction inside catalysts in chemical reactors4,5 and reactive subsurface flows6). Chemical reactions within the porous medium can have a significant effect on the transport properties of the medium since they can enhance (during solid dissolution7,8) or reduce (during solid precipitation9) the flow capabilities of the medium by changing its porosity and permeability. To illustrate the importance of changing the transport properties of a porous medium through chemical reactions, consider the following two examples: (i) In an effort to mitigate the release of CO2 into the atmosphere (CO2 is considered a greenhouse gas10) a sequestration scheme in geologic media (e.g., saline aquifer, depleted oil reservoir, unminable coal seam, etc.) is proposed.11 During the injection process in the neighborhood of the injecting well, CO2 can dissolve the solid, which results in the formation of channels with increased permeability (“wormholes”), increasing thus the injectivity of the CO2 into the geologic formation. Within the same geologic formation, further downstream, CO2 can react with other minerals that are present and produce reactants that precipitate, causing a reduction in the porosity of the system. “Wormhole” formation generally occurs predominantly in limestone because of its high reactivity and high permeability. An excellent discussion on “wormhole” formation is presented by Golfier et al.12 (ii) Methane can be produced thermogenically from buried organic material deep under the ocean floor.13 The produced gas phase, due to buoyancy forces, starts © 2012 American Chemical Society

migrating upward, displacing the aqueous phase present inside the sediments. At certain conditions of pressure and temperature [i.e., inside the hydrate stability zone (HSZ)], stable methane gas hydrate can form, resulting in pore blockage. Hydrates are solid, nonstoichiometric, crystalline structures that are formed by hydrogenbonded water molecules arranged to form cavities (cages). The whole hydrate structure is stabilized by physical interactions of van der Waals-type between the water molecules and the methane gas14,15 that is encaged in the cavities. Hydrate formation within pore spaces eventually acts as a barrier to further upward methane migration and the formation of a methane-gas pool beneath the hydrate stability zone. In the absence of the hydrate zone, upward migrating methane could reach the ocean floor, and possibly the atmosphere, if the methane bubbles do not dissolve completely in the ocean water. This scenario could occur if the ocean is not deep enough for the bubbles to dissolve completely during the time they ascend to the ocean surface. The two aforementioned examples clearly demonstrate that small scale phenomena, which occur locally at the pore space (i.e., dissolution/precipitation), can affect the global percolation behavior of the porous medium. In general, large-scale calculations at the continuum level (in that approach phenomena occurring at the small scale are ignored and their effect is lumped into spatially averaged parameters) are essential for performing system-scale simulations. On the other hand, calculations focusing at the pore scale are essential for clarifying how small-scale effects can affect large-scale processes. Pore-network studies have contributed significantly in enhancing the understanding of phenomena occurring within porous media. Pore networks have primarily been used to study two-phase or multiphase displacement processes in porous Received: February 1, 2012 Revised: May 2, 2012 Published: May 8, 2012 3935

dx.doi.org/10.1021/ef300192x | Energy Fuels 2012, 26, 3935−3950

Energy & Fuels

Article

Figure 1. (a) Self-affine front during external drainage indicting an IPSG process. (b) Single finger in external drainage indicating an IPDG process. Invading phase is shown in black and defending phase in white.

media after their introduction in the pioneering work of Fatt.16 Excellent and detailed reviews have been provided by Feder,17 Sahimi,18,19 Celia et al.,20 Berkowitz and Ewing,21 and by Blunt and co-workers.22−24 Primarily, pore networks have been used to study relative permeabilities, electrical resistivities, and capillary pressures, the effects of gravity and viscous forces, and the effects of spatial correlation and wettability.17−25 Pore networks were also used to study the motion of solid particles through porous media.26,27 In recent years, pore networks have also been used to study phenomena where mass/heat transfer are important. In particular, pore networks were applied to study “liquid-to-gas” phase change processes including solution gas drive,28,29 drying,30−35 boiling,36 critical gas saturation during primary oil production37−39 and the “gas-to-liquid” phase change, as occurs in gas condensates.40 Pore networks were also used to study the “solid-to-liquid/gas” phase change that occurs during hydrate dissociation.41−43 Chemical reactions inside pore networks that result in solid-matrix dissolution were studied by Fogler and co-workers.44−47 Li et al.48−50 used pore networks to study CO2 sequestration in the subsurface. The problem of interest in this study can be considered as a combination of two-phase displacement with a phase change (“liquid-to-solid”) process, which can benefit from a porenetwork study. The purpose of this work is to examine the effect of solid formation kinetics on the nonwetting phase invasion patterns observed, and the change of porosity of the system as a result of pore blocking due to solid formation. The paper is organized as follows: First, we give a brief presentation of the theoretical background needed for a better understanding of the problem. Next, we introduce the description of the physical problem and the details of the pore-network simulator. We present the results of several simulations and discuss their implications on flow through a pore network. Finally, the limitations of the model are discussed, along with suggestions for possible improvements and their implication on the results.

porous medium is replaced by a lattice of sites and bonds connecting the sites. Pore bodies are located at the lattice sites and pore throats at the lattice bonds. Both pores and throats can be distributed following an appropriate size distribution and can be either correlated or not. At each time step only one pore throat is invaded. Namely, the one with the least capillary resistance (i.e., invasion threshold), among all throats at the perimeter of the wetting/ nonwetting interface that defines the invading front. In the absence of gravity or viscous forces (i.e., when there is no spatial or time dependence of the invasion threshold), this process is known to produce self-similar fractal patterns that resemble the ordinary percolation cluster.17,53 On the other hand, if gravity and/or viscous forces are important, the invasion threshold will have a gradient in space, which is known to produce self-affine fronts with a characteristic finite width.54,55 For the simpler case of gravity forces, the competition between gravity and capillary forces is given by the dimensionless Bond number (B): B=

g (ρdef − ρinv )k γ

(1)

where g is the component of gravity in the displacement direction and considered positive in the downward direction, ρ is the fluid density, and subscripts def and inv denote defending and invading phases, respectively; k is the permeability of the porous medium, and γ is the interfacial tension between the fluids. Two cases are possible: (i) If B > 0, as in the case of a lighter fluid displacing downward a heavier fluid, then the two fluids are separated by a self-affine front (see Figure 1a) that has a finite characteristic width σi, which scales with the Bond number as follows54,55 (Similar invasion patterns are obtained for the case of a heavier fluid displacing upward a lighter fluid.)

σi ∼ B−ν /(ν + 1)

■

(2)

For two-dimensional (2-D) geometries, the subscript i = f indicates the front. However, for three-dimensional (3-D) geometries, the subscript i = ft indicates the front tail. Rosso et al.56 and Gouyet et al.57 give a detailed discussion for the difference of σ f and σ ft in the two geometries. In eq 2, ν is the correlation length exponent for percolation (ν = 4/3 for 2-D invasion percolation with trapping, and ν ≃ 0.88 for 3-D

THEORETICAL BACKGROUND Invasion percolation51 (IP) is a stochastic approach that has been used extensively to model the slow immiscible displacement, controlled by capillary forces, of two phases within a porous medium.52 The case of a nonwetting phase displacing a wetting phase is also known as drainage. In this approach the 3936


Energy & Fuels

Article

invasion percolation with or without trapping53). If we introduce the values for ν into eq 2, the exponent becomes −0.57 for 2-D and −0.47 for 3-D. Such a process is also known as an invasion percolation in a stabilizing gradient (IPSG). Upstream of the front the pattern is compact, while around the front the pattern is a fractal with characteristics similar to the percolation cluster. In addition, one can define a mean front position, Xft, as

electrolytes, dissolved gases, etc. However, their presence and their effects, as well as the possible chemical reactions that lead to solid dissolution are ignored in the current study. It is known that, during immiscible displacement of a wetting fluid (e.g., water) by a nonwetting fluid (e.g., oil, supercritical carbon dioxide) inside a circular capillary, all wetting fluid will be displaced with the exception of a thin film covering the solid surface of the capillary. The resulting film thickness depends on the velocity of the advancing fluid (e.g., see the pioneering work of Bretherton65 for the case of displacement in a cylindrical tube). Thus, during the displacement of the wetting phase by an invading nonwetting phase, the wetting phase will be wetting the walls of the porous medium while the displacing nonwetting phase will occupy the center of the pores/throats. To simulate immiscible displacement in a porous medium, a pore-network representation in the form of a regular periodic lattice of size N × M is assumed, a special case of which is N = M. The lattice sites (pores) provide the storage capacity of the network, while the lattice bonds (throats) connecting neighboring sites, provide the resistance to flow. A typical schematic of such a network is depicted in Figure 2. Both pores

∞

X ft =

∫0 zpf (z) dz ∞

∫0 pf (z) dz

(3)

where pf(z) is the probability of finding one site of the front at z. The mean front position and the width of the front are related through ∞

2

σf =

∫0 (X ft − z)2 pf (z) dz ∞

∫0 pf (z) dz

(4)

The case of IPSG was examined experimentally by Clement et al.58,59 and Hulin et al.60 who injected Wood’s metal into nonconsolidated crushed glass to obtain measurements of the structure of the invasion fronts. (ii) If B < 0, as in the case of a heavier fluid displacing downward a lighter fluid, then the displacement takes the form of capillary fingering, where the heavy fluid invades in the form of fingers (see Figure 1b). Similar invasion patterns are obtained for the case of a lighter fluid displacing upward a heavier fluid. Each of these fingers is made up of a collection of (more or less) isotropic fractal blobs that retain the structure of an invasion percolation cluster at short length scales.61 The mean width of the fingers, ξ, is described by an equation that is similar to eq 2. The characteristic width ξ scales with the gravity Bond number as follows: ξ ∼ |B|−ν /(ν + 1)

(5)

Such a process is also known as an invasion percolation in a destabilizing gradient (IPDG). This process was examined in detail by Meakin et al.61 and Frette et al.62 The more complicated case of viscous forces will not be discussed in detail here. An extensive discussion can be found in the work of Yortsos and co-workers.63,64 However, we note that at small scales where capillary forces are important as well, the obtained patterns are similar to the case of gravity forces. In particular, depending on the value of the ratio M = M = μw/μnw, where μ is the viscosity and subscripts w and nw denote wetting and nonwetting phases, patterns similar to IPSG can be obtained for M < 1, and patterns similar to IPDG can be obtained for M > 1.63,64

Figure 2. Schematic of a square (N × M) site-bond lattice.

and throats can be randomly distributed, and in this study, they are also uncorrelated. Usually, pores are taken to be spherical while throats are cylindrical. Application of the methodology to a more realistic geometry of the porous medium such as using square,66−68 triangular69 or biconical70,71 throats is straightforward; however, it is more computationally demanding. The existence of corners within the pores/throats is usually associated with the presence of “thick” or “corner” films33,72,73 that can result during the immiscible displacement of a wetting phase by a nonwetting phase. Such films can facilitate viscous flows that are driven by the so-called “capillary pumping” mechanism74 and are responsible for producing interesting results, as is clearly shown by Chauvet et al.75 for the case of evaporation within a single square capillary tube. If pore/throat geometry that includes corners is considered, then liquid flow along the corners can be taken into account as well. That will result in minimizing the residual wetting saturation since the wetting phase can not be completely trapped (i.e., surrounded by nonwetting phase) as happens in the case of spherical pores or cylindrical throats, but the wetting phase can “leak out” through the corners.

■

DESCRIPTION OF THE PHYSICAL PROBLEM Consider a porous medium that is initially fully saturated with a wetting (w) defending phase (e.g., water or brine). The wetting phase is subsequently displaced immiscibly by a nonwetting (nw) invading phase. Typical examples of such a system include the injection of air into a contaminated soil during a remediation process (e.g., air sparging), the upward migration of methane displacing connate water inside the suboceanic sediments, and the injection of supercritical carbon dioxide in a saline aquifer during the process of geologic sequestration. Note that during a real geologic sequestration case scenario, the wetting phase could contain a variety of species including, 3937


Energy & Fuels

Article

depth of 3 km (i.e., pressure ≈ 300 bar) is a characteristic depth value were the densities of carbon dioxide and brine become approximately equal.77 The density of carbon dioxide is lower than the density of brine for depths less than 3 km, and becomes larger for depths over 3 km. Therefore, there is a narrow zone around the depth of 3 km were the Bond number can be considered B ≈ 0 for all practical purposes. The simulation algorithm consists of two distinct parts: the invasion part (displacement) and the reaction part (solid formation and precipitation). So far, the discussion given pertains to the invasion part of the algorithm. The description of the reaction part follows. We assume that, after a certain depth, zR within the pore network, such that, 2 < zR < M, the reaction zone (RZ) begins. In this study, we are interested only in the case of reactions that result in products that are solid and therefore precipitate in the pore space where the chemical reactions take place. Namely, we ignore the case where the solid formed from the reactions can migrate through the pores and block smaller pores located further downstream. Thus, once the nonwetting phase enters the RZ, solid formation can occur, which could block the pore where it occurs completely. Of importance to the process are the solid formation reaction rate and the characteristic time for reaction compared to the characteristic time for invasion. Alternatively, one could consider the case where the reactions can result in dissolving the original porous medium and therefore create channels of high porosity/permeability (“wormholes”). We examine the solid formation process parametrically in the following simple way. A dimensionless number λ, that is related to the reaction rate constant kr [LT−1], is defined. The dimensionless number λ is equal to “the number of invasion steps performed before one reaction step is performed”. The aforementioned definition produces integer values for λ (i.e., 2, 3, 4, etc.). A more generalized definition can be introduced as follows: for a given time period t,

In this work, trapping rules are considered76 (i.e., a phase is considered trapped if it is completely surrounded by a different phase). One can use reported experimental values of pores/ throats that correspond to oceanic sediments (see for example Tsimpanogiannis and Lichtner41 for the case of methane hydrate dissociation) to reconstruct a more realistic porous medium. However, in this work we consider that all pores have the same pore volume. In addition we consider that the throats have an invasion threshold πj given in the general case by πj = τj + Bzj

(6)

where B is the Bond number defined by eq 1, zj is the elevation of site j in the direction where gravity acts, and τj is a random number drawn from a uniform distribution on the unit interval. The introduction of a random term reflects the randomness of capillary thresholds corresponding to the throats, which is a result of their random size. This approach is followed in order to capture the disorder of the porous medium. The numerical simulation proceeds as follows: initially all the lattice sites are filled with the wetting phase, which is the defending phase during the process, while the nonwetting phase is the invading phase. Periodic boundary conditions are used in the directions perpendicular to the gravity gradient. A line injection scheme is considered along the entire top boundary of the network that is invaded by the nonwetting phase (i.e., in all pores/sites with coordinates (1, i), namely all sites at depth (zj = 1), the wetting phase is replaced by the nonwetting phase). Alternatively, one could also consider the point injection scheme, which would resemble closer the case of injecting through a well. In that case, we can choose a single, particular site [e.g., located at the center of the network, with coordinates (M/2, N/2)] which is initially invaded by the nonwetting phase. In either case, an interface develops within the pore network that separates sites that are occupied by the wetting phase (denoted as W-sites) from sites occupied by the nonwetting phase (denoted as NW-sites). A wetting site belongs to the interface (denoted as I-site) if at least one next-nearing neighbor site is occupied with nonwetting phase. As a first step in the process, all wetting sites that belong to that interface (invasion front) are identified along with the corresponding throats. Once the invasion front is identified, the throat (and the respective connected pore to that particular throat) with the least invasion threshold is selected to be invaded next. In the case Bzj > 0, the nonwetting phase will tend to grow toward throats with smaller zj, therefore limiting the width of the invading front to a narrow zone perpendicular to the direction of the imposed gradient. Eventually, invasion patterns that are consistent with IPSG are produced. In the case Bzj < 0, the nonwetting phase will tend to grow toward throats with larger zj. As a result, this will favor the invasion of throats that are further advanced in the direction parallel to the gravity vector. Therefore, finger-like invasion patterns are expected. Eventually, invasion patterns are produced that are consistent with IPDG. In this work, we are interested in the case of Bz < 0 (i.e., fingering invasion patterns). If we consider that zj > 0 for values of zj increasing downward, then there are two possible cases that can result in Bz < 0: (i) when ρdef > ρinv, therefore B > 0, and the invading-phase is injected upward (zj < 0) (e.g., methane gas moving upward toward the hydrate stability zone), and (ii) when ρdef < ρinv, (B < 0), and the invading-phase is injected downward zj > 0 (e.g., carbon dioxide injected downward into brine at depth larger than 3000 m). The

λ=

no. sites that invasion took place no. sites that reaction took place

(7)

In that case, noninteger values can be obtained [e.g., for λ = 1.5 we take λ = 3/2 (i.e., first three invasion steps occur, followed by two reaction steps)]. During the simulation of the reaction part, we first identify all pores saturated with nonwetting phase, which belong at the RZ, and then we select, at random, one of the pores to be occupied completely with solid. In this work, we identify two different solid-forming mechanisms (referred to as M-I, M-II), through which the decision is made where precipitation occurs, and they are based primarily on the location of the pore. In the case of mechanism M-I, all invaded pores that belong to the RZ are considered and one of them is chosen randomly to be occupied by solid. In case of mechanism M-II, from all invaded pores that belong to the RZ, those pores that have a first-neighboring pore occupied with wetting phase, are first identified. Then, a pore belonging to that subgroup of invaded pores in the RZ is chosen randomly to be occupied by solid. In the M-I case, the inherent assumption is that there is relatively fast interdiffusion between the two phases (compared to the reaction time). In the M-II case, interdiffusion is considered slower, and thus, there is preferential formation of solid at the wetting/nonwetting interface. In this work, emphasis is placed on the M-I case. As a simple example for demonstration purposes, consider the case where λ = 2. In that case, the nonwetting phase first invades two sites (pores). Then, one invaded site can be 3938


Energy & Fuels

Article

numbers (i.e., large λ values) correspond to slow reaction kinetics.

occupied with solid if the available invaded sites are present inside the RZ. Note that smaller λ implies faster reaction kinetics, while larger λ implies slower reaction kinetics. The upper limiting case when λ → ∞ implies that no reaction occurs in the system and therefore no solid formation and precipitation. The lower limiting case when λ = 1 implies that, in every pore in the RZ that is invaded, instantaneous reaction and solid formation and precipitation takes place. In that particular limiting case, at the end of the process the pores that have reacted and were blocked belong on a straight line with depth of a single pore site. Solid formation and precipitation inside a pore results in complete blocking of that pore and therefore cannot accommodate any flow after it is blocked. Consequently, the porosity and permeability of the system evolve as a function of time. The end of the simulation process is defined as either one of the following two cases: (i) during the nonwetting phase breakthrough at the opposite end (in case of small λ values), or (ii) when the invasion stops as a result of the complete blockage of a thin layer at the beginning of the reaction zone (in case of large λ values). In a process in which flow and reaction occurs, the following dimensionless parameters are utilized to describe the process. They are the Peclet number Pe =

UL D

RESULTS AND DISCUSSION Invasion Patterns. Figures 3−6 show typical invasion patterns obtained from 200 × 200 pore networks and for various different cases examined in this study. In the results shown, we consider the effects of both the Bond number B (B = 0.1, 0.01, 0.001, 0.0) and the reaction rate parameter λ (λ = 2,6,10,∞). The results range from very fast (small λ values, Figure 3), to moderate (medium λ values, Figure 4), to very slow (large λ values, Figure 5) solid formation and precipitation kinetics. Ultimately, when λ = ∞, no reaction occurs (Figure 6). Note that for cases with fast reaction kinetics, as the nw− phase enters the RZ, it reacts immediately and forms a solid phase, which subsequently precipitates. This process results in blocking the pore space. The production and precipitation of a solid phase results in the formation of a thin layer (with a width up to several lattice sites) of solid which does not allow any further penetration of the nw−phase into the RZ. As a result of the pore blocking, a pool of trapped nw−phase is formed above the RZ (e.g., Figure 3b−d, and Figure 4d). For fast or moderate reaction kinetics, the invading phase can reach the opposite end, namely, can percolate, only when the gradient is strong. In particular, the faster the kinetics are, the stronger the gradient needs to be for percolation to occur. This behavior is clearly shown in Table 1. For slow reaction kinetics the invading fluid phase percolates during all Bond values examined. For the case that no solid-formation reactions occur in the system, a single dominant finger appears in the presence of strong gravity gradients (e.g., Figure 6a−c). This finger has a characteristic width ξ that scales with the gravity Bond number as described by eq 5. However, in the case of chemical reactions with solid formation and precipitation, multiple fingers can develop as a result of pore blocking. Essentially, the growth of a finger can stop due to pore blocking and fluid invasion continues through the creation of a new finger (e.g., Figure 4a−c and Figure 5c). This process is repeated until the simulation is terminated. The analysis of multiple-finger development in such systems is currently under study and will be the subject of a future paper. Reaction Front. Figure 7a shows the effect of λ on the width σf of the reaction zone that has formed, and which is defined through eq 3−eq 4. The results shown in Figure 7 correspond to pore networks with lattice size 100 × 100 and four different Bond numbers (B = 0.1, 0.01, 0.001, 0.0). The results for the width σf of the reaction zone that has formed (for each Bond number) are obtained after averaging the resulting values from 100 realizations. Values for the standard deviations for all the various (λ, B)-systems considered are given in Table 2. For this series of simulations, the RZ begins at zj = zR = 10. As is shown in the Figure 7a, σf increases with λ (for a given Bond number) up to a critical value of λ where σf starts approaching a constant value. This behavior indicates that the invading fluid reached the opposite exit of the network and is a result of the finite size of the pore network used. Using larger networks is expected to result in shifting the critical value for λ to higher values. For a given value of λ, larger value of B implies deeper penetration of the invasion front and opportunity to react

(8)

and the first or second Damkohler number DaI =

■

kr DaII = U Pe

and

DaII =

k rL D

78

(9)

where U is the average steady-state velocity of the fluid over the entire domain, L is a characteristic length of the system, and D is the molecular diffusivity. The Pe number describes the effect of convection relative to that of molecular diffusion on solute transport and the DaI and DaII describe the effect of reaction relative to that of convection and diffusion, respectively. The dimensionless numbers, as defined in eqs 8 and 9, pertain to the global description of the system. In a similar manner, we can define the corresponding local (at the level of a single pore) dimensionless numbers where L is replaced by Lconv * and Lreact * , which are characteristic lengths for convection and reaction, respectively, at the local level. Next, we define two characteristic times, tconv and treact, as follows: tconv is the time required for the interface to advance to the unit distance (i.e., one site), and treact is the time required for the reaction and complete blocking of one site due to solid formation and precipitation. Note that tconv ∼ L*conv/U and treact ∼ L*react/kr. From the definition of λ, we have treact = λtconv

(10)

After substituting into eq 10 we obtain, ⎛ L * ⎞⎛ U ⎞ λ ∼ ⎜ react ⎟⎜ ⎟ * ⎠⎝ k r ⎠ ⎝ Lconv

(11)

and therefore, λ ∼ DaI−1 =

Pe DaII

(12)

Recall that large Damkohler numbers (i.e., small λ values) correspond to fast reaction kinetics and small Damkohler 3939


Energy & Fuels

Article

Figure 3. Effect of the Bond number, B, on the obtained invasion patterns of the nonwetting phase in a 200 × 200 lattice for λ = 2 and various B: (a) B = 0.1; (b) B = 0.01; (c) B = 0.001; and (d) B = 0.0. Color code: dark blue (nonwetting phase), light blue (wetting phase), red (trapped wetting phase), and yellow (solid). The yellow solid line denotes the beginning of the RZ (zR = 20).

Figure 4. Effect of the Bond number, B, on the obtained invasion patterns of the nonwetting phase in a 200 × 200 lattice for λ = 6 and various B: (a) B = 0.1; (b) B = 0.01; (c) B = 0.001; and (d) B = 0.0. Color code: dark blue (nonwetting phase), light blue (wetting phase), red (trapped wetting phase), and yellow (solid). The yellow solid line denotes the beginning of the RZ (zR = 20). 3940


Energy & Fuels

Article

Figure 5. Effect of the Bond number, B, on the obtained invasion patterns of the nonwetting phase in a 200 × 200 lattice for λ = 10 and various B: (a) B = 0.1; (b) B = 0.01; (c) B = 0.001; and (d) B = 0.0. Color code: dark blue (nonwetting phase), light blue (wetting phase), red (trapped wetting phase), and yellow (solid). The yellow solid line denotes the beginning of the RZ (zR = 20).

Figure 6. Effect of the Bond number, B, on the obtained invasion patterns of the nonwetting phase in a 200 × 200 lattice for λ = ∞ and various B: (a) B = 0.1; (b) B = 0.01; (c) B = 0.001; and (d) B = 0.0. Color code: dark blue (nonwetting phase), light blue (wetting phase), red (trapped wetting phase), and yellow (solid). The yellow solid line denotes the beginning of the RZ (zR = 20). 3941


Energy & Fuels

Article

Table 1. Type of Patterns for the Invading Phase at the End of the Simulation Inside 200 × 200 Pore Networks invasion pattern B

λ=2

λ = 3, 4, 5

λ = 6, 7, 8

λ = 9, 10

λ=∞

0.1 0.01 0.001 0.

percolated blocked blocked blocked

percolated percolated blocked blocked

percolated percolated percolated blocked

percolated percolated percolated percolated

percolated percolated percolated percolated

Figure 7. (a) Effect of λ on the width σf of the reaction zone that has formed. (b) Effect of λ on the mean position of the reaction zone. Results are averaged over 100 realizations and are for 100 × 100 pore networks and four different Bond numbers, B.

Table 2. Statistics for the Width of the Reaction Zone that Has Formed, σfa

position is limited to small values, located close to the beginning of the RZ. (ii) When the gravity gradient is strong (i.e., large B), the invasion front can eventually percolate and the reaction zone extends to the entire network within the RZ with higher solid saturation close to the beginning of the RZ. At large λ (i.e., slow kinetics), the following two cases exist: (i) When the gravity gradient is weak (i.e., small B), the invasion front can percolate, and the limiting value corresponds to

standard deviation for σf B

λ=2

λ=3

λ=4

λ=5

λ=6

λ=7

λ=8

λ=9

0.1 0.01 0.001 0.

6.25 2.55 1.80 1.57

1.87 4.84 3.90 3.80

1.83 4.05 5.46 4.76

1.85 3.84 5.58 5.99

1.51 3.40 5.81 5.51

1.57 2.08 5.28 6.09

1.50 2.01 5.33 5.69

1.48 2.25 4.41 5.87

Obtained from averaging the results of 100 realizations using 100 × 100-sized pore networks.

a

X ft = z R +

deeper inside the pore network and, therefore, larger σf values, as clearly shown in Figure 7a. The effect of λ on the mean position of the reaction zone, Xft, defined through eq 3, is shown in Figure 7b. Values for the standard deviations for all the various (λ, B)-systems considered are given in Table 3. At small λ (i.e., fast kinetics), the following two cases exist: (i) When the gravity gradient is weak (i.e., small B), the invasion process is interrupted early, as a result of pore clogging limited to a narrow zone that does not allow any further fluid invasion. Therefore, the mean reaction front

⎛ N − zR ⎞ ⎜ ⎟ ⎝ 2 ⎠

(13)

where zR is the location where the RZ begins and N is the size of the square pore network. For the results shown in Figure 7b the following values apply: zR = 10 and N = 100, therefore, Xft = 55 (approximately in the middle of the pore network). (ii) When the gravity gradient is strong (i.e., large B), the invasion front can percolate following fingering-type patterns and the mean position of the reaction zone, Xft, has values less than N/2 (e.g., ≈40). As a result of the strong gradient, the invading fluid percolates soon after the process begins and therefore there is not enough

Table 3. Statistics for the Mean Position of the Reaction Zone, Xfta standard deviation for Xft

a

B

λ=2

λ=3

λ=4

λ=5

λ=6

λ=7

λ=8

λ=9

0.1 0.01 0.001 0.

9.18 4.35 3.29 3.48

3.71 8.59 8.59 9.15

3.56 6.82 10.94 12.01

3.69 5.63 10.16 11.17

2.75 6.68 10.55 11.85

3.11 4.66 8.29 10.80

3.15 4.60 9.89 9.98

2.47 4.62 8.47 11.51

Obtained from averaging the results of 100 realizations using 100 × 100-sized pore networks. 3942


Energy & Fuels

Article

Figure 8. Effect of the Bond number, B, on the overall averaged saturations as a function of λ: (a) wetting phase saturation; (b) trapped wetting phase saturation; (c) nonwetting phase saturation; and (d) produced solid phase saturation. Results are averaged over 100 realizations and are for 100 × 100 pore networks.

Table 4. Statistics for the Produced Solid Phase Saturation, Ssa standard deviation for Ss

a

B

λ=2

λ=3

λ=4

λ=5

λ=6

λ=7

λ=8

λ=9

0.1 0.01 0.001 0.

0.024 0.013 0.010 0.011

0.009 0.016 0.017 0.018

0.006 0.010 0.016 0.017

0.004 0.008 0.012 0.012

0.003 0.007 0.010 0.011

0.002 0.006 0.007 0.009

0.002 0.005 0.007 0.007

0.001 0.004 0.006 0.008


Table 5. Statistics for the Wetting Phase Saturation, Swa standard deviation for Sw

a

B

λ=2

λ=3

λ=4

λ=5

λ=6

λ=7

λ=8

λ=9

0.1 0.01 0.001 0.

0.073 0.038 0.032 0.034

0.045 0.081 0.082 0.090

0.042 0.096 0.119 0.124

0.047 0.095 0.123 0.122

0.027 0.100 0.137 0.135

0.027 0.089 0.119 0.132

0.029 0.079 0.138 0.130

0.022 0.072 0.129 0.154


saturation, Sw, Figure 8b shows the trapped wetting phase saturation, Swt, Figure 8c shows the nonwetting phase saturation, Snw, and Figure 8d shows the produced solid phase saturation, Ss. The various overall averaged saturations are measured at the end of the process, as it was defined in the previous section. Values for the standard deviations of all

time for solid-formation reactions and subsequent precipitation of the solid phase to occur close to the opposite end. Overall Averaging. The effect of the Bond number, B, on the overall averaged saturations as a function of λ is shown in Figure 8. In particular, Figure 8a shows the wetting phase 3943


Energy & Fuels

Article

Table 6. Statistics for the Trapped Wetting Phase Saturation, Swta standard deviation for Swt

a

B

λ=2

λ=3

λ=4

λ=5

λ=6

λ=7

λ=8

λ=9

0.1 0.01 0.001 0.

0.031 0.012 0.011 0.011

0.021 0.039 0.034 0.038

0.022 0.063 0.060 0.063

0.029 0.062 0.070 0.065

0.015 0.064 0.083 0.075

0.014 0.056 0.077 0.079

0.016 0.045 0.091 0.083

0.012 0.048 0.082 0.087


Table 7. Statistics for the Non-Wetting Phase Saturation, Snwa standard deviation for Snw

a

B

λ=2

λ=3

λ=4

λ=5

λ=6

λ=7

λ=8

λ=9

0.1 0.01 0.001 0.

0.023 0.014 0.012 0.013

0.019 0.032 0.035 0.038

0.018 0.031 0.049 0.051

0.017 0.034 0.049 0.051

0.013 0.038 0.052 0.058

0.013 0.036 0.042 0.055

0.014 0.035 0.050 0.051

0.011 0.032 0.048 0.066


Figure 9. (a) Effect of λ on the overall available pore-space of the system. (b) Effect of λ on the overall available space for flow of the non wetting phase. Results are averaged over 100 realizations and are for 100 × 100 pore networks and four different Bond numbers, B.

Transverse Averaging. When studying systems where a reacting front propagates through the system, thus affecting the porosity/permeability of the system, it is useful to put emphasis on the study of phenomena occurring at the vicinity of the reacting front. Thus, it is useful to define transverse averages (see also Appendix A for definitions). Figure 10 shows the transversely averaged (over 100 realizations) pore-space available for nw-phase flow as a function of lattice position, z, for a pore network. An important observation from Figures 9 and 10 is that even though the total porosity of the network does not drop significantly due to solid formation and precipitation (ϕI/ϕo in the range 0.92−0.98), this drop occurs in a narrow zone, located close to the beginning of the RZ, and thus blocks the pore network, inhibiting any further fluid invasion. Smaller values of λ result to narrower zones. Results obtained in this study are in qualitative agreement with the results from Lattice Boltzmann simulations in fractures, reported by Kang et al.9 For λ = 2, at the beginning of the reaction zone (i.e., zR = 10), less than 10% of the pores are open to fluid invasion at the end of the simulation. The remaining pores are completely blocked as a result of solid formation and precipitation. This behavior is repeated to the subsequent transverse layers (having higher % of open pores, however). The combination of several

the averaged saturations that were considered are given in Tables 4−7. Figure 9 shows the changes in the overall available porespace of the system. In particular, Figure 9a depicts the overall change in porosity defined as ϕI/ϕo = 1 − Ss, where ϕo is the initial porosity. Essentially, in this calculation, we subtract the newly formed solid fraction from the initial porosity. Figure 9b depicts the overall change in the available space for flow of the non wetting phase, defined as ϕII/ϕo = 1 − Ss − Swt = ϕI/ϕo − Swt. Namely, from the overall available pore-space of the system, we subtract the pore-space occupied by the trapped wetting phase, since we ignore any mobilization of the trapped phase. Both ϕI/ϕo and ϕII/ϕo exhibit similar behavior as a function of λ. In the presence of strong gravity gradients, the curves are monotonically increasing as λ increases (i.e., less pore space is blocked and therefore more space is available to flow). As gradients become weaker, the curves initially decrease, and after going through a minimum value, they increase as λ increases. The first (decreasing) part of the curve is a result of the complete blocking of a thin layer of the invasion process that results in the premature stopping of the process (see also Table 1). 3944


Energy & Fuels

Article

Figure 10. Effect of the Bond number, B, on the transversely averaged available porosity for the invading fluid to flow as a function of the lattice position, j, for a 100 × 100 pore network (averaged over 100 realizations).

from invasion percolation. The model combines aspects of fluid displacement and chemical reaction accompanied with solid precipitation inside the porous medium. A number of simplifications was considered during the development of the model, which result in certain model limitations. Relaxing some of the assumptions can result in obtaining more realistic behavior for the model. However, such an approach would require a significantly higher computational effort, which then would limit the size of systems that can be studied. Various modifications of invasion percolation have been used extensively to study immiscible displacement in porous media. The modification used in this work can not be used to study miscible fluid displacement. To this purpose, a more complicated approach should be taken. Instead of using random thresholds for pores/throats, as done here, one could use experimental information for pore/throat size distributions in order to reconstruct a more realistic porous medium. Furthermore, pore/throat geometries can be used that include the presence of corners, resulting thus in producing more realistic flow fields. By applying a pressure gradient at the opposite ends of the model, the flow field and the pressure distribution can be obtained within the porous medium. Such an approach can be used to study the case of miscible flow. However, this approach is more computationally demanding (i.e., the pressure field needs to be calculated every time there is a change in the fluid distribution) when compared to applying

such layers results in the complete stopping of the invasion as a result of pore blocking. For higher λ values, less pore blockage occurs at each transverse layer, therefore, wider zones are required to achieve the formation of such blocked-pore distributions that are capable to stop the invasion process. Figures 11 and 12 show the transversely averaged saturations (denoted by subscript T) as a function of lattice position (depth of invasion), z, for a pore network. Averaging is performed over 100 realizations. Figure 11 shows the results for the transversely averaged solid saturation, ⟨Ss⟩T, and as expected, higher solid saturations are observed for smaller λ values and closer to the beginning of the RZ. Results for the transversely averaged nw-phase saturation, ⟨Snw⟩ T are shown in Figure 12. Note that for the case of λ = 2 the nw-phase saturation drops to almost zero after invading 10−11 lattice units from the beginning of the RZ (recall that the RZ begins at zR = 10). That behavior clearly shows that the process ends, due to solid blocking, before the invading nw−phase front reaches the opposite end at z = N. In order to observe similar behavior for higher values of λ, larger lattice sizes are required. The lattice position where the RZ begins is located at z = 10 and is denoted by the vertical line in the figures. Model Limitations and Future Outlook. In this study, we developed a simulation model that was based on principles 3945


Energy & Fuels

Article

Figure 11. Effect of the Bond number, B, on the transversely averaged saturations for the solid phase formed as a function of the lattice position, j, for a 100 × 100 pore network (averaged over 100 realizations).

blocked at every step. Furthermore, in a real porous system partial blocking of pores/throats can occur as a result of solid precipitation. A more detailed simulator, which also solves for the concentration field, could capture more accurately such behavior. Such capability is not available at the current model which follows the rule of “complete blockage” of a single pore body at each reaction step. In a real porous medium solid precipitation can happen under various conditions. Namely, precipitation results from reactions, either occurring during two-phase immiscible flow, or occurring during single-phase miscible conditions. Of relevance to the problem under consideration are the recent studies of Tartakovsky et al.,79 Zhang et al.,80 and Ghezzehei et al.81 In the current study, two mechanisms were identified (i.e., M-I, M-II) for selecting the pore, where precipitation occurred. Mechanism M-I considers all invaded (by the nonwetting phase) pores in the reaction zone as candidates for precipitation. This approach can be seen as a primitive way to describe precipitation that occurs under miscible conditions. On the other hand, mechanism M-II considers only the invaded pores that are at the wetting/nonwetting interface and are also in the reaction zone, as candidates for precipitation. Namely, it corresponds to the case of precipitation under immiscible conditions. Earlier experimental and numerical studies61,62,82,83 have shown that the invading phase, under the influence of a

the simple invasion percolation rules as implemented in this study. The chemical reaction in the porous medium that results in solid formation and subsequent precipitation within the pores is described using a very simple mechanistic model in this work. This approach is mainly driven by the objective to qualitatively study the interplay between fluid displacement and reaction in a very broad range of conditions including the presence/absence of gradients affecting the displacement and the rate of chemical reactions (e.g., fast or slow chemical kinetics). In order to obtain a more quantitative result, the introduction of a kinetic rate law and an equilibrium constant are required. For example, a first order linear rate form such as R = k(c − ceq), where k is the local rate constant, c is the concentration, and the subscript eq denotes equilibrium value, could be used. In this approach, the concentration field needs to be calculated in addition to the pressure field. The mechanistic model for chemical reactions, which was used in the current study, considers that precipitation occurs only inside the pores. Alternatively, precipitation can occur inside the throats only or inside both pores and throats. Limiting the precipitation to pores only implies that when a pore is blocked at every step, as a result of precipitation, all the next-nearing throats are also excluded from flow. Therefore, the invasion process can be interrupted earlier due to complete blockage, compared to the case when only a single throat is 3946


Energy & Fuels

Article

Figure 12. Effect of the Bond number, B, on the transversely averaged saturations for the invading fluid phase as a function of the lattice position, j, for a 100 × 100 pore network (averaged over 100 realizations).

gradient, can fragment into several “blobs” that can migrate independently from the main part of the invading phase. The phenomenon of the creation of the disconnected parts occurs as a result of fluid “snap-offs”. This process can be repeated until the blobs become smaller than a critical size to be mobilized. The current model does not capture the “snap-off” mechanism. However, such a mechanism, if included in the simulations, could have an effect on the results, particularly when considering the case of mechanism M-II for reactions. The “snap-off” mechanism results in creating additional wetting/nonwetting interfaces and, therefore, could affect the precipitation patterns inside the reaction zone. The current study is focused only in 2-D geometries. The effect of considering 3-D geometries on the results needs to be examined as well. For example, in 3-D geometries, the effect of the wetting-phase trapping during the displacement is less significant than in the case of 2-D geometries. This could have an effect on the regions were solid formation and precipitation occurs.

and precipitation) occurring when the nonwetting phase entered the reaction zone. Even though the developed pore network model used simplistic rules/mechanisms to describe the displacement and reactive processes, the model was able to capture, at least qualitatively, the rich interplay between fluid displacement and chemical reaction in a wide range of conditions, including the presence of strong or weak gravity gradients, and the fast or slow chemical reaction rates. Comparing the simulation results with appropriate experiments would be a desirable future objective. There is significant room for improvement of the pore network simulator and a number of suggestions was offered, focusing especially on enhancing the quantitative aspects of the simulator. However, following such a path would result in increasing significantly the computational effort. In this study, we examined the obtained fluid invasion patterns, the solid deposition patterns, and the change in porosity due to solid formation and how they depend on various parameters including the solid formation kinetics and the strength of the imposed gravity gradient. A rich variety of solid deposition patterns was observed that ranged from compact to single-finger or multiple-finger dominated patterns. The percolation of the nonwetting invading phase was found to depend on the strength of the gradients and the reaction kinetics. Fast kinetics when coupled with weak gravity gradients acted against the percolation of the invading phase. On the

■

CONCLUSION A 2-D pore network model based on concepts from invasion percolation in a gradient was developed to simulate the invasion of nonwetting phase into porous media and the subsequent immiscible displacement of the wetting phase originally present in the porous medium. The invasion process was also coupled with a reaction process (i.e., solid formation 3947


Energy & Fuels

Article

Laboratory, Los Alamos, NM. Partial financial support by LDRD-DR 20030059DR and LDRD-DR 20040042DR projects, funded by Los Alamos National Laboratory, is gratefully acknowledged. Partial financial support by the European Commission EC Grant PERL (Contract No. REGPOT-20081-229773) is gratefully acknowledged.

other hand, slow kinetics when coupled with medium/strong gravity gradients acted in favor of the percolation of the invading phase. We performed a parametric study of the formation kinetics and it was shown that in the case of fast reaction kinetics the local porosity at the vicinity of the invading front decreases significantly, due to solid formation and precipitation, essentially inhibiting any further invasion of the nonwetting phase into the porous medium. In the case of slower kinetics, further penetration of the nonwetting phase was possible. It was clearly shown that during the presence of fast reaction kinetics the transversely averaged porosity has a more pronounced effect on the transport properties of the pore network. Namely, even very small changes in the overall porosity may have a very detrimental effect on the percolation ability of the pore network, especially when the changes in the porosity, due to pore blocking, are focused in a very narrow zone. In that case, the transversely averaged porosity changes very drastically.

■

■

APPENDIX A Here, we introduce the definition of certain parameters used in the discussion of the results. We denote with subscript p the different phases present in the pore network. We further differentiate between nontrapped water (W) and trapped water (WT). So, in general, we can define the overall average (denoted by ⟨Sp⟩) and the jth-row transverse average (denoted by ⟨Sp⟩T(j)) for the saturation of phase p, respectively, as follows: j=M

⟨Sp⟩ =

i=N

∑j = 1 ∑i = 1 site(i , j) j=M

=

i=N

∑j = 1 ∑i = 1 Wp(i , j)site(i , j) j=M

i=N

∑j = 1 ∑i = 1 Wp(i , j)site(i , j) M×N

(14)

and i=N

⟨Sp⟩T (j) =

∑i = 1 Wp(i , j) site(i , j) i=N

∑i = 1 site(i , j) i=N

=

∑i = 1 Wp(i , j) site(i , j) N

(15)

where site(i,j) = 1 for all sites that belong to the network and Wp(i,j) is a phase function defined as follows:

■

⎧1 Wp(i , j) = ⎨ ⎩0

if site(i , j) contains phase p otherwise

REFERENCES

(1) Bear, J. Hydraulics of Groundwater; McGraw Hill: New York, 1979. (2) Dake, L. P. The Practice of Reservoir Engineering; Elsevier: Amsterdam, 1994. (3) Bedient, P. B.; Rifai, H. S.; Newell, C. J. Ground Water Contamination: Transport and Remediation; Prentice-Hall: Englewood Cliffs, NJ, 1994. (4) Alvarado, V.; Davis, H. T.; Scriven, L. E. Effects of pore-level reaction on dispersion in porous media. Chem. Eng. Sci. 1997, 52, 2865−2881. (5) Rieckmann, C.; Keil, F. J. Multicomponent diffusion and reaction in three-dimensional networks: General kinetics. Ind. Eng. Chem. Res. 1997, 36, 3275−3281. (6) Steefel, C. I.; DePaolo, D. J.; Lichtner, P. C. Reactive transport modeling: An essential tool and a new research approach for the earth sciences. Earth Planet. Sci. Lett. 2005, 240, 539−558. (7) Ortoleva, P.; Chadam, J.; Merino, E.; Sen, A. Geochemical selforganization II. The reactive infiltration instability. Am. J. Sci. 1987, 287, 1008−1040. (8) Daccord, G. Chemical dissolution of a porous medium by a reactive fluid. Phys. Rev. Lett. 1987, 58, 479−482. (9) Kang, Q.; Tsimpanogiannis, I. N.; Zhang, D.; Lichtner, P. C. Numerical modeling of pore-scale phenomena during CO2 sequestration in oceanic sediments. Fuel Process. Technol. 2005, 86, 1647−1665. (10) IPCC Special Report: Carbon Dioxide Capture and Storage; Cambridge University Press: New York, 2005. (11) Herzog, H. J.; Drake, E. M. Carbon dioxide recovery and disposal from large energy systems. Annu. Rev. Energy Environ. 1996, 21, 145−166. (12) Golfier, F.; Zarcone, C.; Bazin, B.; Lenormand, R.; Lasseux, D.; Quintard, M. On the ability of a Darcy-scale model to capture wormhole formation during the dissolution of a porous medium. J. Fluid Mech. 2002, 457, 213−254. (13) Hyndman, R. D.; Davies, E. E. A mechanism for the formation of methane hydrate and seafloor bottom-simulating reflectors by vertical fluid expansion. J. Geophys. Res. B 1992, 97, 6683. (14) Sloan, E. D. Clathrate Hydrates of Natural Gases, 2nd ed.; Marcel Dekker: New York, 1998. (15) Sloan, E. D. Fundamental principles and applications of natural gas hydrates. Nature 2003, 426, 353−359. (16) Fatt, I. The network model of porous media. I. Capillary pressure characteristics. Trans. AIME 1956, 207, 144−159. (17) Feder, J. Fractals; Plenum; New York, 1988. (18) Sahimi, M. Flow phenomena in rocks: from continuum models to fractals, percolation, cellular automata, and simulated annealing. Rev. Mod. Phys. 1993, 65, 1393−1534. (19) Sahimi, M. Flow and Transport in Porous Media and Fractured Rock: From Classical Methods to Modern Approach; VHC: Weinheim, 1995. (20) Celia, M. A.; Reeves, P. C.; Ferrand, L. A. Recent advances in pore scale models for multiphase flow in porous media. Rev. Geophys. 1995, 33 (S1), 1049−1058. (21) Berkowitz, B.; Ewing, R. P. Percolation theory and network modeling applications in soil physics. Surveys Geophys. 1998, 19, 23− 72. (22) Blunt, M. J. Effects of heterogeneity and wetting on relative permeability using pore level modeling. SPE J. 1997, 2, 70−87. (23) Blunt, M. J. Flow in media − pore-network models and multiphase flow. Curr. Opin. Colloid Interface Sci. 2001, 6, 197−207.

(16)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. § Visiting Scientist.

■

ACKNOWLEDGMENTS Part of this work was completed when I.N.T. was at the Hydrology, Geochemistry, and Geology Group (EES-6), Earth and Environmental Sciences Division, Los Alamos National 3948


Energy & Fuels

Article

(24) Blunt, M. J.; Jackson, M. D.; Piri, M.; Valvante, P. H. Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow. Adv. Water Resour. 2002, 25, 1069−1089. (25) Singh, M.; Mohanty, K. K. Dynamic modeling of drainage through three-dimensional porous materials. Chem. Eng. Sci. 2003, 58, 1−18. (26) Rege, S. D.; Fogler, H. S. A network model for deep bed filtration of solid particles and emulsion drops. AIChE J. 1988, 34, 1761−1772. (27) Imdakm, A. O.; Sahimi, M. Computer simulation of particle transport processes in flow through porous media. Chem. Eng. Sci. 1991, 46, 1977−1993. (28) Li, X.; Yortsos, Y. C. Visualization and simulation of bubble growth in pore networks. AIChE J. 1995, 41, 214−222. (29) Dominguez, A.; Bories, S.; Prat, M. Gas cluster growth by solute diffusion in porous media. Int. J. Multiphase Flow 2001, 26, 1951− 1979. (30) Prat, M. Isothermal drying of non-hygroscopic capillary porous materials as an invasion percolation process. Int. J. Multiphase Flow 1995, 21, 875−892. (31) Prat, M.; Bouleux, F. Drying of capillary porous media with a stabilized front in two dimensions. Phys. Rev. E 1999, 60, 5647−5656. (32) Yiotis, A. G.; Stubos, A. K.; Boudouvis, A. G.; Yortsos, Y. C. A 2D pore-network model of the drying of single-component liquids in porous media. Adv. Water Resour. 2001, 24, 439−460. (33) Yiotis, A. G.; Boudouvis, A. G.; Stubos, A. K.; Tsimpanogiannis, I. N.; Yortsos, Y. C. Effect of liquid films on the isothermal drying of porous media. Phys. Rev. E 2003, 68, 037303. (34) Yiotis, A. G.; Tsimpanogiannis, I. N.; Stubos, A. K.; Yortsos, Y. C. Pore-network study of the characteristic periods in the drying of porous materials. J. Colloid Interface Sci. 2006, 297, 738−748. (35) Yiotis, A. G.; Tsimpanogiannis, I. N.; Stubos, A. K.; Yortsos, Y. C. Coupling between external and internal mass transfer during drying of a porous medium. Water Resour. Res. 2007, 43, W06403. (36) Satik, C.; Yortsos, Y. C. A pore network study of bubble growth in porous media driven by heat transfer. ASME J. Heat Transfer 1996, 118, 455−562. (37) Du, C.; Yortsos, Y. C. A numerical study of the critical gas saturation. Transp. Porous Media 1999, 35, 205−225. (38) Mc Dougall, S. R.; Sorbie, K. S. Estimation of critical gas saturation during pressure depletion in virgin and waterflooded reservoirs. Pet. Geosci. 1999, 5, 229−233. (39) Tsimpanogiannis, I. N.; Yortsos, Y. C. The critical gas saturation in a porous medium in the presence of gravity. J. Colloid Interface Sci. 2004, 270, 388−395. (40) Wang, X.; Mohanty, K. K. Critical condensate saturation in porous media. J. Colloid Interface Sci. 1999, 214, 416−426. (41) Tsimpanogiannis, I. N.; Lichtner, P. C. Pore-network study of methane hydrate dissociation. Phys. Rev. E 2006, 74, 056303. (42) Liang, H.; Song, Y.; Liu, Y.; Yang, M.; Huang, X. Study of the permeability characteristics of porous media with methane hydrate by pore network model. J. Nat. Gas Chem. 2010, 19, 255−260. (43) Jang, J.; Narsilio, G. A.; Santamarina, J. C. Hydraulic conductivity in spatially varying media. A pore-scale investigation. Geophys. J. Int. 2011, 184, 1167−1179. (44) Hoefner, M. L.; Fogler, H. S. Pore evolution and channel formation during flow and reaction in porous media. AIChE J. 1988, 34, 45−54. (45) Rege, S. D.; Fogler, H. S. Competition among flow, dissolution, and precipitation in porous media. AIChE J. 1989, 35, 1177−1185. (46) Thompson, K. E.; Fogler, H. S. Pore-level mechanisms for altering multiphase permeability with gels. SPE J. 1997, 2, 350−362. (47) Fredd, C. N.; Fogler, H. S. Influence of transport and reaction on wormhole formation in porous media. AIChE J. 1998, 44, 1933− 1949. (48) Li, L.; Peters, C. A.; Celia, M. A. Upscaling geochemical reaction rates using pore-scale network modeling. Adv. Water Resour. 2005, 29, 1351−1370.

(49) Li, L.; Peters, C. A.; Celia, M. A. Effects of mineral spatial distribution on reaction rates in porous media. Water Resour. Res. 2007, 43, W01419. (50) Li, L.; Peters, C. A.; Celia, M. A. Applicability of averaged concentrations in determining geochemical reaction rates in heterogeneous porous media. Am. J. Sci. 2007, 307, 1146−1166. (51) Wilkinson, D.; Willemsen, J. F. Invasion Percolation: A new form of percolation theory. J. Phys. A: Math. Gen. 1983, 16, 3365− 3376. (52) Ferer, M.; Ji, C.; Bromhal, G. S.; Cook, J.; Ahmadi, G.; Smith, D. H. Crossover from capillary fingering to viscous fingering for immiscible unstable flow: Experiment and modeling. Phys. Rev. E 2004, 70, 016303. (53) Stauffer, D.; Aharony, A. Introduction to Percolation Theory; Taylor & Francis: London, 1992. (54) Wilkinson, D. Percolation model of immiscible displacement in the presence of buoyancy forces. Phys. Rev. A 1984, 30, 520−531. (55) Sapoval, B.; Rosso, M.; Gouyet, J. F. The fractal nature of the diffusion front and the relation to percolation. J. Phys. (France) Lett. 1985, 46, L149−L156. (56) Rosso, M.; Gouyet, J. F.; Sapoval, B. Gradient percolation in three dimensions and relation to diffusion fronts. Phys. Rev. Lett. 1986, 57, 3195−3198. (57) Gouyet, J. F.; Rosso, M.; Sapoval, B. Fractal structure of diffusion and invasion fronts in three-dimensional lattices through the gradient percolation approach. Phys. Rev. B 1988, 37, 1832−1838. (58) Clement, E.; Baudet, C.; Hulin, J. P. Multiple scale structure on nonwetting fluid invasion fronts in 3D model porous media. J. Phys. (France) Lett. 1985, 46, L1163−L1171. (59) Clement, E.; Baudet, C.; Guyon, E.; Hulin, J. P. Invasion front structure in a 3D model porous medium under a hydrostatic pressure gradient. J. Phys. D: Appl. Phys. 1987, 20, 608−615. (60) Hulin, J. P.; Clement, E.; Baudet, C.; Gouyet, J. F.; Rosso, M. Quantitative analysis of an invading-fluid invasion front under gravity. Phys. Rev. Lett. 1988, 61, 333−336. (61) Meakin, P.; Feder, J.; Frette, V.; Jossang, T. Invasion percolation in a destabilizing gradient. Phys. Rev. A 1992, 46, 3357−3368. (62) Frette, V.; Feder, J.; Jossang, T.; Meakin, P. Buoyancy-driven fluid migration in porous media. Phys. Rev. Lett. 1992, 68, 3164−3167. (63) Yortsos, Y. C.; Xu, B.; Salin, D. Phase diagram of fully developed drainage in porous media. Phys. Rev. Lett. 1997, 79, 4581−4584. (64) Xu, B.; Yortsos, Y. C.; Salin, D. Invasion percolation with viscous forces. Phys. Rev. E 1998, 57, 739−751. (65) Bretherton, F. P. The motion of long bubbles in tubes. J. Fluid Mech. 1961, 10, 166−188. (66) Dong, M.; Chatzis, I. The imbibition and flow of a wetting liquid along the corners of a square capillary tube. J. Colloid Interface Sci. 1995, 172, 278−288. (67) Dillard, L. A.; Blunt, M. J. Development of a pore network simulation model to study nonaqueous phase liquid dissolution. Water Resour. Res. 2000, 36, 439−454. (68) Zhao, W.; Ioannidis, M. A. Effect of NAPL film stability on the dissolution of residual wetting NAPL in porous media: A pore-scale modeling study. Adv. Water Resour. 2007, 30, 171−181. (69) Tsakiroglou, C. D.; Bourganos, V. N.; Jacobsen, J. Porestructure analysis by using nitrogen sorption and mercury intrusion data. AIChE J. 2004, 50, 489−510. (70) Reeves, P. C.; Celia, M. A. A functional relationship between capillary pressure, saturation, and interfacial area as revealed by a porescale network model. Water Resour. Res. 1996, 32, 2345−2358. (71) Acharya, R. C.; Van der Zee, S. E. A. T. M.; Leijnse, A. Approaches for modeling longitudinal dispersion in pore-networks. Adv. Water Resour. 2007, 30, 261−272. (72) Yiotis, A. G.; Boudouvis, A. G.; Stubos, A. K.; Tsimpanogiannis, I. N.; Yortsos, Y. C. Effect of liquid films on the drying of porous media. AIChE J. 2004, 50, 2721−2737. (73) Prat, M. On the influence of pore shape, contact angle and film flows on drying of capillary porous media. Int. J. Heat Mass Transfer 2007, 50, 1455−1468. 3949


Energy & Fuels

Article

(74) Tsimpanogiannis, I. N.; Yortsos, Y. C.; Poulou, S.; Kanellopoulos, N.; Stubos, A. K. Scaling theory of drying in porous media. Phys. Rev. E 1999, 59, 4353−4365. (75) Chauvet, F.; Duru, P.; Geoffroy, S.; Prat, M. Three periods of drying of a single square capillary tube. Phys. Rev. Lett. 2009, 103, 124502. (76) Dias, M. M.; Wilkinson, D. Percolation with trapping. J. Phys. A: Math. Gen. 1986, 19, 3131−3146. (77) Liro, C. R.; Adams, E. E.; Herzog, H. J. Modeling the release of CO2 in deep ocean. Energy Convers. Management 1992, 33, 667−674. (78) Damkohler, G. Z. Elektrochem. Phys. Chem. 1936, 42, 846−862. (79) Tartakovsky, A. M.; Redden, G.; Lichtner, P. C.; Scheibe, T. D.; Meakin, P. Mixing-induced precipitation: Experimental study and multiscale numerical analysis. Water Resour. Res. 2008, 44, W06S04. (80) Zhang, C.; Dehoff, K.; Hess, N.; Oostrom, M.; Wietsma, T. W.; Valocchi, A.; Fouke, B. W.; Werth, C. J. Pore-scale study of transverse mixing induced CaCO3 precipitation and permeability reduction in a model subsurface sedimentary system. Environ. Sci. Technol. 2010, 44, 7833−7838. (81) Ghezzehei, T. A. Linking sub-pore scale heterogeneity of biological and geochemical deposits with change in permeability. Adv. Water Resour. 2012, 39, 1−6. (82) Wagner, G.; Birovljev, A.; Meakin, P.; Feder, J.; Jossang, T. Fragmentation and migration of invasion percolation clusters: Experiments and simulations. Phys. Rev. E 1997, 55, 7015−7029. (83) Vedvik, A.; Wagner, G.; Oxaal, U.; Feder, J.; Meakin, P.; Jossang, T. Fragmentation transition for invasion percolation in hydraulic gradients. Phys. Rev. Lett. 1998, 80, 3065−3068.

3950


Fluid Displacement and Solid Formation in a Porous Medium Using

Recommend Documents