Scaling Analysis for the Tracer Flow Problem in Self-Similar Permeability Fields

The spatial variations in porous media (aquifers and petroleum reservoirs) occur at all length scales (from the pore to the reservoir scale) and are incorporated into the governing equations for multiphase flow problems on the basis of random fields (geostatistical models). As a consequence, the velocity field is a random function of space. The randomness of the velocity field gives rise to a mixing region between fluids, which can be characterized by a mixing length = (t). Here we focus on the scale-up problem for tracer flows. Under very general conditions, in the limit of small heterogeneity strengths it has been derived by perturbation theories that the scaling behavior of the mixing region is related to the scaling properties of the self-similar (or fractal) geological heterogeneity through the scaling law (t) ∼ tγ , where γ = max{1/2, 1 − β/2}; β is the scaling exponent that controls the relative importance of short vs. large scales in the geology. The goals of this work are the following: (i) The derivation of a new, mathematically rigorous scaling analysis for the tracer flow problem subject to self-similar heterogeneities. This theoretical development relates the large strength to the small strength heterogeneity regime by a simple scaling of solutions. It follows from this analysis that the scaling law derived by perturbation theory is valid for any strength of the underlying geology, thereby extending the current available results. To the best of the knowledge of the authors this is the only rigorous result available in the literature for the large strength heterogeneity regime. (ii) The presentation of a Monte Carlo study of highly resolved simulations, which are in excellent agreement with the predictions of our new theory. This indicates that our Monte Carlo results are accurate and can be applied to other models for stochastic geology.


Introduction.
The dispersive mixing induced by a random velocity field is of considerable interest in a number of practical applications involving porous media flows.Examples include the remediation of contaminated aquifers and the management of enhanced oil recovery processes.In such flows the growth of the region where the macroscopic mixing of fluids occurs (the mixing region) is determined by multilength-scale rock heterogeneity as well as by fluid instabilities.We focus on the mixing induced by heterogeneity alone and therefore consider unit mobility miscible flow (tracer flow), which is governed by a linear transport equation.Tracer flow corresponds physically to flow of a fluid that carries a passive solute through a porous medium.
The variety of geological mechanisms acting on various length scales support the idea that subsurface formations are heterogeneous at all length scales.It has become increasingly apparent that such heterogeneities, particularly in the permeability field, can have a significant impact on large scale flow.Due to the difficulty in complete and certain characterization of these heterogeneities, stochastic representations of subsurface geologic properties (permeability, porosity) have become commonplace.
We consider self-similar (or fractal) models for geologic variability, which were introduced in earlier theoretical works [17,3,19,16].Fractal models are attractive for two reasons.First, they are at least qualitatively compatible with dispersivity data obtained from a large number of tracer floods [7] in which the dispersivity is observed to increase systematically with travel distance, perhaps in a power law relationship.Second, they are simple enough so that analytic information concerning the associated fluid flow can be obtained at least in some limiting cases.Since the numerical solution of stochastic partial differential equations, including convergence of ensemble averages as the number of realizations is increased, is a barely explored area of numerical analysis, even limited analytic information is extremely valuable as a check on the computational procedures employed.
We remark that as geologic data appears to be nonuniversal and slowly varying as a function of length scale, a multifractal, rather than pure fractal, approach has been introduced in [24].The basic assumption of the multifractal approach is that data varies smoothly in log-log variables, having well-defined tangents with nonuniversal (scale dependent) slopes.Such an assumption provides a very general description of variability, in contrast to the self-similarity assumption, and may provide better models for specific reservoirs or aquifers.
A question of major importance in applications concerns the growth of the mixing region.A characterization of the long-time asymptotic growth of the mixing region was obtained in [16,24] within distinct analytic approximations at the level of perturbation theory.These theories, developed in the Eulerian picture, give the time-asymptotic growth rate of the mixing region as a function of the distance-asymptotic exponent of the permeability field.Finite-time transient effects were also studied in [24].These theoretical predictions are restricted to the regime of small fluctuations of either the permeability [16] or the velocity field [24].We note that the aforesaid results were obtained under ergodic conditions (i.e., assuming infinitely wide tracer plumes).
The main result of the present work is the derivation of a new, mathematically rigorous scaling analysis for the tracer flow problem subject to self-similar permeability heterogeneities.It follows from this analysis that the scaling law derived by perturbation theory is valid for any strength of the underlying heterogeneity, thereby extending the current available results.Such a result for moderately large heterogeneity strengths had been obtained earlier through high-resolution Monte Carlo studies (see [12] and the references therein).To the best of our knowledge our scaling analysis is the only rigorous result available in the literature for the large strength heterogeneity regime.Moreover, we present new Monte Carlo studies employing highly resolved simulations which are in excellent agreement with the predictions of our new theory.In this study both numerical and statistical convergence have been considered for large computational regions (to reduce boundary effects, which have not been included in the theory).Our results indicate that our Monte Carlo approach is accurate and can be applied to other models for geologic variability.
We note that the linear tracer flow problem has been investigated from a number of points of view.We refer the reader to [2] for a renormalization group approach, to [23] for a perturbation expansion analysis, and to [1] for a homogenization treatment.See also [7,6,8,9,4,13,20,25,26] and the references cited therein.The modeling of fluid flow in fractal porous media has also received considerable attention recently; see [21].
The rest of this paper is organized as follows.The model problem along with some known results from perturbation theory for this model are discussed in section 2. Then in section 3 a definition for the mixing length is proposed.The new scaling analysis for the tracer flow problem appears in section 4. Some computational studies which indicate that our Monte Carlo calculations are accurate appear in section 5.In section 6 we compare predictions of our scaling analysis against highly resolved simulations.
2. The stochastic flow model.

Governing equations.
We consider a linear transport equation where s is the saturation of the tagged fluid and the random velocity field v is determined by Darcy's law and the incompressibility condition: Here k ρ is a log-normal permeability field, whose properties are discussed below, μ is the fluid viscosity, and p is the pressure.In formulating (2.1), (2.2) we neglected the effects of gravity, compressibility, and pore scale dispersion and set the porosity equal to a constant (it has been removed from the transport equation by rescaling the time variable).It has been suggested [11], on a somewhat speculative basis due to the very modest availability of field data, that the porosity is linearly correlated with the log-permeability.However, porosity variability is much less than that of permeability, and its effect on fluid mixing is customarily assumed to be negligible as compared with that of permeability variability.
In our theoretical development we consider a steady uniform mean flow, which we assume to be in the x-direction, in an unbounded domain and an initial interface separating the fluid carrying the passive tracer (tagged, s = 1) from that which is not (untagged, s = 0) normal to the mean flow: (2.3) s( x, t 0 ) = 1 when x < 0, 0 when x > 0.
For the numerical simulations we consider (2.1) and (2.2) in a rectangular domain Ω = [0, L x ] × [0, L y ] (Figure 2.1), with boundary conditions where n is the outward-pointing normal vector to ∂Ω.The boundary conditions (2.4) simulate a flow predominantly parallel to the x-axis (left-to-right).The tagged fluid is injected uniformly (at a constant rate q) through the left vertical boundary (x = 0) of Ω, which initially is saturated with untagged fluid.No flow conditions are imposed along the horizontal boundaries y = 0 and y = L y .Fluid is produced from a well kept at constant (zero) pressure at the right vertical boundary (x = L x ).

Self-similar geology.
We consider scalar, log-normal permeability fields so that ξ is Gaussian and its distribution is determined by its mean and covariance function.We assume the distribution to be stationary, isotropic, and, to introduce variability over all length scales, self-similar (or fractal) [15,17].Thus the mean is an absolute constant, which we take to be zero, and the covariance is given by a power law [15]: (Angle brackets denote ensemble averaging.)The scaling exponent β controls the degree of multiscale heterogeneity: As it decreases, the heterogeneities concentrated in the larger length scales are emphasized and the field becomes more regular (locally).The parameter ρ in (2.5), which we dub the heterogeneity strength, sets the overall strength of the fluctuations of the random field k ρ .
This scaling law was consistently derived from leading order results of primitive and renormalized perturbation theory.Moreover, perturbation theory [24,25] was used to derive an effective equation for the tracer flow problem: t 0 is the initial time, and δ v denotes the velocity field fluctuation around its mean value.
The effective equation (2.8) is an approximation to more complex alternatives based on nonlocal diffusion operators; see, for example, [26,5].This approximation is valid in the limit of small velocity fluctuations, and provides quantitative agreement with numerical simulations for a range of moderately large permeability heterogeneities, but it should not be sufficient for sufficiently strong heterogeneities.We use it only to motivate our definition of mixing length in the next section.

The mixing length.
In this section, following [12] we derive an expression for the mixing length in terms of the one-dimensional convection-diffusion model problem Equation (3.1) is a simplified form of the effective equation (2.8).The solution s ν of (3.1)-(3.2) is given by the complementary error function Differentiation of (3.4) followed by integration yields Next, if we notice that the last integral in (3.5) involves the diffusive flux (−ν(t)∂s/∂x), which transports "mass" from the region x < vt to the region x > vt, then we can conclude that 4. Scaling analysis.Here we discuss new scaling relations which allow us to relate the behavior of the tracer flow problem for small and large heterogeneity strengths.In particular, we show that effectively the parameter ρ in (2.5), which sets the heterogeneity strength, rescales the time and spatial variables as well as (t).

Permeability scaling.
Consider the log-normal permeability field defined by (2.5).We note that, for any ρ > 0, where the above equality is in the probabilistic sense, that is, an equality of probability law.Thus, the permeability fields (2.5) satisfy the following scaling relation: where . The scaling property (4.1) in turn entails that flow problems with random permeabilities k ρ1 and k ρ2 are related, as we discuss next.

Flow scaling.
Let s ρ , v ρ , and p ρ solve the stochastic flow equations (2.1)-(2.2) with random permeability k ρ .Note that the initial conditions (2.3) and the uniform mean flow, which is assumed to be in the x-direction in this theoretical development, are both invariant under the spatial scaling (4.1).Then, if In particular, if sρ (x, t) = s ρ ( x, t) denotes the mean saturation profile, which is one dimensional under the assumptions above (i.e., isotropy and stationarity of k ρ ; scale-invariant boundary and initial conditions), then it follows that Equation (4.5) allows us to derive a relation between the mixing lengths corresponding to flow problems with random permeabilities k ρ1 and k ρ2 .
Here L(t) = v t gives the mean distance traveled by the heterogeneous fronts over the time period t.Then the scaling relation (4.5) implies that (4.7) BORGES, FURTADO, PEREIRA, AND AMARAL SOUTO Equation (4.7) elucidates the significance of the heterogeneity strength ρ in the macroscopic diffusion process: The mixing length and the time are scaled by the parameters σ 1 and σ 2 (or, more fundamentally, by ρ 1 and ρ 2 ).In particular, this result shows that the growth of the mixing length for any value of the heterogeneity strength can be obtained from the mixing length growth for a single fixed strength.
We point out the importance of (4.7).The scaling law (2.7) was derived by perturbation theory applied to the fluid transport equation (3.1).As such, it is expected to be valid only under the assumption of weak velocity fluctuations.However, highresolution numerical simulations of the tracer flow problem (2.1)-(2.2) have suggested that these scaling relations are valid for velocity fields produced from permeability fields with very large heterogeneity strengths [15,12].The theory just presented provides an explanation for these numerical results.
The expression (4.7) shows that the mixing length growth for any two values of the heterogeneity strength can be related to each other by a simple scaling.Since (4.7) is mathematically rigorous, it can be used to validate Monte Carlo studies of fluid flow.

Computational results.
We employ Monte Carlo simulations in this study.We are then confronted with two possibilities to define the mixing length.One possibility is to define it for each realization, or sample, in the ensemble of stochastically generated permeability fields; the mixing length of each realization is a random variable, and a statistical moment of this variable would define the mixing length for the ensemble.An alternative possibility is to base the definition of mixing length on the dispersion of the ensemble.It is shown in [12] that the mixing lengths computed using each of the two viewpoints are indistinguishable in our case, and we adopt the definition in terms of the dispersion of the ensemble.
For a given heterogeneity strength ρ we compute numerically sρ (x, t), the mean saturation profile, and ρ (t), the corresponding mixing length.Realizations of the Gaussian random field ξ are generated numerically with the convolution method described in [15].In Figure 5.1 we display one realization of the permeability field k ρ .
To compute sρ , (2.1)-(2.2) are solved numerically for an ensemble of realizations of the permeability field k ρ .These numerical solutions are constructed with an accurate simulator, which combines the mixed finite element method [10] for the solution of the velocity-pressure equation with a second order, nonoscillatory central finite difference scheme [18] for the solution of the saturation equation.The computed saturation profiles s ρ = s ρ ( x, t) are then averaged in the direction transverse to the mean flow (the y-direction) and across the ensemble to define sρ = sρ (x, t) (see Figure 5.2), which is used to compute numerically the mixing length (see (4.6)): (5.1) Here L x is the length of the computational domain in the horizontal direction x.We test our Monte Carlo approach using the scaling relations (4.1), (4.3), (4.5), and (4.7).In order to compare the predictions of such relations with Monte Carlo calculations the following issues have to be considered: • The theoretical results assume an unbounded region.Thus, in our simulations we must consider large computational regions to reduce the boundary effects on the fluid mixing process.• Numerical convergence.The computational grid has to be refined with respect to the grid where the permeability field is given (we refer to the latter grid as the geological grid).We point out that the permeability field is piecewise constant on a uniform grid.
• Statistical convergence.The number of realizations of the permeability field has to be increased until statistical convergence is achieved.The realization of the permeability fields k ρ on a finite lattice provides a short distance regularization of the singular statistics (2.6) due to the lattice cutoff.Then one can define the coefficient of variation where φ k is the standard deviation of the k field.We use the coefficient of variation as a dimensionless measure of the heterogeneity of the numerically generated permeability field.Varying ρ changes the coefficient of variation CV k .The numerical simulations discussed in this section were performed for different values of ρ, which are given in Table 5.1 along with the corresponding values for σ and CV k .The permeability scaling exponent used in our computational studies is β = 0.5; according to (2.7), this value of β corresponds to an anomalous growth rate for the mixing region (γ = 0.75 in (2.7)).

Mesh refinement study.
To verify the effect of the level of refinement of the computational grid on the computed mixing length we performed a mesh BORGES, FURTADO, PEREIRA, AND AMARAL SOUTO refinement study.The computational grids used were 512 × 512, 1024 × 1024, and 2048 × 2048; the geological grid was kept fixed: 512 × 512 grid cells for all cases considered.
We used 64 numerically generated realizations of the permeability field.The value ρ = 0.596 was used for the computational study.The results displayed in Figure 5.3 show the mixing length vs. travel distance in log-log scale.We believe that the discrepancy among the mixing length curves observed at earlier times (short travel distances) is due to numerical diffusion.As the grid is refined, the size of the mixing region is reduced at the beginning of such simulations.This study indicates the numerical convergence of our computed solutions.
We further illustrate the numerical convergence of our computed solutions in Figure 5.4.The frames in this figure show, from top to bottom, saturation surface plots for a given fixed time for the computational grids 512 × 512, 1024 × 1024, and 2048 × 2048, respectively.As the mesh is refined, more detailed fronts are captured.However, this added resolution has no significant influence on the size of the mixing region for later times, as seen in the results reported in Figure 5.3.
In the rest of this work we shall consider computational grids which are twice as fine in each direction as the geological grid.

Statistical convergence study.
In order to determine the number of realizations needed for the statistical convergence of the results for the mixing length growth we compare ensembles with 8, 16, 32, and 64 realizations of a permeability field defined on a 512 × 512 geological grid.A computational grid of size 1024 × 1024 and the value ρ = 1.0 have been considered in this study.The results of this investigation appear in Figure 5.5.This figure displays log-log plots of the mixing length as a function of travel distance for ensembles with increasing size (8, 16, 32, and 64 realizations).This study indicates that statistical convergence has been obtained by averaging over 16 permeability field realizations.

Boundary effects.
The theories for the growth of the mixing length developed in [16,24,14,15] and the scaling analysis developed in this work consider unbounded domains.Aiming at a comparison between Monte Carlo and theoretical results, we investigate the possibility of boundary effects on the computed results.
We consider the tracer flow problem in rectangles having aspect ratios (L x : L y ) 1 : 1, 1 : 1/2, 1 : 1/4, and 1 : 1/8.The corresponding computational (geological) grids used in the numerical studies have 1024 × 1024 (512 × 512), 1024 × 512 (512 × 256), 1024 × 256 (512 × 128), and 1024 × 128 (512 × 64) mesh blocks, respectively.For each one of the aspect ratios we performed an ensemble average study with 16 realizations; the value ρ = 1 was considered.The results of this study appear in Figure 5.6, where the mixing length is shown as a function of the travel distance in log-log scale.We interpret these results as follows: • For the regions with aspect ratios 1 : 1/4 and 1 : 1/8 the results exhibit weak finite domain effects, especially at large length scales.• As the size of the computational region increases in the transverse direction (y-direction) the mixing length curves seem to converge, indicating that the influence of the boundary on the computed mixing length becomes negligible.• The mixing length curves corresponding to regions with aspect ratios 1 : 1 and 1 : 1/2 are very close to each other, suggesting that the influence of the boundary on the computed mixing length is negligible.In the remaining studies of this work we consider regions with aspect ratio 1 : 1; we shall compare the theoretical results with Monte Carlo simulations.

Theory and simulations.
In order to compare the theoretical results given by (4.1), (4.3), (4.5), and (4.7) with the results of Monte Carlo simulations we perform a study for four values of the heterogeneity strength ρ.Table 5.1 shows the values used in our simulations.The geological grid used was 512 × 512 and the computational grid was 1024 × 1024 for all studies reported here.We considered 64 realizations within each ensemble.The covariance for the permeability and velocity fields were computed using a statistical estimator described in [22].
6.1.Permeability field scaling.Figure 6.1 shows the computed covariance function of the permeability fields k ρ characterized by the exponent β = 0.5 for four different values of σ (see Table 5.1).The displayed covariance, in log-log variables, is the result of an average over 64 realizations, each having 512 × 512 elements.We observe good agreement between the computed covariance and the theoretical value (2.6).In Figure 6.2 we see these four curves rescaled according to (4.1), indicating that the statistical structure of our numerically generated realizations of the permeability fields is correct.beyond which the covariance function is dominated by sampling size error.For each of the σ values considered we found good agreement between the computed covariance and the theoretical prediction (see [24]) over a distance of approximately 40 geological mesh blocks.The scaling law for v ρ ( x) that follows from (4.3), namely, (6.1) is observed over a distance of approximately 30 geological mesh blocks in the results reported in Figure 6.4.

Saturation scaling.
In Figure 6.5, we plot the mean saturation profiles (s (x, t)) at times t, 2t, 4t, and 8t for σ = 1/8, 1/4, 1/2, and 1, respectively.These profiles are rescaled according to (4.5) in Figure 6.6.The four distinct profiles collapse  to a single curve.This indicates a very good agreement between the Monte Carlo calculations and the scaling analysis.

Mixing scaling.
We now turn to an investigation of the mixing length growth.Figure 6.7 displays the mixing length ( ρ (t)) in terms of the travel distance, in log-log scale, for the four distinct values of ρ given in Table 5.1.The dotted straight lines in Figure 6.7 have slope 0.75, which is the value of the mixing scaling exponent predicted by the perturbation theory of [16].
In Figure 6.8, the mixing length curves of Figure 6.7 were scaled according to (4.7), taking for reference the ensemble average result corresponding to ρ = 0.596.This result shows that our numerical results are in excellent agreement with the theoretical predictions.The discrepancies at short distances possibly result from the nonfractal behavior of the permeability fields at short distances, due to the finite lattice cutoff, and from numerical dissipation effects in the flow computations.

Asymptotic behavior.
Perturbation theory results [16,24] and numerical simulations [15,12] suggest that the large-time behavior of the mixing length is given by a power law ρ (t) ∼ a(ρ)t γ as t → ∞.
In view of (6.3), we have (6.5)R ρ (t) = ρ in the large-time asymptotic regime.Figure 6.9 shows R ρ (t) for distinct values of ρ.These numerical results are, again, in excellent agreement with our analysis.

Final remarks.
To conclude, we remark that the derivation of the scaling relations in section 4 hinges entirely on the self-similar structure of the permeability fields k ρ and the scaling structure of the flow and transport equations (2.1), (2.2).The latter structure is also observed in the flow equations governing multiphase flow in porous media when dissipative effects (e.g., capillarity effects) are neglected.Therefore we expect these scaling relations to also be valid for multiphase flows.Indeed, preliminary results by the authors and collaborators seem to confirm this prediction; we expect to report our results in the near future.

Fig. 5 . 5 .
Fig. 5.5.Log-log plot of mixing length as a function of travel distance for ensembles with 8, 16, 32, and 64 realizations of permeability fields.Statistical convergence has been obtained by averaging over 16 realizations.

Table 5 . 1
Values of the ρ and corresponding values for σ and CV k .Fig.5.3.Log-log plot of mixing length as a function of travel distance for ρ = 0.596.Results for computational grid sizes 512 × 512, 1024 × 1024, and 2048 × 2048.The geological grid is 512 × 512 for all cases.