Lesson 8: Flow and Transport Processes in 2D Heterogeneous Porous Media

Lesson 8: Flow and Transport Processes in 2D Heterogeneous Porous Media sxr133 Mon, 04/15/2013 - 14:56

Overview

This unit introduces general concepts of spatial heterogeneities, principles of physical flow and transport processes in heterogeneous porous media, as well as how to set up simulations for flow and non-reactive transport in a 2D heterogeneous domain in CrunchFlow.

Learning Outcomes

By the end of this lesson, you should be able to:

Understand fundamental principles and geostatistical measures of spatial heterogeneities
Set up flow and transport simulation in 2D heterogeneous porous medium domain in CrunchFlow

Lesson Roadmap

Lesson Roadmap
Recommended Readings	Knudby and Carrera (2005) Koltermann and Gorelick (1996) (optional) CrunchFlow manual, section on Flow and Transport
To Do	Read Li et al. (2014) for the example question. Heidari and Li (2014). Run tracer tests in the sandboxes with the permeability of high contrast.

Questions?

If you have any questions, please post them to our Questions? discussion forum (not e-mail), located in Canvas. The TA and I will check that discussion forum daily to respond. While you are there, feel free to post your own responses if you, too, are able to help out a classmate.

8.1 Flow and Transport in Multiple Dimensions

8.1 Flow and Transport in Multiple Dimensions ksc17 Tue, 04/23/2013 - 10:21

Flow and transport processes, including advection, dispersion and diffusion are described in the lesson Flow and Transport Processes in 1D System. In natural system we often need to consider these processes in multiple dimensions and in heterogeneous systems where the physical and geochemical properties of the subsurface are not evenly distributed. As an example, in Figure 1, we show that the distribution of permeability dictates the spatial distribution of injected tracer from injection wells during a biostimulation experiments at Old Rifle, Colorado (Li et al., 2010).

Perm distribution on the left and Br distribution on the right. See image description in caption.

Figure 1. A. The spatial distribution of permeability in an biostimulation field experiment plot at Old Rifle, Colorado. Circles are injection wells and squares are monitoring well. B. The spatial distribution of injected tracer (bromide). The injected tracer tends to bypass the low permeability zone (modified from (Li et al., 2010))

The general Advection-Dispersion Equation (ADE) in multiple dimensions is as follows:

\frac{\partial C}{\partial t} + \nabla \cdot (- D \nabla C + v C) = 0

Here C is the tracer concentration (mol/m³ pore volume), t is time (s), D is the combined dispersion–diffusion tensor (m²/s), v (m/s) is the flow velocity vector and can be decomposed into v_L and v_T in the directions longitudinal and transverse to the main flow in a 2D system. The dispersion-diffusion tensor D is defined as the sum of the mechanical dispersion coefficient and the effective diffusion coefficient in porous media D^e(m²/s). At any particular location (grid block), their corresponding diffusion / dispersion coefficients D_L (m²/s) and D_T (m²/s) are calculated as follows:

D_{L} = D_{e} + α_{L} v_{L}

D_{T} = D_{e} + α_{T} v_{T}

Here $α_{L}$ and $α_{T}$ are the longitudinal and transverse dispersivity (m). The dispersion coefficients vary spatially due to the non-uniform permeability distribution. Values of $α_{T}$ are typically at least one order of magnitude smaller than $α_{L}$ (Gelhar et al., 1992; Olsson and Grathwohl, 2007).

8.2 Homogeneity vs. Heterogeneity

8.2 Homogeneity vs. Heterogeneity jls164 Wed, 02/26/2014 - 12:01

The terms homogeneity and heterogeneity are used to describe the uniformity and regularity in spatial distribution of geomaterial properties in natural subsurface systems. Homogeneity means spatially uniform-distributed properties. In natural systems, geological media are almost always heterogeneous. That is, their physical and (bio)geochemical properties vary spatially. Spatial heterogeneity can refer to both physical and geochemical properties. Physical properties include mineral grain size, porosity, and permeability / conductivity. Geochemical properties include, for example, mineral types, lithology, mineral surface area, and cation exchange capacity. Figure 2 shows a picture of an outcrop with layers of different types of geomaterials from the Macrodispersion site (MADE) in Columbia, Mississippi (Zheng and Gorelick, 2003).

Image of layered dirt showing poorly sorted mixture of gravel, sand and silt, coarse sand, and clay matrix

Figure 2. Illustration of spatial heterogeneity in the ground surface where a 1- to 2- meter thick, silty-clay layer is underlain by sand, sandy gravel, and gravelly sand layers.

Credit: Zheng and Gorelick, 2003. Use with permission.

Spatial heterogeneities can lead to significantly different chemical reactions and physical transport processes from their homogeneous counterparts. As a result, heterogeneities play an important role in determining processes and applications in subsurface systems. Examples include contaminant reactive transport (Li et al., 2011), oil and gas production (Chen et al., 2014; Hewett, 1986), and environmental bioremediation (Murphy et al., 1997; Song et al., 2014). For example, as shown in Figure 3, the different geomaterials in Figure 2 leads to orders of magnitude in the spatial variation of hydraulic conductivity (Figure 3A), and abnormal solute transport (Figure 3B) (Zheng et al., 2011).

a 3D visualization of the horizontal hydraulic conductivity, see image caption

Figure 3. A. Three-dimensional visualization of the horizontal hydraulic conductivity distribution based on ordinary kriging of the flowmeter data for the MADE site (Zheng et al., 2011). B. observed MADE-2 tritium plume at 328 days after injection

Zheng and Gorelick 2003).

8.3 Quantitative Measures of Spatial Heterogeneity

8.3 Quantitative Measures of Spatial Heterogeneity jls164 Wed, 02/26/2014 - 12:02

Geostatistics has been well developed to quantitatively describe the characteristics of spatial physical heterogeneities, especially conductivity / permeability, in the natural subsurface. Here we briefly introduce a few geostatistical measures. Interested readers are referred to geostatistical books and software tools for more detailed information (Deutsch and Journel, 1992; Remy and Wu., 2009).

Mean of permeability

It measures the average water-conducting capability of porous media. There are different definitions widely used in subsurface hydrology:

Arithmetic mean $κ_{T}$ is defined as:

κ_{a} = \frac{1}{n} \sum_{i = 1}^{n} κ_{i}

(1)

Here n is the total number of zones or subsystems; $κ_{i}$ is the local permeability of subsystem.

Harmonic mean $κ_{h}$ is defined as:

κ_{h} = n {(\sum_{i = 1}^{n} \frac{1}{κ_{i}})}^{- 1}

(2)

Geometric mean $κ_{g}$ is defined as:

κ_{g} = {(\prod_{i = 1}^{n} κ_{i})}^{\frac{1}{n}}

(3)

Effective permeability $κ_{e}$ is derived from the Darcy’s Law:

κ_{e} = \frac{u μ Δ L}{Δ P}

(4)

where u is the average flow velocity (m/s) of the field; $μ$ is the flow viscosity (Pa·s); $Δ L$ is the length (m) along the main flow direction; $Δ P$ is the pressure differential (Pa) along the main flow direction. The effective permeability describes the real water-conducting capability of the porous media. The upper bound of the effective permeability is often the arithmetic mean, and the lower bound is the harmonic mean. They correspond to flow through a perfectly layered systems, parallel and perpendicular to the layering, respectively (Heidari and Li, 2014; Zinn and Harvey, 2003).

Variance of permeability

The extent of deviation from the mean is calculated as follows:

Var (κ) = \frac{1}{n} \sum_{i = 1}^{n} {(κ_{i} - κ_{a})}^{2}

(5)

where $Var (κ)$ is the permeability variance (m⁴). Larger variance indicates larger extent of spatial heterogeneity. Permeability of geological media differs significantly in natural subsurface systems. Conductivity (K) is also commonly used as a measure of conducting capability of porous media. The Macrodispersion Experiment (MADE) site in Mississippi has both small silt and sand regions. It has a centimeter scale ln(K) variance as large as 20, while its variance at the meter-scale measured by borehole flowmeters is approximately 4.0, which is considerably smaller because small-scale variability is averaged out (Feehley et al., 2000; Harvey and Gorelick, 2000; Zinn and Harvey, 2003). The Wilcox aquifer in Texas also has high ln(K) variance of about 10, as does the Livermore site (ln(K) variance >20) (Fogg, 1986; LaBolle and Fogg, 2002). The Bordon site at Borden, Ontario, however, is relatively homogeneous with ln(K) variance around 0.20 (Mackay et al., 1986; Sudicky, 1986). Note that ln(K) mentioned here is the natural logarithm of hydraulic conductivity K (m/s).

Correlation Length

Correlation length is a measure of correlation in spatial variation. Two points that are separated by a distance larger than the correlation length have fluctuations that are relatively independent, or uncorrelated. In contrast, properties of two points that are within correlation length are correlated. For detailed calculation and application of correlation length in subsurface fields, readers are referred to (Bruderer-Weng et al., 2004; Gelhar, 1986; Heidari and Li, 2014; Salamon et al., 2007). Correlation lengths differ significantly in natural subsurface systems, as shown in Table 1. Figure 4 shows the patterns of permeability fields (m²) with different correlation length generated from Gaussian Sequential Simulation Method.

Figure 4. Representative permeability fields (m²) generated using Gaussian Sequential Simulation Method with different correlation lengths.

Credit: Deutsch and Journel, 1992

Table 1. Correlation length of *ln(K*) of hydraulic conductivity (modified from [*Gelhar*, 1986])
Source	Geological Media	Correlation Length (m): Horizontal	Correlation Length (m): Vertical	Overall Length (m): Horizontal	Overall Length (m): Vertical
Socorro, New Mexico	fluvial sand	>3	0.1	14	5
Rio Grande Valley, New Mexico	fluvial sand	7.6		760
Cape Cod, Massachusetts	glacial outwash sand	5	0.26	20	5
East central Illinois	sandstone aquifer	4.5×10⁴	--	5×10⁵	--
Las Vruces, New Mexico	alluvial silty-clay loam soil	0.1	--	6	--

Hydraulic connectivity

A series of hydraulic connectivity have been defined as integrative measures of spatial heterogeneity characteristics. Static connectivity can be calculated by the connective structure of the hydraulic conductivity or geological facies. Dynamic connectivity directly relates to physical flow or transport processes. For interested readers, we refer them to (Knudby and Carrera, 2005; Knudby and Carrera, 2006; Renard and Allard, 2013; Siirila-Woodburn and Maxwell, 2014; Willmann et al., 2008).

8.4 Setting Up 2D Heterogeneous Domain and Transport Process in CrunchFlow

8.4 Setting Up 2D Heterogeneous Domain and Transport Process in CrunchFlow jls164 Wed, 02/26/2014 - 12:02

Setting up 2D heterogeneous domain for flow and transport calculation generally includes three primary steps. First, the 2D domain needs to be defined with the targeted size, dimension, and resolution. The relevant keyword blocks for setting up the domain include DISCRETIZATION, INITIAL_CONDITION, and BOUNDARY_CONDITION. The second step involves the set up for the calculation of flow, which involves the use of various keywords in the FLOW key block, including constant_flow, calculate_flow, read_permeabilityFile, and pressure. For 2D heterogeneous systems we always use “calculate_flow” to calculate flow by specifying a pressure gradient using “pressure” at the two main boundaries. The permeability distribution can be either defined in INITIAL_CONDITIONS with relevant permeability values in specific grid blocks defined under GEOCHEMICAL_CONDITIONS, or by reading in permeability distribution using the keyword read_permeabilityFile. Detailed format of the permeability file is provided in the manual. After that, the transport keywords can be specified in the TRANSPORT keyword block, similar to the ones discussed in the 1D homogeneous setup.

Example 8.1

Click here for CrunchFlow Files for Example 8.1 and other related files for Li et al. 2014

The following gives an example of setting up flow and transport calculations in 2D homogeneous and heterogeneous domains. Here we use the physical set up of 2D cross-sections of the Mixed and One-zone column in (Li et al., 2014). The authors packed 4 columns of 2.5 cm in diameter and 10.0 cm in length with relatively similar amount of magnesite and quartz distributed in different spatial patterns. Here we focus on two columns that represent two extreme cases: the Mixed column with uniform distributed magnesite and quartz, as well as the One-zone column with magnesite all gathered in one cylindrical zone of diameter of 1.0 centimeter. To keep it simple, here we will do the calculation for 2D cross-sections, instead of following the steps in section 3.2.4 in the paper to convert 2D to 3D. So our numbers might be different from the paper. We are also assuming that in the middle one zone magnesite is 100% of the solid phase. The 2D system has a size of 25 mm by 100 mm. A constant differential pressure is imposed at the boundaries in the z (vertical) direction, leading the primary water flow direction in the z direction from the bottom to the top. No flux boundaries are specified in the X direction (Li et al., 2014).

Table 1. Physical properties of the columns*
Columns	Mg zone	Qtz zone	$^{2} a_{L}$ (cm)	$^{4} a_{T}$ (cm)	$κ_{e}$ (×10^-13m²)	$ϕ_{avg}$
Mixed			0.05	0.005	8.26	0.44
One-zone	Width: 1.0 cm $Φ_{M g :} 0.54$	Width: 1.5 cm $Φ_{Q t z :} 0.38$	0.07	0.004	10.74	0.44

*The permeability of the pure sand columns of the same grain size was measured to be $8.7 \times 10^{- 13} m^{2}$ .

Schematic figures of the 2D cross sections of the Mixed and One-zone columns with the left being the homogeneous column (mixed) and right being the One-zone column having withal magnesite in the mid white zone.

Figure 5. Schematic figures of the 2D cross sections of the Mixed and One-zone columns. The left is for the homogeneous column (mixed) and the right being the One-zone column having magnesite in the mid white zone. Modified based on Figure 1 in Li et al. (2014).

Please do the following:

Given the effective permeability of the One-zone column and the pure sand column, calculate the permeability of the magnesite zone based on the arithmetic average of the equation;
Calculate the volume fraction of the magnesite and quartz assuming the two columns have the same total mass and volume fraction of magnesite and quartz;
Following descriptions in Li et al. (2014) Fig. 1, calculate the pressure field needed for each column in order to have a flow velocity of 3.6 m/day;
What is the residence time when the flow velocity is 3.6 m/day?

Run the simulation at the grid block resolution of 1 mm by 1 mm.

Generate the flow field of the two columns. What is flow velocity in the homogeneous column? What is the flow velocity in the magnesite and sand zone of the One-zone column? What is the ratio of the flow velocity between the two zones? Are they consistent with the permeability ratio of the two zones?
If we inject Br at the concentration of 1.2×10^-4 mol/L at the inlet, please generate the bromide breakthrough curves for the two columns. Note that here you will assume that you only one outlet that integrate all the flow and concentration coming out the outlet grid blocks using the equation
$C_{a v g} = \sum_{i = 1}^{n} C_{i} q_{i} / \sum_{i = 1}^{n} q_{i}$
(1)
where C_i and q_i are the concentration and flow rate from the grid block i in the outlet cross-section. Please plot the overall breakthrough curves from the two columns and discuss their differences.

Note: in plotting 2D figures, you can use softwares such as Tecplot or Matlab. Penn State provides free Matlab web access on WebLabs at Penn State. You can access matlab from there. In addition, you can also get video tutorials from LinkedInLearning at Penn State about using matlab.

Lesson 8 - Mixed 2D.1 Video (56:05 minutes)

Click here for a transcript of the Lesson 2 Ex CrunchFlow video.

Li Li: All right, so this is the video for lesson eight, continued. We just finished doing the homogeneous column. 2D, we call that-- mixed 2D. So you should have watched that first, before you start this one, because the other one has a lot of backdoor information that I provided.

And on top of that, we made the mistake of missing that keyword, zone, in specified permeability. So hopefully, that helps you remember. This is very important.

And then, you have-- I guess, that's the value of making mistakes from my side. And then, you see how things are fixed, and all that. So here, we are going to-- in this one, I'm going to talk about how to do the one zone case. And again, we'll be specifying all the different keyword blocks starting from specify the domain, and then flow, and then transport, and then the reaction-- and then, the initial condition, boundary condition, and all that. So it will be followed in order.

And first, let's also start with discretization. And again, we're using 25 by 100 grid blocks-- 2,500 grid blocks. So essentially, you would have the same discretization. For example, just like what we did before for the homogeneous case-- distance_units. And you would have millimeters. And then, you have xzones-- 25 and space in 1.0 millimeter. Y zones-- you have 100, again 1.0 millimeter as a space.

And then, you are going to specify your domain. But now, we can not specify porosity directly, because there's average porosity, but there's also porosity in different zones, as they are different. The middle zone it's magnesite. The offsite zone is quartz.

And they have different porosity. So we are now going to specify fixed porosity. So here, we are going to just do a porosity update, whatever it is.

And then, we do need to specify the minerals, because different minerals are your different parts of the column. So you actually need magnesite. You also need quartz. And then, you will be specifying different mineral in different parts of the column.

So let's do flow first, just like what we did before, because you do need to calculate its flow, for the flow units are all-- everything you need to calculate the flow field. So again, we are doing time units being seconds. And then, you have distance units being meters. And then, you have calculate flow, true.

So then, you will be specifying permeability. So when you specify permeability, when you do default, you do need to really specify your zone. So I'm going to specify for everything to be 8.74E-13, which is a permeability of the cross-zone first. So when I do default, that means a whole domain have this permeability.

Why do I do that? Because then, I don't need to specify separately. Later on, we'll override. I will specify the magnesite side part, and then override that. So this is the permeability of-- you can specify for the whole domain default. You don't specify zones.

But then I write, later on, permeability_x for the magnesite zone. We calculate that in that, before they're doing this simulation. It's for the middle zone, right? So you have one centimeter with magnesite zone in the middle, and then 0.75 choosing the two other side there.

So into my grid blocks, you would have 8 to 18. Fluctuation, be 9 to 18, is-- sorry-- is the middle magnesite zone.

You would also have this. 13 points-- so that's for the magnesite. Or this would be y. So that specifies the magnesite zone permeability.

OK, so by doing that-- so first of all, here, I specify everywhere is have this value. Now, later, I specify in these zones, you would have these values. So this part of a zone will be overwritten by these values. And the original ones that is not specified is going to remain the same. So you still have-- quartz zone have this value, but the middle of magnesite is going to be overwritten by these values.

And then again, we're going to specify permeability_x like the no-flow boundaries. And you would have zone 0-0, 26-26, as these other ghost cells from-- I'm sorry, 0-0, 1-100, 1-1. And they are permeability_x, 0. And then, again, this zone is very important-- 26, 26.

Actually, I might made a mistake in the previous [INAUDIBLE]. It should be-- really it's a ghost cell outside of the grid block one for x, grid block 25 for x, because these are the no-flow boundary. Again, this will be no flow. So if I made a mistake, make sure you know that was a mistake-- no-flow boundary.

All right, that's permeability. And then, we need to specify pressure. Now see, again, we put the value upon one everywhere. This doesn't really matter. The code we essentially use Darcy's Law given the heterogeneous distribution of permeability and then discretization. So whatever initial value is going to be erased, but we don't want to leave it undefined.

Pressure in that will be 3890. We do the calculation in that before we do the simulation. But again, specified zone will be in the 1-25 in the ghost cell just before you have grid block one. One fix, and then you have pressure, zero again. Zone 1-25-- that's the x. And then, y is 101-101 fix.

So that's a flow. So you are supposed to have everything you need to calculate the flow field. Then, let's look at transport.

Transport will be very similar to what we had before. We need all these dispersion diffusion-related distance units, again, with the centimeters. Time units-- seconds. And then, I'm going to do the same fix diffusion as what we did in the previous one-- 6.342E-6. And then, you have cementation exponent-- cementation, dispersivity.

Dispersivity-- again, the first one is for x, and the second one is for y. So according to the number in the table, you have 0.004, 0.07 centimeters for the two dispersivities. So that's for the transport. And then we need to specify [INAUDIBLE] of the column. Now in here, you have two different zones. One is quart zone. The other is magnesite zone. So we actually specify the two minerals. First of all, we need to specify the primary species.

Now I talk about before, this is a tracer test. So we definitely would need bromide. But also [INAUDIBLE] you put different minerals, you need all the elements, building blocks for that mineral. So for magnesite, you need Mg2 plus because it's MgCO3, so you also need CO2 aq. You need carbonate species.

For quartz, you would need the aqueous species SiO2 aq. But also it is in water. So let's just put H plus there. Just in case some [INAUDIBLE], for example, magnesite is written in terms of H plus. And then second species, you will need OH minus. These species don't really matter in the tracer test. But it will be important later on. But since we are going to specify different mineral zones, even where only for physical process, we do need all the chemical building blocks.

So CO2 aq, then you need bicarbonate, carbonate species. I believe that's it. So these are primary species. And then you need to specify your inlet. Let's just copy all the primary species then, because we know we need that. Inlet, just like what we had before. Units, let's still do molar. Just to be consistent, let's still do molar.

Temperature, bromide, OK, inlet, you have bromide concentration of 1.210 e minus 4. Magnesium, we consider this clean water. So let's just put very low concentration there. Here we also put very low concentration in SiO2, we have low concentration. Because this doesn't really matter at this point. And maybe we put a pH of 7.0. Shouldn't really matter because it's not going to calculate that much anyway.

So for the condition, in the column we have quartz zone and we have magnesite zone. So essentially, we will need both of these. And then primary species, this would be initial in the quartz part. So let's see, still do low concentration bromide. But here, notice these are specified in the aqueous phase. Now we are also going to specify the solid phase.

And for quartz zone, we talk about the positive quartz zone is 0.38. So that means the solid phase is 0.62. So you would have quartz is essentially volume fraction of 0.62. And you put magnesite because you don't have magnesite there. You are going to put very low value.

So essentially, that will give you-- the code will calculate, OK, you have this much solids with a total solid phase of 0.62, then you will calculate porosity being 0.38. I'm just making notes here.

So for magnesite zone, we will still having the whole thing accept the solid phase is different. Because for magnesite, you still have initial condition, that much of these species. But now your quartz will be very small number.

And then your magnesite, remember your magnesite porosity is 0.54. So your magnesite mole volume fraction should be 0.46. So that adding up to be 1. You have the different conditions, just to avoid confusion. So that's the three conditions-- inlet, quartz, magnesite.

Let's specify for the initial conditions. At first, let's specify the whole domain. So there is more quartz [INAUDIBLE] let's specify the whole domain being quartz. So you're putting like one to 25, and then 1 to 100, [INAUDIBLE] for quartz. So that's saying, everywhere you have quartz.

But then, I'm going to override that, the middle part which is nine to 18 will be the magnesite. We're taking the magnesite initial conditions. So that's it gives you, specify initial condition. Then we need a boundary condition. Where is boundary condition?

For the boundary condition, it will be very similar, x_begin, y_begin, x would be inlet. And then you have flux. x_end. Let's say you will do quartz. It doesn't matter much. It's inlet that it is important. And y_begin, y_end. So again, it will take whatever it is coming out of the column. It shouldn't really matter. Boundary condition.

And then you have output. Again, you need to specify how frequently you want things to-- how long do you want to run, how frequent you want output, where do you want output. Spatial profiles, specify, let's do it the same as what we had before-- 0.0011, 0.1, 0.15, and maybe then 0.2. So run the same time as the previous one.

And let's do the time series output, maybe also for the middle. But if you want to look at every outlet block, you're welcome to do that. You don't have to have only one.

For example, you can specify 25 different time series for all the alternate block. Let's say I'm just doing another one time series. Breakthrough, let's say I'm looking at 20. And then [INAUDIBLE] 20, 100, 1. So now you have two different groups. You're looking at the 13 square blocks in the outlet cross section and then 20 square blocks in the outlet.

And then you will be output since this is tracer test. So all that matters to us is bromide. And then we have times series interval should be whatever, 4, 10, or whatever works. All right.

So that seems all we have. Just kind of quickly looking through. I usually try to be clear about what I put in and what are keywords and not making it too crowded. So I'm going to separate that because a lot of times it's easier to see. Line up boundary condition, transport, porosity, discretization, minerals.

For minerals, it's important that you realize that for the name of the minerals, you need the capital in the first letter. Otherwise, the code doesn't recognize it. And then flow, and then you define all the calculated flow, permitted value, permitted value, make sure I have all the zones there when I define specific zone for certain numbers. Zone, zone, OK.

All right. Primary species, secondary species, initial condition, quartz, magnesite. That's consistent always. Let me see. Initial condition, let me just make sure. Yeah, it's consistent with [? permitted that I ?] set up.

Condition inlet, you're putting high concentration bromide and then condition quartz, condition magnesite. You're putting specific volume fraction quartz in various fraction of magnesite, for them to be adding up to be 1 with porosity. All right. So it looks like all right. Let's try it.

So it looks like it's taking this in. And we didn't make mistake again this time, which is good. So first timestamp, it spit it out. And I want to point you to the velocity. So the reasons there's velocity and velocity evolved is because I'm putting porosity to update. So the code will try to update the porosity, although there's no reaction going on or anything. But that's fine.

And then you have porosity. So let's look at velocity. You will realize the velocity is different for different zones, as you can see. And that's more or less because we have different-- permitting different zone. So for this 2D, you will need either product. So spatial profile in 2D like [? Techprod ?] or Matlab, whatever software you use.

The x velocity is not very meaningful. So I would say you just plot the y velocity, which is in direction of the flow. And you realize, OK, y is 1.07. They all have the same units. You're going to go back and check with the corresponding [INAUDIBLE] But in any case, this is 1.07. This is 1.68. Because we just have two zones. And you could imagine this would be the one in the middle zones. And magnesite has a little bit higher probability.

In any case, if you divide 1.68 to 1.07, you get about 1.6 or something. And you realize that when you put the input file, the permitted value of the magnesite is about 1.6 times of the quartz zone. So essentially, the flow is distributed according to the permitted values.

And if you [INAUDIBLE] it should come up with, for example, flow velocity distribution somewhat like the figure one that are we showed with different flow velocity and everything in different places. It's still calculating. It seems really slow for some reason. It's to the 6.35 times 10 to minus 2. OK, still have a bit of time to go.

Because we have more species. And even though we don't really care about these species, we just have to put them there. And the code will try to calculate everything, I'm sure. So it's slower than the homogeneous column.

But in any case, the breakthrough, so for the velocity-- and actually for even let's look at total concentration. So you have x and y. And then you have all different species. But all your K is bromide. So you will be looking at, there's a third column. When you try to calculate the overall breakthrough, imagine you have the 2D. And then in each of the grid block in outlet, you have, for example, let me just go back. It's probably easier to explain with a figure.

This really should be x and y, not x and z. I will change that. In any case, this is your x direction. And this is your y direction. So the [INAUDIBLE] will come out flow velocity different when y equal to 100. This will be y equal to 1, y equal to 100. And this is x.

So you will have different flow velocity coming out from each of the outlet grid blocks. You will also have different concentration coming out because of different [INAUDIBLE]. So what I'm asking you to do, for example, I'll specify a little bit, the concentration average, imagine when we take sample, we actually only have one outlet, not like everywhere you have a sampling point. So that one outlet is essentially integrating all the concentration and flow velocity out of it. So when you calculate average breakthrough, you actually use a flow rate weighted average to do the calculation.

So it would need to be like CI coming out from each of these times, VI and cross section I of each out of this. That way, you get flow rate times concentration, which is mass per time for integrating summation [INAUDIBLE] grid block and divided by the overall flow rate. That will give you the average breakthrough.

So what you need to do when you look at these would be you look at bromide. And you will be looking at y equal to 100. In the example, you get each value, the flow velocity. Cross section area is area of the grid block times porosity. Velocity times that cross section area times the concentration is what you get for each grid block. And then you add them up, divided by the overall flow rate. So that's what you are going to do when you do the breakthrough curve.

And you are going to see a bit of difference in the mix column and in the 2D column. I'm not going to tell you at this point. But in any case, this is what you are going to do. Now we're at 0.16. So we have a little bit. I think we specified to 0.2. So you still have a little bit of time to go. So that's total concentration.

Let me just leave it open. One more. The permitted will give you the distribution of the permitted in different parts of the grid block. These have the units of log 10 meters squared. So these are in log units.

Porosity, [INAUDIBLE] you have xy. And then you have porosity somewhere. For the magnesite zone, you have 54% porosity. And 38 is 38%. So these are the output files. Let me see. Rate, we're not really interested in that.

So it has velocity. If you do have [INAUDIBLE] dissolving out, it will be updating the porosity, permeability, and everything, according to the porosity-permeability relation in the code. And that actually will give you different distributional porosity and permeability from the initial values over time.

I think we are done. Everything is out. And all you need to do is plot thing out and look at the difference in the breakthrough curve, think about why they are different. And you would do-- be interesting to think about how they are different, why they are different, and all that. And I think this is the one zone example with heterogeneity. And you can specify all that.

There is another option. Let's see if you have kind of random distribution of permeability. We can specify in another separate file and [INAUDIBLE] code to read permeability for like irregularly shaped heterogeneity distribution. And this, we're not going to cover. It's a little bit too complicated for now.

So I'm going to stop here. And I'm sure you will have fun. It'll take a little bit more work. But it's also very interesting and very useful when you think about heterogeneity of system, how physical processes are different.

And later on, there will be another session. In Unit Three, we'll be taking about 2D reactive. So this is really building the physical process for that unit later on for the reactive one. All right. I'm going to stop here. I think this is a bit too long.

Click here for the solution.

Solution

Effective permeability of the One-zone column:
$κ_{eff-onezone} = 10.74 \times 10^{- 13} m^{2}$
Pure sand column permeability:
$κ_{sand} = 8.74 \times 10^{- 13} m^{2}$
This indicates that the permeability of sand zone in the One-zone column is $8.74 \times 10^{- 13} m^{2}$ . As the sand and magnesite zones are layered parallel to the flow direction, we can apply Equation (4) to solve the magnesite zone permeability:
$κ_{a} = \frac{1}{n} \sum_{i = 1}^{n} κ_{i}$
Here, the volume fraction for sand zone is:
$(D_{column} - D_{magnesite-zone}) / D_{column} = (2.5 - 1.0) cm / 2.5 cm = 0.6$
The volume fraction for magnesite zone is:
$1 - 0.6 = 0.4$
Hence,
$\begin{matrix} κ_{eff-onezone} = 0.6 κ_{sand-zone} + 0.4 κ_{magnesite-zone} \\ κ_{magnesite-zone} = \frac{κ_{eff-onezone - 0.6 κ_{sand-zone}}}{0.4} \\ = \frac{10.74 \times 10^{- 13} - 0.6 \times 8.74 \times 10^{- 13}}{0.4} m^{2} \\ = 13.74 \times 10^{- 13} m^{2} \end{matrix}$
In the One-zone column, the magnesite volume fraction = volume fraction of magnesite zone * $(1 - f_{M q}) = {0.4}^{*} (1 - 0.54) = 0.184$ . In the Mixed column, the average porosity is 0.44 so the solid phase volume fraction is 0.56. The Mixed column should have the same total amount of magnesite volume as in One-zone column. Therefore The volume fraction of magnesite in the Mixed column = 0.184 and the volume fraction of quartz is the Mixed column = 0.56 – 0.184 = 0.376.
We can apply Darcy's Law, Equation (7) to calculate the pressure gradient:
$κ_{e} = \frac{u μ Δ L}{Δ P}$
Here,
$\begin{matrix} u = 3.6 m / d a y = 4.17 \times 10^{- 5} m / s \\ μ = 1.002 \times 10^{- 3} Pa \cdot s \\ Δ L = 10.0 cm = 0.10 m \end{matrix}$
For the Mixed column,
$κ_{e f f - m i x} = 8.26 \times 10^{- 13} m^{2}$
Hence,
$\begin{matrix} Δ P_{mix} & = \frac{u μ Δ L}{κ_{eff-mix}} \\ = \frac{4.17 \times 10^{- 5} \times 1.002 \times 10^{- 3} \times 0.10}{8.26 \times 10^{- 13}} P a \\ = 5.06 \times 10^{3} Pa \end{matrix}$
For the One-zone column,
$κ_{eff-onezone} = 10.74 \times 10^{- 13} m^{2}$
Hence,
$\begin{matrix} Δ P_{one-zone} & = \frac{u μ Δ L}{κ_{eff-onezone}} \\ = \frac{4.17 \times 10^{- 5} \times 1.002 \times 10^{- 3} \times 0.10}{10.74 \times 10^{- 13}} Pa \\ = 3.89 \times 10^{3} Pa \\ τ & = \frac{L ϕ}{u} = \frac{0.1 \times 0.44}{3.6} = 0.0122 day \end{matrix}$
Insert a video for setting up flow and transport calculation in CrunchFlow here.
Figure 6: Flow fields of the Mixed and One-zone columns
(figure from (Li et al., 2014))
In the homogeneous case, the flow velocity is 3.6 m/day for all the locations. For the One-zone column, the flow velocity in the magnesite zone is 4.5 m/day while the flow velocity in the sand zone is 2.9 m/day. The ratio of flow velocity between the magnesite zone and sand zone is 1.6, which is constant with the permeability ratio of the two zones.
For the 2D homogeneous domain, results should be similar to the 1D example in the lesson Flow and Transport Processes in 1D system except there is an additional parameter transverse dispersivity a_T, a quantitative measure of dispersion in the direction transverse to the flow. Note that the a_T value should not change the shapes of breakthrough curve as a function of residence time in the 2D homogeneous case because there is no transport limitation in the direction transverse to the main flow. Readers can analyze the sensitivity of differential pressure, dispersivity, molecular diffusion coefficient, and cementation factor.
Figure 7. Breakthrough curves of the two columns (figure from (Li et al., 2014)).

Credit: Li Li @ Penn State University is licensed under CC BY-NC-SA 4.0

8.5 Homework Assignment

8.5 Homework Assignment jls164 Wed, 02/26/2014 - 12:02

In Heidari and Li (2014), 2D flow experiments and modeling were used to understand non-reactive solute transport in heterogeneous porous media with different spatial patterns. There are three two-dimensional (2D) sandboxes (21.9 cm × 20.6 cm) packed with the same 20% (v/v) fine and 80% (v/v) coarse sands in three patterns that differ in correlation length: Mixed, Four-zone, and One-zone. The Mixed cases contain uniformly distributed fine and coarse grains. The Four-zone and One-zone cases have four and one square fine zones, respectively (Figure 8).

schematic diagrams of spatial distribution of coarse and fine sand grains, see caption for image description

Figure 8: Schematic diagrams of spatial distribution of coarse and fine sand grains for (a) Mixed case, (b) Four-zone case, (c) One-zone case; (d) schematics of the experimental setup. All three patterns have the same 80% volume of coarse sand and 20% volume of fine sand.

Credit: Heidari and Li, 2014.

Read the paper carefully and follow the example in this lesson to generate the 2D flow fields and tracer breakthrough curves for the three sandboxes of high permeability contrast (HC). All parameters and properties of the three flow cells are in Table 1 and Table 2 in the paper. You can compare your modeling output of the flow field and breakthrough curves with Figure 5 (Flow field and local breakthrough curves) and Figure 6 (overall breakthrough curves from 2D model) in the paper. Note that one pore volume is the same as one residence time.

Note:

Apologies for no original CrunchFlow files for this hw. you can follow Example 8.1 and similarly set it up for this homework.

8.6 Summary and Final Tasks

8.6 Summary and Final Tasks jls164 Wed, 02/26/2014 - 12:03

Summary

In this lesson we studied general concepts and geostatistical measures of spatial heterogeneities, as well as how to set them up in 2D domains.

Reminder - Complete all of the Lesson 8 tasks!

You have reached the end of Lesson 8! Please make sure you have completed all of the lesson activities before you begin Lesson 7.

References

References jls164 Wed, 02/26/2014 - 12:03

Boulanger, R. and Duncan, J. (2000) Geotechnical Engineering Photo Album: University of California at Davis, accessed March 15, 2006.

Bruderer-Weng, C., Cowie, P., Bernabé, Y. and Main, I. (2004) Relating flow channelling to tracer dispersion in heterogeneous networks. Adv. Water Resour. 27, 843-855.

Chen, C., Balhoff, M.T. and Mohanty, K.K. (2014) Effect of Reservoir Heterogeneity on Primary Recovery and CO2 Huff'n'Puff Recovery in Shale-Oil Reservoirs. SPE Reservoir Evaluation & Engineering.

Deutsch, C.V. and Journel, A.G. (1992) Geostatistical software library and user’s guide. New York.

Feehley, C.E., Zheng, C. and Molz, F.J. (2000) A dual‐domain mass transfer approach for modeling solute transport in heterogeneous aquifers: Application to the Macrodispersion Experiment (MADE) site. Water Resour Res 36, 2501-2515.

Fogg, G.E. (1986) Groundwater Flow and Sand Body Interconnectedness in a Thick, Multiple‐Aquifer System. Water Resour Res 22, 679-694.

Gelhar, L.W. (1986) Stochastic subsurface hydrology from theory to applications. Water Resour Res 22, 135S-145S.

Gelhar, L.W., Welty, C. and Rehfeldt, K.R. (1992) A critical review of data on field‐scale dispersion in aquifers. Water Resour Res 28, 1955-1974.

Harvey, C. and Gorelick, S.M. (2000) Rate‐limited mass transfer or macrodispersion: Which dominates plume evolution at the macrodispersion experiment (MADE) site? Water Resour Res 36, 637-650.

Heidari, P. and Li, L. (2014) Solute transport in low‐heterogeneity sandboxes: The role of correlation length and permeability variance. Water Resour Res 50, 8240-8264.

Hewett, T. (1986) Fractal distributions of reservoir heterogeneity and their influence on fluid transport, SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers.

Knudby, C. and Carrera, J. (2005) On the relationship between indicators of geostatistical, flow and transport connectivity. Adv Water Resour 28, 405-421.

Knudby, C. and Carrera, J. (2006) On the use of apparent hydraulic diffusivity as an indicator of connectivity. J Hydrol 329, 377-389.

Koltermann, C.E. and Gorelick, S.M. (1996) Heterogeneity in sedimentary deposits: A review of structure-imitating, process-imitating, and descriptive approaches. Water Resour Res 32, 2617-2658.

LaBolle, E.M. and Fogg, G.E. (2002) Role of molecular diffusion in contaminant migration and recovery in an alluvial aquifer system, Dispersion in Heterogeneous Geological Formations. Springer, pp. 155-179.

Li, L., Gawande, N., Kowalsky, M.B., Steefel, C.I. and Hubbard, S.S. (2011) Physicochemical heterogeneity controls on uranium bioreduction rates at the field scale. Environ Sci Technol 45, 9959-9966.

Li, L., Salehikhoo, F., Brantley, S.L. and Heidari, P. (2014) Spatial zonation limits magnesite dissolution in porous media. Geochimica Et Cosmochimica Acta 126, 555-573.

Li, L., Steefel, C., Kowalsky, M., Englert, A. and Hubbard, S. (2010) Effects of physical and geochemical heterogeneities on mineral transformation and biomass accumulation during a biostimulation experiment at Rifle, Colorado. J Contamin Hydrol 112, 45 - 63.

Mackay, D.M., Freyberg, D., Roberts, P. and Cherry, J. (1986) A natural gradient experiment on solute transport in a sand aquifer: 1. Approach and overview of plume movement. Water Resour Res 22, 2017-2029.

Murphy, E.M., Ginn, T.R., Chilakapati, A., Resch, C.T., Phillips, J.L., Wietsma, T.W. and Spadoni, C.M. (1997) The influence of physical heterogeneity on microbial degradation and distribution in porous media. Water Resour Res 33, 1087-1103.

Olsson, Å. and Grathwohl, P. (2007) Transverse dispersion of non-reactive tracers in porous media: a new nonlinear relationship to predict dispersion coefficients. J Contam Hydrol 92, 149-161.

Remy, N., Alexandre Boucher, and Wu., J. (2009) Applied Geostatistics with SGeMS. Cambridge University Press.

Renard, P. and Allard, D. (2013) Connectivity metrics for subsurface flow and transport. Adv Water Resour 51, 168-196.

Salamon, P., Fernàndez‐Garcia, D. and Gómez‐Hernández, J. (2007) Modeling tracer transport at the MADE site: the importance of heterogeneity. Water Resour Res 43.

Siirila-Woodburn, E.R. and Maxwell, R.M. (2014) A heterogeneity model comparison of highly resolved statistically anisotropic aquifers. Adv Water Resour.

Song, X., Hong, E. and Seagren, E.A. (2014) Laboratory-scale in situ bioremediation in heterogeneous porous media: Biokinetics-limited scenario. J Contam Hydrol 158, 78-92.

Sudicky, E.A. (1986) A natural gradient experiment on solute transport in a sand aquifer: Spatial variability of hydraulic conductivity and its role in the dispersion process. Water Resour Res 22, 2069-2082.

Willmann, M., Carrera, J. and Sánchez‐Vila, X. (2008) Transport upscaling in heterogeneous aquifers: What physical parameters control memory functions? Water Resour Res 44.

Zheng, C., Bianchi, M. and Gorelick, S.M. (2011) Lessons Learned from 25 Years of Research at the MADE Site. Ground Water 49, 649-662.

Zheng, C.M. and Gorelick, S.M. (2003) Analysis of solute transport in flow fields influenced by preferential flowpaths at the decimeter scale. Ground Water 41, 142-155.

Zinn, B. and Harvey, C.F. (2003) When good statistical models of aquifer heterogeneity go bad: A comparison of flow, dispersion, and mass transfer in connected and multivariate Gaussian hydraulic conductivity fields. Water Resour Res 39.