Next Article in Journal
Optimal Transient Growth in an Incompressible Flow past a Backward-Slanted Step
Previous Article in Journal
DNS Study of the Bending Effect Due to Smoothing Mechanism
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Ionic Fracture Fluid Leak-Off

1
Lavrentyev Institute of Hydrodynamics, Novosibirsk State University, Novosibirsk 630090, Russia
2
Trofimuk Institute of Oil and Gas Geology and Geophysics, Novosibirsk State University, Novosibirsk 630090, Russia
*
Author to whom correspondence should be addressed.
Fluids 2019, 4(1), 32; https://doi.org/10.3390/fluids4010032
Submission received: 10 January 2019 / Revised: 5 February 2019 / Accepted: 15 February 2019 / Published: 19 February 2019

Abstract

:
The study is motivated by monitoring the space orientation of a hydrolic fracture used in oil production. Streaming potential arises due to the leakage of ionic fracking fluid under the rock elastic forces which make the fracture disclosure disappear after pumping stops. The vector of electric field correlates with the fracture space orientation since the fluid leakage is directed normally to the fracture surfaces. We develop a mathematical model for the numerical evaluation of the streaming potential magnitude. To this end, we perform an asymptotic analysis taking advantage of scale separation between the fracture disclosure and its length. The contrast between the virgin rock fluid and the fluid invading from the fracture is proved to be crucial in a build up of a net charge at the invasion front. Calculations reveal that an increase of the viscosity and resistivity contrast parameters results in an increase of the streaming potential magnitude. Such a conclusion agrees with laboratory experiments.

1. Introduction

In applying hydraulic fracturing treatment of rocks, one should prepare for execting the size of a fracture and its space orientation. One of the reasons is to avoid hydraulic connectivity of a production well with an injection one, since during water-flood of an initially oil-filled reservoir, encroaching water spoils hydrocarbon recovery. Directional fracturing is also of importance in petrothermal power production while designing geothermal circulation systems within the dry rocks [1]. Typical circulation systems consist of two wells connected by a fracture. By injecting water into the hot, deep, crystalline rocks, a huge amount of heat can be harnessed [2].
The pressure transient data during the injection fall-off test provide a way to determine the dimensions of an induced fracture but not its direction [3]. Microseismic measurements allow us to find the direction of the fracture [4]. However, the reliability of such measurements is not evident due to the extremely low energy of microseismic events and a high level of acoustic noise. Therefore, alternative approaches are of interest.
In this study, we estimate the electric field induced near the hydrofracture by leakage of a contrast fluid from the fracture under elastic forces in the rock which make the fracture disclosure disappear after pumping stops. The electric potential, known as the streaming potential, arises due to electrokinetic phenomena when ionic fluid moves through a porous rock [5]. The vector of the induced electric field is directed normally to the fracture surfaces, while the electric potential at the fracture surfaces strongly depends on the location of the invasion front of the fracture fluid. Polarization occurs due to charge concentration at the invasion front. The charge density can be evaluated via the resistivity jump across the front.
Streaming potential measurement is based on response (DC voltage) associated with the excess of electrical charge that is transported by flow of conductive fluids through porous media under a pressure gradient or an electrical conductivity change. Such a method presents a key electrokinetic mechanism to study fluid flows in porous media. An electrical double layer (EDL) forms at the rock–fluid interfaces. The onset of fluid flow perturbs the distribution of the ions present in that layer and induces a current, which is the source of the streaming potential. This electrokinetic mechanism has been tested and employed in the laboratory to monitor the evolution of an oil–water encroachment front over time during hydrocarbon recovery Earlier, the streaming potential dynamics were studied near the production well during water-flood of an initially oil-filled [6,7]. Calculations reveal [8] that encroaching water causes changes in the streaming potential at the production well which could be resolved above background electrical noise; water approaching the well could be detected at several 10s to 100s of meters away. The magnitude of the measured potential depends upon the production rate, the coupling between fluid and streaming potentials, and the nature of the front between the displaced oil and the displacing water. Similar results were obtained via simulation for the streaming potential near the injection borehole during drilling [9].
The fact that the induced electric field is orthogonal to the fracture surfaces allows one to tell the fracture direction starting from the electric field testing. If such a field can be strong enough, one can pass to the next step trying to invent a method of measuring the electric field close to the fracture. In this study, we do not address measurements.
In our mathematical approach, we take advantage of the scale difference between fracture disclosure and its length. With fracture disclosure being small, one can reduce fracture–rock interaction onto the central plane between the fracture surfaces. To evaluate the electric field near the fracture, we apply an asymptotic technique assuming that the contrast between the invasion zone and the virgin rock is small. Agreement with laboratory experiments on correlations between the variations of the streaming potential and fluid flow in the core-sample is provided [10]; it follows from our calculations that electric field in the core-sample can be increased in two ways: when the viscosity of the invading fluid decreases or when the electric conductivity of the invading fluid decreases. Thus, the proposed mathematical model explains the laboratory experiments [10] and set theoretical tools for the streaming potential method applied both in the tracking the water front during production operations and the hydro-fracture development.

2. Basic Equations

We treat the fluid-saturated rock as a Bio’s poroelastic medium enjoying quasi-static deformations and incompressible flows which are explained within mathematical models developed in [3,11]. According to this model, the fluid mass balance and the momentum law of the poroelastic medium are given by the equations
( ρ f ϕ ) t + div ( ρ f Q ) = 0 , 0 = div τ + g ρ .
Here,
τ = ( λ div u α * p ) · I + 2 μ E ( u ) , ρ = ϕ ρ f + ( 1 ϕ ) ρ s ,
where ϕ is the porosity, Q is the Darcy velocity, ρ f is the fluid density, ρ s is the solid density, τ is the effective stress tensor, p is the pore pressure, E is the deformation tensor, u is the rock displacement vector, α * is the Biot number, g is the gravitational acceleration vector, μ and λ are the elasticity moduli. We use the notations
( div τ ) i = τ i j x j , div u = u j x j , E i j = 1 2 u i x j + u j x i .
An equilibrium state before fracturing is described by the equations
0 = div τ * + g ρ , τ * = ( λ div u * α * p * ) · I + 2 μ E ( u * ) .
Here, I is the identity tensor, τ * is the equilibrium stress tensor in the rock before fracturing, p = ρ f g h is the equilibrium pore pressure, u * is the equilibrium rock displacement. The equilibrium state is determined by regional stresses as follows. Let us consider a vertical cylinder of the thickness 2 h 0 which is in equilibrium state and which is compressed from above by a poroelastic medium of thickness h, h 0 h . It is proved in [3] that vertical load is transmitted to horizontal compression by the formula
n · τ * n = λ l σ h σ ,
where
σ h = g h ρ , λ l = ν 1 ν + ( 1 2 ν ) p α * σ h ( 1 ν ) , ( τ n ) i = τ i j n j ,
n is the unit external normal vector at the lateral cylinder surface, the point · stands for the scalar product, ν is the Poisson coefficient, λ l is the lateral stress coefficient.
Generally, porosity ϕ depends on pore pressure and rock compressibility, i.e., ϕ = ϕ ( p , J ) , J = tr E ( u ) . By definition of the Biot number α * , we have the equality ϕ / J = α * [12].
Let us introduce the compressibility coefficient S = ϕ / p and the deviation quantities
u = u u * , p = p p , τ = τ τ * .
Then Equation (1) become
α * t div u + S p t + div ( ρ f Q ) = 0 , 0 = div τ ,
where
τ = ( λ div u α * p ) · I + 2 μ E ( u ) .
The fluid flux Q obeys the generalized Darcy law which governs electrolyte flows in a porous medium [13,14,15]:
Q = λ 11 p λ 12 ψ , J = λ 21 p λ 22 ψ , div J = 0 ,
where J is the density of electric current, ψ is the streaming potential, λ i j are electro-kinetic coefficients. It is the double electric layer (DEL) theory which lays behind the first two equations in (4). Particularly, these equations explain the electro-osmosis effect discovered by F.F. Reuss in 1808, Figure 1. The latter equation in (4) is known as the charge conservation law. By the two-scale homogenization theory, it is proved in [16] that λ i j are given by the following representation formulas:
λ 11 = k r η r , λ 22 = σ r , λ 12 = λ 21 = F λ 11 λ 22 ,
where k r is the rock permeability, η r is the dynamical viscosity of the pore fluid, σ r is the effective density of the electric conductivity, F is the dimensionless tortuosity. For sandstones, F 0.001 [8,16].
In [17,18], a mathematical model and a computer code are developed for calculation of the electric conductivity σ r of a saturated rock. The model allows one to determine an optimal Archie-like law
σ r σ f = s * q ϕ ϕ p 1 ϕ p m ,
where σ f is the density of the pore fluid electric conductivity, m is the cementation factor, s * water saturation, q mineralization factor, and ϕ p is the percolation limit.
We assume that a vertical fracture of height 2 r and length 2 R is extended along the horizontal x-axis and its disclosure along another horizontal y-axis does not depend on the vertical variable z, Figure 2a. We also assume that the fracture is symmetrical with respect to the y-axis, and the rock stress-strain state is symmetrical with respect to the x-axis and y-axis. Thus, in what follows we study a two-dimensional problem with all the vectors lying in the horizontal plane ( x , y ) .
Let us write the mass conservation law for fracture fluid. As in the lubrication theory, we have (Batchelor, 1967)
w t + x ( w V ) = q .
Here, w stands for one-half the fracture aperture. By continuity of displacements,
w ( x , t ) = u · e y | y = 0 ,
where e y is the unite vector along the y-axis, q is the fluid leakage, V ( t , x ) is the fluid velocity along the x-axis. The velocity V enjoys the representation
V = ( 2 w ) 2 12 η c p x ,
where η c is the viscosity of the fracture fluid. It is remarked in [3] that one can use the following alternative formula:
V = k c ϕ c η c p x ,
where k c is the fracture permeability and ϕ c is the fracture porosity.
Leakage from the fracture satisfies the continuity condition
q = ( λ 11 p λ 12 ψ ) · e y | y = 0 .
It was proved in [19,20] that the invasion front can be determined by the method of characteristics from the following equation for the position vector ξ = ( ξ 1 , ξ 2 ) at the ( x , y ) -plane:
ϕ d d t ξ = Q ( x , y , t ) | x = ξ 1 , y = ξ 2 .
It means that the infinitesimal displacement d ξ of the front point ( ξ 1 , ξ 2 ) is equal to Q ( ξ 1 , ξ 2 , t ) d t / ϕ , Figure 2b.
Let us comment on the front since it plays a crucial role in the study. Mathematically, the front can be associated with the interface F ( x , y , t ) = 0 , with F ( x , y , t ) being a scalar function. It means that the solution ξ ( t ) of (9) satisfies the equation F ( ξ ( t ) , t ) = 0 . Propagation of the interface F along the normal vector n = F / | F | can be determined from Equation (9) as follows:
n · Q + F t | F | = 0 .
Let us denote by [ f ] ξ the jump of the function f ( x , y ) across the interface F along the normal vector n at the point ξ :
[ f ] ξ = lim δ 0 f ( ξ + δ n ) f ( ξ δ n .
By continuity,
[ p ] ξ = 0 , [ ψ ] ξ = 0 , [ J · n ] ξ = 0 , [ Q · n ] ξ = 0 .
On the other hand, [ λ i j ] ξ 0 . Hence, [ ψ · n ] ξ 0 . The latter condition implies that there is a charge concentration at the interface F. In short, our goal is to determine dynamics of the jump [ ψ · n ] ξ .
Let us formulate boundary-value conditions for the fracture lying in the rectangular
Ω 1 = { ( x , y ) : | x | < a , | y | < H } .
Since the fracture aperture is small compared to its length we reduce fracture-rock interactions to its footprint, Figure 1b,
F c = { ( x , y ) : | x | < R , y = 0 } .
Because of symmetry we study flows and deformations in the half domain Ω 1 + = Ω 1 { y > 0 } only.
We assume that the pressure and stress at the external boundary Ω 1 ( y > 0 ) coincide with the regional values. This is why in terms of deviated quantities we impose the following conditions:
Ω 1 : n · τ n = 0 , s · τ n = 0 , p = 0 , ψ = ψ ,
where s is the unite vector tangential to the boundary Ω 1 , ψ is the prescribed streaming potential. Stress continuity at the rock–fracture interface reduces to the footprint F c as follows
y = 0 , | x | < L : n · τ n = σ p p , s · τ n = 0 , ψ = ψ c ,
where ψ c is the streaming potential inside the fracture. Pressure continuity at the rock–fracture interface implies that the pore pressure and the fracture pressure coincide at y = 0 . Outside the fracture, we impose the following symmetry conditions
y ( u · e x ) = 0 , u · e y = 0 , p y = 0 , ψ y = 0 .
To complete the picture, one should set initial conditions. We shall do it in what follows.

3. Very Long Fracture

Our goal is determine correlations between streaming potential dynamics and invasion front propagation. To perform asymptotic analysis, we pass to a simplified problem under the following hypotheses. We assume that the fracture aperture is the same along the fracture, with the latter being infinite. In this case a = L = and the solution depends on time and the space variable y only, 0 < y < H .
We denote
v ( t , y ) = u · e y , ξ = ξ 2 .
Omitting the prime in the pressure notation p , we arrive at the following equations for the functions v , p , ψ , ξ in the domain 0 < y < H :
0 = y ( λ + 2 μ ) v y α * p ,
α * v y t + S p t = y λ 11 p y + λ 12 ψ y
0 = y λ 21 p y + λ 22 ψ y
ϕ ξ t = λ 11 p y + λ 12 ψ y | y = ξ .
Boundary conditions become
y = 0 : v t = λ 11 p y + λ 12 ψ y , ( λ + 2 μ ) v y + ( 1 α * p ) = σ p , ψ = ψ c ,
y = H : ( λ + 2 μ ) v y α * p = 0 , p = 0 , ψ = ψ .
It is crucial in the present study that the kinetic coefficients λ i j undergo jumps at the invasion front:
λ i j = λ i j , , 0 < y < ξ ( t ) , λ i j + , ξ ( t ) < y < H ,
with λ i j and λ i j + being different constants.
Let us analyse system (12)–(17). Equation (12) implies that the function β 1 ( λ + 2 μ ) v y α * p does not depend on the variable y. Now, it follows from the boundary condition (17) that β 1 = 0 . Hence, one can exclude the function v by the formula
v y = α * λ + 2 μ p .
One can derive from (17) the following boundary condition for pressure at y = 0 :
y = 0 : p = σ p
By Equation (14), the function
β λ 21 p y λ 22 ψ y
does not depend on y. Clearly,
ψ y = β + λ 21 p y λ 22 .
We conclude from (16) and (21) that the fracture disclosure can be determined from the equation
y = 0 : v t = λ 11 λ 12 λ 21 λ 22 p y λ 12 λ 22 β .
Let us integrate equality (20) paying attention to jumps of λ i j :
β ξ = λ 21 ( p | ξ p | 0 ) λ 22 ( ψ | ξ ψ | 0 ) ,
β ( H ξ ) = λ 21 + ( p | H p | ξ ) λ 22 + ( ψ | H ψ | ξ ) .
By excluding ψ | ξ , we find the following representation for β ( t ) :
β = n 1 + n 2 p | ξ n 3 ξ ( t ) + n 4 ( H ξ ( t ) ) , p | ξ p ( t , y ) | y = ξ ( t ) ,
where
n 1 = ψ c ψ + λ 21 λ 22 ( σ p ) , n 2 = λ 21 λ 22 ξ , n 3 = 1 λ 22 , n 4 = 1 λ 22 + .
Here and in what follows, [ f ] ξ stands for jump of the function f ( y ) at the point ξ :
[ f ] ξ = f ( ξ + ) f ( ξ ) lim δ 0 f ( ξ + δ ) f ( ξ δ ) .
Thus, by eliminating displacement and potential, we reduce the problem (12)–(17) to the following mathematical model. We look for unknown functions p ( t , y ) , ξ ( t ) , satisfying the equations:
S + α * 2 λ + 2 μ p t + Q y = 0 , Q Δ λ λ 22 p y + λ 12 λ 22 β ( t ) , Δ λ λ 11 λ 22 λ 12 λ 21 ,
ϕ d ξ d t = Q ( t , y ) | y = ξ ( t ) ,
p | y = 0 = σ p , p | y = H = 0 .
Observe, that solution should obey the no-jump restrictions:
[ p ] ξ = 0 , [ Q ] ξ = 0 .
Let us pass to dimensionless variables
z = y H , p ˜ = p p * , τ = t T , ψ ˜ = ψ ψ * , ζ = ξ H .
Now, the functions p ˜ ( τ , z ) , ζ ( τ ) satisfy the equations
A 1 p ˜ τ = z A 2 λ 11 λ 11 p ˜ z A 3 β ˜ ( τ ) λ 12 λ 22 λ 22 λ 12 , p ˜ | z = 0 = σ p p * , p ˜ | z = 1 = 0 ,
ϕ ζ τ = A 6 β ˜ ( τ ) A 5 p ˜ z | z = ζ ,
on the dimensionless interval 0 < z < 1 . Here, A i are dimensionless parameters:
A 1 = 1 + α * 2 S ( λ + 2 μ ) , A 2 = T H 2 S Δ λ λ 22 = λ 11 ( 1 F 2 ) T H 2 S , A 3 = ψ * T λ 12 H 2 S p * ,
A 4 = p * λ 21 ψ * λ 22 , A 5 = p * T H 2 Δ λ λ 22 = p * T λ 11 ( 1 F 2 ) H 2 , A 6 = ψ * T λ 12 H 2 .
The kinetic coefficients λ i j are stepwise functions:
λ i j = λ i j , 0 < z < ζ ( τ ) , λ i j + , ζ ( τ ) < z < 1 .
We denote
Δ ψ ˜ = ψ ˜ | z = 1 ψ ˜ | z = 0 , Δ p ˜ = p ˜ | z = 1 p ˜ | z = 0 .
The function β ˜ ( τ ) is given by the formula
β ˜ = Δ ψ ˜ A 4 Δ p ˜ + p ˜ ( ζ ( τ ) , τ ) A 4 λ 21 + λ 22 / ( λ 21 λ 22 + ) 1 ζ + ( 1 ζ ) λ 22 / λ 22 + .
By definition, β ˜ is dimensionless electric current density since β ˜ = β H ψ * λ 22 and β J is dimension electric current density.
We formulate initial conditions as follows:
p ˜ | τ = 0 = p ˜ 0 ( z ) , ζ | τ = 0 = 0 .
Choosing p * = σ p , we obtain that p ˜ | z = 0 = 1 and Δ p ˜ = 1 . Observe that we study electric field induced by the fracture closing. Hence, initial fracture aperture is assumed known and it will appear in what follows.
The solution obeys the following no-jump restrictions:
[ p ˜ ] ζ = 0 , A 2 λ 11 λ 11 p ˜ z ζ = β ˜ ( τ ) A 3 λ 12 / λ 12 λ 22 / λ 22 ζ .
Given the functions p ˜ ( z , τ ) and ζ ( τ ) , one can determine dimensionless electric field by the formula
E = ψ ˜ z = β ˜ + ( λ 12 / λ 12 ) A 4 p ˜ z λ 22 / λ 22 .

4. Asymptotic Analysis

Observe, that (27) is not a differential equation in the common sense since it contains not only partial derivatives of p ˜ ( z , τ ) at the running point ( z , τ ) but at the point ( ζ ( τ ) , τ ) also, with ζ ( τ ) being unknown function. It is the main mathematical difficulty of the study. Up to now, no tools are developed to solve the problem numerically. This is why we try to solve it analitically applying an asymptotic analysis. There are publications on differential equations which contain partial derivatives at the running point ( z , τ ) and a point ( z 0 , τ ) , with z 0 being fixed. Such equations are known as loaded differential equations.
We analyse system (27)–(29) under the assumption that jumps of the coefficients λ i j are small in the sense that there is a small ε such that λ i j = k i j ε where the step-wise functions k i j do not depend on ε . We perform an asymptotic analysis assuming that ε 0 . Clearly, both the functions p ˜ ( τ , z ) and ζ ( τ ) depend on ε . But there is a difficulty in addressing the difference p ˜ ε 1 ( τ , z ) p ˜ ε 2 ( τ , z ) starting from Equation (27) for p ˜ ε ( τ , z ) . The problem is that the coefficients λ i j depend on ε also and undergo jumps across the moving unknown line z = ζ ε ( τ ) . To emphasize that the functions λ i j depend on ε , we write λ i j ε ( x , τ ) instead of λ i j ( x , τ ) .
The way to avoid this difficulty is to pass to new variables ( x , τ ) such that all the coefficients λ i j ε ( x , τ ) are defined on the same fixed domain. We define the change of variables ( z , τ ) ( x , τ ) as follows:
x = z ζ ( τ ) 1 ζ ( τ ) , ζ ( τ ) < z < 1 , z ζ ( τ ) ζ ( τ ) , 0 < z < ζ ( τ ) .
In new variables, the function p ˜ ε ( x , τ ) is defined on the fixed interval 1 < x < 1 and the functions p ˜ ε ( x , τ ) , ζ ε ( τ ) satisfy the equations
A 1 p ˜ t ε + x τ p ˜ x ε = x z x A 2 x z a 1 ( x ) p ˜ x ε A 3 β ˜ ( τ ) a 2 ( x ) ,
p ˜ ε | x = 1 = 1 , p ˜ ε | x = 1 = 0 ,
ϕ ζ τ ε = A 6 β ˜ ( τ ) A 5 ζ ε p ˜ x ε | x = 0 ,
where
β ˜ = Δ ψ ˜ A 4 Δ p ˜ ε + p ˜ ε ( 0 , τ ) A 4 [ a 2 ] 0 ζ ε + ( 1 ζ ) a 3 .
The functions a i ( x ) are defined in Appendix A.
The no-jump conditions at x = 0 become
x = 0 : [ p ˜ ε ] 0 = 0 , A 2 a 1 ( x ) x z p ˜ x ε 0 = β ˜ ( τ ) A 3 a 2 ( x ) 0 .
We look for p ˜ ε ( x , τ ) , ζ ε ( τ ) via the expansion series
p ˜ ( x , τ ) = p 0 ( x , τ ) + ε p 1 ( x , τ ) + , ζ ( τ ) = ζ 0 ( τ ) + ε ζ 1 ( τ ) + .
We prove in Appendix A, that the function p 0 ( z , τ ) in the dimensionless physical variables is given by the formula
p 0 = 1 z γ π e A 2 π 2 τ / A 1 sin π z .
We derive from (A1) and (A2) that the function ζ 0 ( τ ) can be determine by the Cauchy problem
ϕ d ζ 0 d τ = A 6 β 0 + A 5 1 + γ e A 2 π 2 τ / A 1 cos ( π ζ 0 ) , ζ 0 | τ = 0 = 0 .
Clearly, ζ 0 ( τ ) increases in time even if Δ ψ ˜ = 0 .
Let us determine approximate dimensionless electric current β a ( τ ) = β 0 + ε β 1 , making use the formulas
ε k i j λ i j = λ i j + λ i j λ i j , ε a 3 1 = 1 λ 22 + λ 22 ,
p 0 | x = 0 = p 0 | z = ζ 0 = 1 ζ 0 γ π e A 2 π 2 τ / A 1 sin π ζ 0 .
We find that
β a ( τ ) = Δ ψ ˜ + A 4 A 4 ( 1 ζ 0 ) 1 λ 12 + λ 12 + ( 1 ζ 0 ) 1 λ 22 + λ 22 Δ ψ ˜
γ A 4 λ 12 + λ 12 λ 22 + λ 22 sin π ζ 0 π e A 2 π 2 τ / A 1 .
In dimension variables the density of the electric current is given by the formula J a = β a ψ * λ 22 / H [Vm 1 s 1 ].
Starting from the definition (31), we calculate dimensionless gradient of potential using the formula
Ψ z ( z , τ ) = β a ( τ ) + ( λ 12 ( z , τ ) / λ 12 ) A 4 p z 0 ( z , τ ) λ 22 ( z , τ ) / λ 22 .
Let us introduce the contrast parameters
a = λ 11 + λ 11 , b = λ 22 + λ 22 .
When Δ ψ ˜ = 0 , it follows from (34) that
Ψ z ( z , τ ) = A 4 ( 1 a b ) [ 1 + sign ( z ζ 0 ( τ ) ) ] / 2 ( 1 ζ 0 ( τ ) ) 1 + ( b 2 1 ) [ 1 + sign ( z ζ 0 ( τ ) ) ] / 2 + F ( z , τ ) e A 2 π 2 τ / A 1 ,
where
F ( z , τ ) = γ A 4 π 1 ( a b + b ) sin ( π ζ 0 ( τ ) ) + cos ( π z ) 1 + ( b 2 1 ) [ 1 + s i g n ( z ζ 0 ( τ ) ) ] / 2 .
We calculated dynamics of Ψ z neglecting the second term in (35) which decreases exponentially in time. Figure 3 depicts the reduced potential gradient Ψ Z ( z , τ ) Ψ z / A 4 across the invasion front z = ζ 0 ( τ ) for the case a = 0.33 and different values of b. Potential gradient jump implies charge concentration at the invasion front z = ζ 0 ( τ ) . The charge concentration is build up as b grows. Observe that the signs of Ψ z are different before and after the front.
Let us calculate Ψ z | τ = 0 within the fracture assuming that Δ ψ ˜ = 0 . We have
Ψ z | z = 0 , τ = 0 = A 4 1 λ 12 + λ 12 + γ A 4
Due to the formula, λ 12 = λ 21 = F λ 11 λ 22 , we derive that
Ψ z | z = 0 , τ = 0 = b a p * F ψ * λ 11 + λ 22 + ( 1 a b + γ )
Let us estimate magnitude of Ψ z | z = 0 , τ = 0 for the typical values of data given in Nomenclature setting γ = 0 for simplicity. In this case, b 2 = 10 4 / 6 . Thus, the dimension initial value of the electric field within the fracture is equal to
ψ * H Ψ z | z = 0 , τ = 0 = 4.53 · 10 4 V / m .
It follows from (36), that electric field in the fracture can be increased in two ways: when the viscosity of the invasion fluid decreases or when the electric conductivity of the invasion fluid decreases. This conclusion agrees with the recent laboratory experiments [10].

5. Attenuation Time

The induced electric field exists until the fracture aperture disappears. To estimate the closure time, we introduce one more dimensionless parameter A 7 = L / H , where L is the one-half the initial fracture aperture. Let v ˜ ( τ ) = v / L stand for dimensionless fracture aperture at the dimensionless moment τ . Due to (22), v ˜ ( τ ) can be determined from the equation
x = 1 : A 7 v ˜ τ = A 5 p ˜ x x z A 6 β ˜ , v ˜ | τ = 0 = v ˜ 0 .
We remind that L v ˜ 0 = d 0 [ s m ] is the dimension initial fracture disclosure. A closure time τ c satisfies the equality
A 7 v ˜ 0 = 0 τ c A 5 p ˜ x x z | x = 1 A 6 β ˜ d τ .
Hence, approximately
A 7 v ˜ 0 = 0 τ c A 5 p x 0 x z 0 | x = 1 A 6 β 0 d τ .
With the initial data given by (A2), the above equation is equivalent to
τ c A 5 + A 4 A 6 Δ ψ ˜ A 6 A 7 v ˜ 0 = γ A 5 A 7 · 1 e A 2 π 2 τ c / A 1 A 2 π 2 / A 1 .
Assuming that Δ ψ ˜ = 0 and using the formula e x 1 x , we derive that
τ c = A 7 v ˜ 0 A 5 A 4 + A 6 + γ A 5 .
The dimension closure time is equal to
T τ c = d 0 H η r p * k r · 1 1 + γ ( 1 F 2 ) d 0 H a 2 p * λ 11 + · 1 1 + γ ( 1 F 2 ) , p * = σ p .
Setting γ 1 / 2 and applying the typical conditions given in Nomenclature, we find that the dimension closure time is equal to T τ c 2 [h], provided d 0 = 0.5 [sm]. One more conclusion is that the closure time decreases as the viscosity of fracture fluid decreases also.

6. Discussions and Conclusions

We formulated a mathematical model to evaluate the streaming potential induced by ionic fracture fluid leak-off after shut-in of a water injection well. Such a potential appears due to electro-kinetic effects relevant to electrolyte flows in a porous medium. The contrast between the virgin fluid and the fluid invading from the fracture is proved to be crucial in a build up of net charge at the invasion front. Calculations reveal that increase of the viscosity and resistivity contrasts results in increase of the streaming potential magnitude. To derive basic equations, we used the scale separation between the fracture disclosure and length. We developed an asymptotic series approach starting from the assumption the contrast parameters is small. The study is motivated by the observation that the vector of electric field correlates with the fracture space orientation since it is directed normally to the fracture surfaces.

Author Contributions

Conceptualization, V.S. and M.E.; methodology, V.S. and M.E.; software, V.S.; validation, V.S. and M.E.; formal analysis, V.S.; investigation, V.S. and M.E.; resources, M.E.; data curation, M.E; visualization, V.S.; supervision, M.E.; funding acquisition, M.E.

Funding

The research is partially supported by Grant of Government of Russian Federation No 14.W03.31.0002.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following notations are used in this manuscript:
α * Biot’s number, dimensionless, 0.7
μ shear elasticity modulus, 8.95 [GPa], 1 [GPa] = 10 g cm · s 2
λ bulk elasticity modulus, λ = 0.5 μ , [GPa]
ν Poisson ratio, dimensionless, 0.25
ρ f fluid density, 1000 kg m 3
ρ s solid density, 2500 kg m 3
ϕ porosity, dimensionless, 0.2
ρ total density, ρ = ϕ ρ f + ( 1 ϕ ) ρ s , kg m 3
ggravitational acceleration, 980 cm s 2
Scompressibility, S 1 =  0.0687 [GPa]
Hdomain’s size of the nearby zone, 10  [m]
Lone-half the initial fracture aperture, 0.5 [cm]
Tcharacteristic time, 3600 [s]
ψ * reference value of streaming potential, 0.1 [mV], 1 [ V 2 ] = 299.79 2 g · cm s 2
hdepth, 2500 [m]
σ h rock’s weight, σ h = ρ h g , [ GPa ]
λ l lateral stress coefficient, dimensionless
σ regional stress, σ = λ l σ h , [ GPa ]
p * reference pressure, p * = σ p , [ GPa ]
k r rock permeability, 1 [ mD ] , 1 [ mD ] = 0.987 · 10 11 cm 2
η r dynamic viscosity of invasion fluid, 0.33 [cp], 1 [cp] = 10 2 g cm · s 2
η r + dynamic viscosity of pore fluid, 3 [cp]
σ f specific electric conductivity of invasion fluid, 0.6 S m , 1 S m = 9 · 10 9 1 s
σ f + specific electric conductivity of pore fluid, 10 5 S m
s * water saturation, dimensionless, s * = 0.7 , s * + = 0.1
mcementation factor, dimensionless, m = 2 , m + = 3
Φ p percolation porosity limit, dimensionless, Φ p = 0.008 , Φ p + = 0.003
σ r effective electric conductivity of saturated rock, S m ,
σ r = s * q σ f ( Φ Φ p ) m , q = 2
Ftortuosity, dimensionless, 10 3
λ i j electrokinetic coefficients
λ i j electrokinetic coefficients in invasion zone
λ i j + electrokinetic coefficients in virgin zone
λ 11 = k r / η r
λ 22 = σ r
λ 12 = λ 21 = F λ 11 λ 22
A 1 = 1 + α * 2 S ( λ + 2 μ ) , dimensionless
A 2 = T H 2 S λ 11 λ 22 λ 12 λ 21 λ 22 = λ 11 ( 1 F 2 ) T H 2 S , dimensionless
A 3 = ψ * T λ 12 H 2 S p * , dimensionless
A 4 = p * λ 21 ψ * λ 22 , dimensionless
A 5 = p * T H 2 λ 11 λ 22 λ 12 λ 21 λ 22 = p * T λ 11 ( 1 F 2 ) H 2 , dimensionless
A 6 = ψ * T λ 12 H 2 , dimensionless
A 7 = L H , dimensionless

Appendix A

It follows from (32) that
x z = 1 1 ζ ( τ ) , 0 < x < 1 , 1 ζ ( τ ) , 1 < x < 0 , x τ = ζ τ ( τ ) ( x 1 ) 1 ζ ( τ ) , 0 < x < 1 , ζ τ ( τ ) ( 1 + x ) ζ ( τ ) , 1 < x < 0 .
Let us introduce the functions
a 1 ( x ) = λ 11 + / λ 11 , 0 < x < 1 , 1 , 1 < x < 0 ,
a 2 ( x ) = λ 12 + λ 22 / ( λ 12 λ 22 + ) , 0 < x < 1 , 1 , 1 < x < 0 ,
and denote a 3 = λ 22 / λ 22 + . Observe that the parameter a 3 and the functions a i ( x ) , β ˜ ( τ ) , x ( z , τ ) depend on ε .
Let us write expansion series for a i ( x ) :
a i = a 0 0 + ε a i 1 +
Starting from the definition of the functions a i ( x ) , we find that a i 0 = 1 , a 3 1 = k 22 / λ 22 , and
a 1 1 = k 11 / λ 11 , 0 < x < 1 , 0 , 1 < x < 0 , a 2 1 = k 12 / λ 12 k 22 / λ 22 , 0 < x < 1 , 0 , 1 < x < 0 ,
Similarly, we find that
x z ( x , τ ) = x y 0 ( x , τ ) + ε x z 1 ( x , τ ) + , x τ ( x , τ ) = x τ 0 ( x , τ ) + ε x τ 1 ( x , τ ) + ,
where
x z 0 = ( 1 ζ 0 ) 1 , 0 < x < 1 , 1 / ζ 0 , 1 < x < 0 ,
x τ 0 = ζ τ 0 ( x 1 ) / ( 1 ζ 0 ) , 0 < x < 1 , ζ τ 0 ( 1 + x ) / ζ 0 , 1 < x < 0 ,
x z 1 = ζ 1 / ( 1 ζ 0 ) 2 , 0 < x < 1 , ζ 1 / ( ζ 0 ) 2 , 1 < x < 0 ,
x τ 1 = [ ζ τ 1 ( x 1 ) ( 1 ζ 0 ) + ζ t 0 ( x 1 ) ζ 1 ] / ( 1 ζ 0 ) 2 , 0 < x < 1 , ζ τ 1 ( 1 + x ) ζ 0 + ζ τ 0 ( 1 + x ) ζ 1 / ( ζ 0 ) 2 , 1 < x < 0 ,
Writing the expansion series
β ˜ ( τ ) = β 0 ( τ ) + ε β 1 ( τ ) + ,
one can derive that
β 0 = Δ ψ ˜ + A 4 , β 1 = A 4 [ a 2 1 ] 0 p 0 ( 0 , τ ) β 0 ( 1 ζ 0 ) a 3 1 .
Setting the expansion series for a i , x ( z , τ ) and β ˜ ( τ ) in Equations (27) and (28), we find that the functions p 0 ( x , τ ) , ζ 0 ( τ ) can be determine from the equations
A 1 a 4 0 ( x ) p τ 0 A 2 p x x 0 A 1 a 5 0 ( x ) p x 0 = 0 , x 0 ,
p 0 | x = 1 = 1 , p 0 | x = 1 = 0 ,
[ p 0 ] 0 = 0 , [ x z 0 p x 0 ] 0 = 0 , ϕ ζ τ 0 = A 6 β 0 A 5 p x 0 | x = 0 / ζ 0 ,
where
a 4 = ( 1 ζ 0 ) 2 , 0 < x < 1 , ( ζ 0 ) 2 , 1 < x < 0 , a 5 = ( x 1 ) ( 1 ζ 0 ) 2 τ / 2 , 0 < x < 1 , ( 1 + x ) ( ζ 0 ) 2 τ / 2 , 1 < x < 0 .
Similarly, we find that the functions p 1 ( x , τ ) and ζ 1 ( τ ) can be determine from the equations
A 1 a 4 0 ( x ) p τ 1 A 2 p x x 1 A 1 a 5 0 ( x ) p x 1 = a 6 0 p τ 0 + a 7 0 p x 0 + a 8 0 p x x 0 , x 0 ,
p 0 | x = 1 = 0 , p 0 | x = 1 = 0 , [ p 1 ] 0 = 0 ,
[ A 2 ( x z 1 p x 0 + x z 0 a 1 1 p x 0 + x z 0 p x 1 ) ] 0 = [ A 3 β 0 a 2 1 ] 0 ,
ϕ ζ τ 1 = A 6 β 1 A 5 p x 1 ζ 0 p x 0 ζ 1 ( ζ 0 ) 2 | x = 0 ,
where a 8 = A 2 a 1 1 and
a 6 = 2 A 1 ( 1 ζ 0 ) ζ 1 , 0 < x < 1 , 2 A 1 ζ 0 ζ 1 , 1 < x < 0 ,
a 7 = A 1 ( ( 1 ζ 0 ) ζ 1 ) τ ( x 1 ) , 0 < x < 1 , A 1 ( ζ 0 ζ 1 ) τ ( 1 + x ) , 1 < x < 0 .
Let us return to the physical dimensionless variables ( z , τ ) , substituting ζ ( τ ) in (32) by ζ 0 ( τ ) . In the variables ( z , τ ) , the functions p 0 ( z , τ ) and ζ 0 ( τ ) satisfy the equations
0 < z < 1 : A 1 p τ 0 = A 2 p z z 0 , p 0 | z = 0 = 1 , p 0 | z = 1 = 0 ,
ϕ d ζ 0 d τ = A 6 β 0 A 5 p z 0 | z = ζ 0 , ζ 0 | τ = 0 = 0 .
To evaluate the electric field, we choose the initial data in (29) as follows
p ˜ | τ = 0 = 1 z γ π sin π z ,
where 0 γ 1 . Under such an assumption, we arrive at Formula (33).

References

  1. Ziagos, J.; Phillips, B.R.; Boyd, L.; Jelacic, A.; Stillman, G.; Hass, E. Technology roadmap for strategic development of enhanced geothermal systems. In Proceedings of the 38-th Workshop on Geothermal Reservoir Engineering, Stanford, CA, USA, 11–13 February 2013. SGP-TR-198. [Google Scholar]
  2. Hunsbedt, A.; Kruger, P.; London, A.L. Energy Extraction from a Laboratory Model Fractured Geothermal Reservoir. J. Petrol. Technol. 1978, 30, 712–718. [Google Scholar] [CrossRef]
  3. Shelukhin, V.V.; Baikov, V.A.; Golovin, S.V.; Davletbaev, A.Y.; Starovoitov, V.N. Fractured water injection wells: Pressure transient analysis. Int. J. Solids Struct. 2014, 51, 2116–2122. [Google Scholar] [CrossRef]
  4. Le Calvez, J.; Malpini, R.; Xu, J.; Stokes, J.; Williams, M. Hydrolic fracturing insights from microseismic monitoring. Oilfeld Rev. 2016, 28, 16–33. [Google Scholar]
  5. Ishio, T.; Mizutani, H. Experimental and the theoretical basis of electrokinetic phenomena in rock-water systems and its application in gheophysics. J. Geophys. Res. 1981, 86, 1763–1775. [Google Scholar] [CrossRef]
  6. Jaafar, M.Z.; Jackson, M.D.; Saunders, J.H.; Vinogradov, J.; Pain, C.C. Measurements of Streaming Potential for Downhole Monitoring on Intelligent Wells; SPE-120460-MS; Society of Petroleum Engineers: Houston, TX, USA, 2009. [Google Scholar]
  7. Jackson, M.D.; Vinogradov, J.; Saunders, J.H.; Jaafar, M.Z. Laboratory measurements and numerical modeling of streaming potential for downhole monitoring in intelligent wells. SPE J. 2011, 16, 625–636. [Google Scholar] [CrossRef]
  8. Saunders, J.H.; Jackson, M.D.; Pain, C.C. A new numerical model of electrokinetic potential response during hydrocarbon recovery. Geophys. Res. Lett. 2006, 33, 1–6. [Google Scholar] [CrossRef]
  9. Eltsov, I.N.; Shelukhin, V.V.; Epov, M.I. Evolution of self-potential near the hole during drilling. Dokl. Earth Sci. 2011, 436, 241–245. [Google Scholar] [CrossRef]
  10. Ghommem, M.; Qiu, X.; Aidagulov, G.; Abbad, M. Streaming potential measurements for downhole monitoring of reservoir fluid flows: A laboratory study. J. Petrol. Sci. 2018, 161, 28–49. [Google Scholar] [CrossRef]
  11. Shelukhin, V.V.; Eltsov, I.N. Geodynamics of the zone close to a borehole during drilling. Dokl. Earth Sci. 2012, 443, 537–542. [Google Scholar] [CrossRef]
  12. Biot, M.A. Theory of elasticity and consolidation for a porous anisotropic solid. J. Appl. Phys. 1955, 26, 182–185. [Google Scholar] [CrossRef]
  13. Amirat, Y.; Shelukhin, V.V. Homogenization of electroosmotic flow equations in porous media. J. Math. Anal. Appl. 2008, 342, 1227–1245. [Google Scholar] [CrossRef]
  14. Hartsen, M.W.; Pride, S.R. Electroseismic waves from point sources in layered media. J. Geophys. Res. 1997, 102, 24745–24769. [Google Scholar] [CrossRef]
  15. Wurmstich, B.; Morgan, F.D. Modeling of streaming potential responses caused by oil pumping. Geophysics 1994, 59, 46–56. [Google Scholar] [CrossRef]
  16. Shelukhin, V.; Eltsov, I.; Paranichev, I. The electrokinetic cross-coupling coefficient: Homogenization approach. World J. Mech. 2011, 1, 127–136. [Google Scholar] [CrossRef]
  17. Shelukhin, V.V.; Terentev, S.A. Frequency dispersion of dielectric permittivity and electric conductivity of rocks via two-scale homogenization of the Maxwell equations. Prog. Electromagn. Res. B 2009, 14, 175–202. [Google Scholar] [CrossRef]
  18. Shelukhin, V.V.; Terentev, S.A. Homogenization of Maxwell equations and the Maxwell-Wagner dispersion. Dokl. Earth Sci. 2009, 424, 155–159. [Google Scholar] [CrossRef]
  19. Shelukhin, V.V. Invasion around a horizontal wellbore. Eur. J. Appl. Math. 2008, 19, 41–60. [Google Scholar] [CrossRef]
  20. Shelukhin, V.V.; Eltsov, I.N. Dynamics of near borehole zone while drilling poroelastic layer. Geophys. J. 2012, 34, 265–272. [Google Scholar]
Figure 1. (a) F.F. Reuss experiment (1808) with water in the U-tube plugged with sandstone sample: applied electric field results in water level change of the hight h (effect of electroosmos). (b) Pressure driven flow in a pore space. (c) Due to the double electric layer (DEL) theory, there is an excess of ions of the same signe in the bulk pore fluid. As a result, a flow occurs if an electric field is applied.
Figure 1. (a) F.F. Reuss experiment (1808) with water in the U-tube plugged with sandstone sample: applied electric field results in water level change of the hight h (effect of electroosmos). (b) Pressure driven flow in a pore space. (c) Due to the double electric layer (DEL) theory, there is an excess of ions of the same signe in the bulk pore fluid. As a result, a flow occurs if an electric field is applied.
Fluids 04 00032 g001
Figure 2. (a) Schematic of a hydrofracture with the height 2 r and the length 2 R . (b) Schematic of an invasion front near the fracture when it is associated with its footprint interval ( R , R ) .
Figure 2. (a) Schematic of a hydrofracture with the height 2 r and the length 2 R . (b) Schematic of an invasion front near the fracture when it is associated with its footprint interval ( R , R ) .
Fluids 04 00032 g002
Figure 3. Reduced potential gradient Ψ Z ( z , τ ) Ψ z / A 4 suffers a jump across the invasion front z = ζ 0 ( τ ) . The contrast parameters are a = 0.33 and b = 0.3 , 0.6 and 1.2 in (ac) respectively.
Figure 3. Reduced potential gradient Ψ Z ( z , τ ) Ψ z / A 4 suffers a jump across the invasion front z = ζ 0 ( τ ) . The contrast parameters are a = 0.33 and b = 0.3 , 0.6 and 1.2 in (ac) respectively.
Fluids 04 00032 g003

Share and Cite

MDPI and ACS Style

Shelukhin, V.; Epov, M. Ionic Fracture Fluid Leak-Off. Fluids 2019, 4, 32. https://doi.org/10.3390/fluids4010032

AMA Style

Shelukhin V, Epov M. Ionic Fracture Fluid Leak-Off. Fluids. 2019; 4(1):32. https://doi.org/10.3390/fluids4010032

Chicago/Turabian Style

Shelukhin, Vladimir, and Mikhail Epov. 2019. "Ionic Fracture Fluid Leak-Off" Fluids 4, no. 1: 32. https://doi.org/10.3390/fluids4010032

Article Metrics

Back to TopTop