Apéndices | Freeze and Cherry Groundwater Book

Apéndices

Traducido por: Cesar Arturo Vera Florez (Colombia)
Revisado por: Marcela Jaramillo (Colombia)

Apéndice I: Elementos de mecánica de fluidos

El análisis del flujo de agua subterránea requiere un conocimiento de los elementos de la mecánica de fluidos. Albertson & Simons (1964) aportan un análisis breve muy útil; Streeter (1962) y Vennard (1961) son textos estándares. Nuestro propósito aquí es simplemente introducir las propiedades básicas del fluido: densidad de masa, densidad de peso, compresibilidad, y viscosidad; y examinar los conceptos de presión de fluido y carga de presión.

Un análisis de los principios de la mecánica de fluidos debe comenzar con un repaso de la mecánica de materiales en general. La Tabla A1.1 proporciona una lista de las propiedades mecánicas básicas de relevancia, junto con sus dimensiones y unidades en el sistema SI métrico. El sistema SI tiene como dimensiones básicas: masa, longitud y tiempo; con unidades SI básicas: kilogramo (kg), metro (m) y segundo (s). Todas las otras propiedades son medidas en unidades derivadas de las básicas. Algunas de estas propiedades son ampliamente usadas y sus dimensiones básicas tan complejas, que nombres especiales en el SI han sido acuñados para sus unidades derivadas. Como se puede apreciar en la Tabla A1.1, fuerza y peso con medidas en Newton (N), presión y tensión en N/m² o pascales (Pa), y trabajo y energía en Joules (J).

Tabla A1.1 Definiciones, dimensiones, y unidades SI para propiedades mecánicas básicas

					Dimensión de la unidad
Propiedad	Símbolo	Definición	Unidad SI	Símbolo SI	Derivado	Básico
Masa			kilogramo	kg		kg
Longitud			metro	m		m
Tiempo			segundo	s		s
Área	A	A = l²				m²
Volumen	V	V = l³				m³
Velocidad	v	v = l/t				m/s
Aceleración	a	a = l/t²				m/s²
Fuerza	F	F = Ma	newton	N	N	kg∙m/s²
Peso	w	W = Mg	newton	N	N	kg∙m/s²
Presión	p	P = F/A	pascal	Pa	N/m²	kg/m∙s²
Trabajo	W	W = Fl	joule	J	N∙m	kg∙m²/s²
Energía		Trabajo hecho	joule	J	N∙m	kg∙m²/s²
Densidad de masa	ρ	ρ = M/V				kg/m³
Densidad de peso	γ	γ = w/V			N/m³	kg/m²∙s²
Tensión	σ, τ	Respuesta interna a una p externa	pascal	Pa	N/m²	kg/m∙s²
Esfuerzo	ϵ	ϵ = ΔV/V				adimensional
Módulo de Young	E	Ley de Hooke			N/m²	kg/m∙s²

La Tabla A1.2 proporciona un análisis según el sistema SI de ciertas propiedades del fluido y términos de aguas subterráneas que se encuentran en este texto. Cada uno está mucho más explicado en el capítulo 2.

Mucha de la tecnología asociada al desarrollo del recurso agua subterránea en Norteamérica está aún vinculada con el sistema de unidades PLS (pies-libras-segundos), FPS (pie-libra-segundo, FPS en inglés).

Tabla A1.2 Definiciones, dimensiones y unidades SI para propiedades del fluido y términos de agua subterránea

					Dimensión de la unidad
Propiedad	Símbolo	Definición	Unidad SI	Símbolo SI	Derivado	Básico
Volumen	V	V = l³	litro (= m³ × 10^-3)	ℓ	ℓ	m³
Caudal	Q	Q = l³/t			/s	m³/s
Presión de fluido	p	p = F/A	pascal	Pa	N/m²	kg/m∙s²
Carga	h					m
Densidad de masa	ρ	ρ = M/V				kg/m³
Viscosidad dinámica	μ	Ley de Newton	centipoise (= N∙s/m² × 10^-6)	cP	cP, N∙s/m²	kg/m∙s
Viscosidad cinemática	v	v = μ/ρ	centistoke (=m²/s × 10^-6)	cSt	cSt	m²/s
Compresibilidad	α, β	α = 1/E			m²/N
Conductividad hidráulica	K	Ley de Darcy			cm/s
Permeabilidad	k	k = Kμ/pg			cm²	m²
Porosidad	n					adimensional
Almacenamiento específico	S_S	S_S = pg(α + nβ)				1/m
Coeficiente de almacenamiento	S	S = S_Sb*				adimensional
Transmisividad	T	T = Kb*				m²/s

*b, espesor del acuífero confinado (ver Sección 2.10).

Tabla A1.3 Factores de conversión del sistema de unidades PLS (pie-libra-segundo) a unidades del SI

	Multiplicar	Por	Para Obtener
Longitud	pie	3.048 × 10^–1	m
	pie	3.048 × 10	cm
	pie	3.048 × 10^–4	km
	milla	1.069 × 10³	m
	milla	1.069	km
Área	pie²	9.290 × 10^–2	m²
	milla²	2.590	km²
	acre	4.407 × 10³	m²
	acre	4.407 × 10^–3	km²
Volumen	pie³	2.832 × 10^–2	m³
	gal U.S.	3.785 × 10^–3	m³
	gal U.K.	4.546 × 10^–3	m³
	pie³	2.832 × 10
	gal U.S.	3.785
	gal U.K.	4.546
Velocidad	pie/s	3.048 × 10^–1	m/s
	pie/s	3.048 × 10	cm/s
	milla/h	4.470 × 10^–1	m/s
	milla/h	1.609	km/h
Aceleración	pie/s²	3.048 × 10^–1	m/s²
Masa	lbm*	4.536 × 10^–1	kg
	slug*	1.459 × 10	kg
	ton	1.016 × 10³	kg
Fuerza y Peso	lbf*	4.448	N
	poundal	1.383 × 10^–1	N
Presión y Tensión	psi	6.895× 10³	Pa o N/m²
	lbf/ft²	4.788 × 10^–1	Pa
	poundal/ft²	1.488	Pa
	atm	1.013 × 10⁵	Pa
	pulg Hg	3.386× 10³	Pa
	mb	1.000 × 10²	Pa
Trabajo y energía	ft-lbf	1.356	J
	ft-poundal	4.214 × 10^–2	J
	Btu	1.055 × 10^–3	J
	caloría	4.187	J
Densidad de Masa	lb/ft³	1.602 × 10	kg/m³
	slug/ft³	5.154 × 10²	kg/m³
Densidad de Peso	lbf/ft³	1.571 × 10²	N/m³
Caudal	pie³/s	2.832 × 10^–2	m³/s
	pie³/s	2.832 × 10	/s
	gal U.S. /min	6.309 × 10^–5	m³/s
	gal U.K. /min	7.576 × 10^–5	m³/s
	U.S. gal/min	6.309 × 10^–2	/s
	U.K. gal/min	7.576 × 10^–2	/s
Conductividad Hidráulica (ver también Tabla 2.3)	ft/s	3.048 × 10^–1	m/s
	U.S. gal/día/pie²	4.720 × 10^–7	m/s
Transmisividad	ft²/s	9.290 × 10^–2	m²/s
	U.S. gal/día/pie	1.438 × 10^–7	m²/s

* Un cuerpo cuya masa es 1 lb masa (lb_m) tiene un peso de 1 lb fuerza (lb_f). 1 lb_f es la fuerza requerida para acelerar un cuerpo de 1 lb_m a una aceleración de g = 32.2 ft/s². Un slug es la unidad de masa que, cuando actúa bajo una fuerza de 1 lb_f adquiere una aceleración de 1 pie/s².

La Tabla A1.3 provee una gama de factores de conversión para convertir unidades PLS a unidades del SI.

La densidad de masa (o simplemente, densidad) ρ de un fluido es definida como masa por unidad de volumen (Tabla A1.1). La densidad de peso (o peso específico, o unidad de peso) γ de un fluido está definida como su peso por unidad de volumen. Los dos parámetros están relacionados por

$\gamma = \rho g$ (A1.1)

Para agua, ρ = 1.0 g/cm³ = 1000 kg/m³; γ = 9.8 × 10³ N/m³. En el Sistema PLS, γ = 62.4 lb_f/ft³.

La gravedad específica G de cualquier material es la relación de su densidad (o peso específico) respecto al agua. Para agua, G = 1.0; para la mayoría de los suelos y rocas G ≈ 2.65.

La viscosidad de un fluido es la propiedad que permite a los fluidos resistir movimientos relativos y deformación de cizalla durante el flujo. Entre más viscoso el fluido, más grande será la deformación de cizalla a cualquier gradiente de velocidad dado. De acuerdo con la ley de Newton de la viscosidad,

$\tau = \mu \frac{dv}{dy}$ (A1.2)

donde τ es la deformación de cizalla, dv/dy el gradiente de velocidad, y μ la viscosidad, o viscosidad dinámica. La viscosidad cinemática v está dada por

$v = \frac{\mu}{\rho}$ (A1.3)

donde ρ es la densidad de fluido.

La compresibilidad de un fluido refleja sus propiedades de tensión-esfuerzo. La tensión es la respuesta interna de un material a una presión externa. Para fluidos, la deformación es medida a través de la presión del fluido. El esfuerzo es una medida de la deformación lineal o volumétrica de un material deformado. Para fluidos, toma la forma de volumen reducido (y densidad incrementada) bajo incremento de la presión del fluido. La compresibilidad del agua β esta ampliamente discutida en la Sección 2.9. Está definida por la Ec. (2.44).

La densidad, viscosidad y compresibilidad del agua son funciones de la temperatura y presión (Dorsey, 1940; Weast, 1972). En general, su variación no es grande, y de acuerdo con el intervalo de presiones y temperaturas que se encuentran en la mayoría de las aplicaciones de agua subterránea, es común considerarlas como constantes. A 15.5°C, ρ = 1.0 g/cm³, μ = 1.124 cP y β = 4.4 × 10^–2 m²/N.

La presión de fluido p en cualquier punto de un cuerpo vertical de agua es la fuerza por unidad de área que actúa en ese punto. Bajo condiciones hidrostáticas, la presión del fluido en un punto refleja el peso de la columna de agua suprayaciendo a una unidad de área transversal alrededor de un punto. Es posible expresar la presión relativa a presión absoluta cero, pero más comúnmente es expresada según la presión atmosférica relativa. En el último caso esto es llamado presión manométrica, en la medida que esta es la lectura de presión que es obtenida en los manómetros calibrados en cero según la presión atmosférica.

La carga de presión ψ en un punto de un fluido es la altura que una columna de agua podría alcanzar en un manómetro puesto en ese punto. En un cuerpo de agua estático ψ equivale a la profundidad del punto en una medida por debajo de la superficie. Si p es expresada como una presión medida, ψ es definida por la relación

$p = \rho g\psi = \gamma\psi$ (A1.4)

En efecto, la carga de presión ψ es una medida de la presión del fluido p.

Las presiones del fluido son también desarrolladas en el agua subterránea a medida que esta fluye a través de formaciones geológicas porosas y suelos. En la Sección 2.2, los elementos de la mecánica de fluidos presentados en este apéndice están aplicados en el desarrollo de la teoría del flujo de agua subterránea.

Apéndice II: Ecuación de flujo para un flujo transitorio a través de un medio saturado deformado

Un desarrollo riguroso de la ecuación de flujo para un flujo transitorio en un medio saturado debe reconocer el hecho de que cambios transitorios en la presión del fluido llevan a la deformación en el esqueleto granular de un medio poroso, y estas deformaciones implican que el medio, así como el agua, está en movimiento. Esta suposición crea la necesidad de dos ajustes a la derivación clásica creada por Jacob (1940) y seguida en la Sección 2.11 de este texto. Primero, como fue reconocido por (Biot, 1955), es necesario proyectar la ley de Darcy en términos de la velocidad relativa del fluido hacia los granos. Y segundo, como fue reconocido por Cooper (1966) es necesario considerar la conservación de la masa para el medio, así como para el fluido, en el volumen de control elemental. Es posible desarrollar la continuidad de las relaciones en una de tres maneras: (1) considerando un volumen elemental deformable en coordenadas deformables, (2) considerando un volumen elemental deformable en unas coordenadas fijas, o (3) considerando un volumen elemental fijo en coordenadas fijas. Siguiendo a (Gambolati & Freeze, 1973) se usará un volumen elemental fijo en unas coordenadas fijas. Esa aproximación requiere el uso de un vector de notación y el material derivativo (total derivativo, sustancial derivativo). Si estos conceptos no son familiares, Aris (1962) y Wills (1958) proporcionan tratamientos introductorios. El desarrollo será presentado aquí para un medio isotrópico homogéneo, con conductividad hidráulica K, porosidad n, y compresibilidad vertical α. La misma aproximación es fácilmente adaptada a medios heterogéneos y anisotrópicos.

En notación de vector, la forma tridimensional de la ley de Darcy [Ec. (2.34)] es

$\vec{v} = -K\nabla h$ (A2.1)

donde $\vec{v} = (v_x, v_y, v_z)$ es la velocidad relativa de fluido a granos y ∇h = (∂h/∂x, ∂h/∂y, ∂h/∂z) es el gradiente hidráulico. Se puede expandir el vector como

$\vec{v} = n(\vec{v}_w - \vec{v}_s)$ (A2.2)

donde $\vec{v}_w$ es la velocidad del fluido y $\vec{v}_s$ la velocidad del medio deformable. La ecuación de estado para el agua [Ec. (2.47)] es

$\rho = \rho_0 e^\beta^p$ (A2.3)

Y aquella para los granos de suelo, que son incomprensibles, es

ρ_s = constante (A2.4)

La ecuación de continuidad para el agua es

$-\nabla \cdot [np\vec{v}_w] = \frac{\partial}{\partial t}[np]$ (A2.5)

y para el suelo es

$-\nabla \cdot [(1 - n)\rho_s\vec{v}_s] = \frac{\partial}{\partial t}[(1 - n)\rho_s]$ (A2.6)

En estas ecuaciones $\nabla \cdot$ es el operador de divergencia:

Expandiendo la Ec. (A2.5), se llega a

$-\rho\nabla \cdot (n\vec{v}_w) - n\vec{v}_w \cdot \nabla\rho = n\frac{\partial\rho}{\partial t} + \rho\frac{\partial n}{\partial t}$ (A2.7)

Si se cancela ρ_s de la Ec. (A2.6) y se reorganiza esa ecuación, se obtiene una expresión para $n\vec{v}_w$ obtenida de la Ec. (A2.2). Luego de dividir por ρ y reorganizar, la Ec. (A2.7) se vuelve

$-\nabla \cdot \vec{v} - (\frac{v}{\rho}) \cdot \nabla\rho = (\frac{n}{\rho})(\frac{\partial\rho}{\partial t} + \vec{v}_s \cdot \nabla\rho) + \nabla \cdot \vec{v}_s$ (A2.8)

Si se usa el material derivado, $D/Dt = \partial/\partial t + \vec{v}_w \cdot \nabla$ , la Ec. (A2.8) puede ser escrita como

$-\nabla \cdot \vec{v} - (\frac{v}{\rho}) \cdot \partial\rho = \frac{n}{\rho}\frac{Dp}{Dt} + \nabla \cdot \vec{v}_s$ (A2.9)

El primer término a mano derecha de la Ec. (A2.9) puede ser relacionado con la compresibilidad del agua β por la relación

$\frac{D\rho}{Dt} = \rho\beta \frac{Dp}{Dt}$ (A2.10)

El material a mano derecha de la Ec. (A2.10) puede ser reemplazado por una derivada parcial solo si la siguiente desigualdad es satisfecha:

$\vec{v}_s \cdot \nabla p \ll \frac{\partial p}{\partial t}$ (A2.11)

Asumiendo luego, a mano izquierda de la Ec. (A2.9), que

$(\frac{\vec{v}}{\rho}) \cdot \nabla\rho \ll \nabla \cdot \vec{v}$ (A2.12)

Luego, substituyendo las Ecua. (A2.1) y (A2.10) en los campos de la Ec. (A2.9)

$\nabla \cdot (K\nabla h) = n\beta\frac{\partial p}{\partial t} + \nabla \cdot \vec{v}_s$ (A2.13)

En un campo tridimensional de tensiones, el vector de velocidad granular $\vec{v}_s = (v_s_x, v_s_y, v_s_z)$ está relacionado con la deformación (o desplazamiento del suelo) vector $\vec{u} = (u_x, u_y, u_z)$ por

$\vec{v}_s = \frac{D\vec{u}}{Dt}$ (A2.14)

En un campo de tensiones unidireccional vertical,

$v_s_x = v_s_y = u_x = u_y = 0$ (A2.15)

Si las condiciones de la Ec. (A2.15) son satisfechas, el término final de la Ec. (A2.13) puede ser expandido (Cooper, 1966; Gambolati & Freeze, 1973) como (Gambolati, 1974)

$\nabla \cdot \vec{v}_s = \frac{\partial v_s_z}{\partial z} = \frac{\partial}{\partial z}(\frac{Du_z}{Dt}) = \frac{D}{Dt}(\frac{\partial u_z}{\partial z}) = \alpha\frac{Dp}{Dt}$ (A2.16)

Donde α es la compresibilidad vertical del medio poroso. El cambio de derivada alrededor de la igualdad central es válido para un vector de posición, pero no en general. El material derivado en la expresión a mano derecha de la Ec. (A2.16) puede ser reemplazado por la derivada parcial si la Ec. (A2.11) es satisfecha. En ese caso, la Ec. (A2.13) se vuelve

$\nabla \cdot (K\nabla h) = n\beta\frac{\partial p}{\partial t} + \alpha\frac{\partial p}{\partial t}$ (A2.17)

Ya que ρ = ρg(h – z), y si K es constante, la Ec. (A2.17) se simplifica a

$\nabla^2h = \frac{\rho g(\alpha +n\beta}{K}\frac{\partial h}{\partial t}$ (A2.18)

O, renombrando a S_s = ρ g(α + nβ) y expandiendo la notación del vector

$\frac{\partial^2h}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} = \frac{S_s}{K}\frac{\partial h}{\partial t}$ (A2.19)

La Ec. (A2.19) es idéntica a la Ec. (2.75) desarrollada por Jacob (1940). El desarrollo más riguroso hace que sea claro que la validez de la ecuación clásica de flujo depende de que se satisfagan las desigualdades de las Ecs. (A2.11) y (A2.12) y las condiciones de tensión de la Ec. (A2.15). Es poco probable que estas condiciones sean siempre satisfechas. Gambolati (1973b) muestra que donde la tasa de consolidación de $\vec{v}_s$ excede la tasa de percolación del fluido $\vec{v}_w$ , así como sucede en capas gruesas de arcilla de baja permeabilidad y alta compresibilidad, las desigualdades pueden no ser satisfechas. Con respecto a las condiciones de tensión, el término $\nabla \cdot \vec{v}_s$ al final de la Ec. (A2.13) es realmente la punta de un iceberg que relaciona el campo de flujo tridimensional con el campo de tensión tridimensional. Biot (1941; 1955) expuso primero las interrelaciones y Verruijut (1969) proporciona una clara derivada. Schiffman et al. (1969) proporcionan una comparación de las aproximaciones clásicas y Gambolati (1974) analiza el intervalo de validez de la ecuación clásica de flujo.

Apéndice III: Ejemplo de una solución analítica a un problema de contorno

Considerando el problema simple de agua subterránea mostrado en la Figura 2.25 (a). La ecuación de flujo para el flujo saturado en estado estacionario en el plano xy es

$\frac{\partial^2 h}{\partial x^2} + \frac{\partial^2 h}{\partial y^2} = 0$ (A3.1)

El resumen matemático de la condición de contorno es

$\frac{\partial h}{\partial y} = 0$ en $y = 0$ y $y = y_L$ (A3.2)

$h = h_0$ en $x = 0$ (A3.3)

$h = h_1$ en $x = x_L$ (A3.4)

Se resolverá para h(x, y) usando la técnica de separación de variables.

Asumiendo que la solución es un producto solución de la forma

$h(x, y) = X(x) \cdot Y(y)$ (A3.5)

La ecuación (A3.1) luego se vuelve

$Y\frac{\partial^2x}{\partial x^2} + X\frac{\partial^2y}{\partial y^2} = 0$ (A3.6)

Dividiendo por los campos XY

$\frac{1}{X}\frac{\partial^2X}{\partial x^2} = \frac{1}{Y}\frac{\partial^2Y}{\partial y^2}$ (A3.7)

El lado izquierdo es independiente de y. Por lo tanto, el lado derecho, a pesar de su apariencia, debe también ser independiente de y, ya que es idénticamente equivalente a la parte izquierda. De la misma manera, el lado derecho es independiente de x, como también lo es el lado izquierdo. Si ambos lados son independientes de x y y, cada lado equivale a una constante. Entonces,

$\frac{1}{X}\frac{\partial^2X}{\partial x^2} = G$ y $\frac{1}{Y}\frac{\partial^2Y}{\partial y^2} = G$ (A3.8)

La constante G puede ser positiva, negativa o cero. Los tres casos llevan a solución de los productos, pero solo el caso G = 0 lleva a una solución que es físicamente significativa para este problema. La ecuación (A3.8) se vuelve

$\frac{1}{X}\frac{\partial^2X}{\partial x^2} = 0$ y $\frac{1}{Y}\frac{\partial^2Y}{\partial y^2} = 0$ (A3.9)

Estas son ecuaciones diferenciales ordinarias cuyas soluciones son muy bien conocidas:

$X = Ax + B$ y $Y = Cy + D$ (A3.10)

El producto de la solución Ec. (A3.5) se vuelve

$h(x, y) = (Ax + B)(Cy +D)$ (A3.11)

Se puede evaluar los coeficientes A, B, C y D invocando las condiciones de contorno. Diferenciando la Ec. (A3.11) con respecto al campo y se tiene que

$\frac{\partial h}{\partial y} = (Ax + B)C$ (A3.12)

Y una apelación a la Ec. (A3.2) implica que C = 0. Desde la ecuación (A3.11) quedamos con

$h(x, y) = (Ax + B)D = Ex + F$ (A3.13)

Apelando a las condiciones de contorno de los campos de las Ecs. (A33) y (A3.4) arrojan que F = h₀ y E = –(h₀ – h₁)/x_L. La solución es, por lo tanto,

$h(x, y) = h_0 - (h_0 - h_1)\frac{x}{x_L}$ (A3.14)

Esta ecuación es idéntica a la Ec. (2.81), presentada en la Sección 2.11sin prueba. Así es inmediatamente claro que la Ec. (A3.14) satisface las condiciones de contorno de las Ecs. (A3.3) y (A3.4). Derivando con respecto a los campos de y arroja cero en satisfacción de la Ec. (A3.2). Derivando dos veces con respecto a x también arroja cero, entonces la solución de la Ec. (A3.14) satisface la ecuación de flujo [Ec. (A3.1)].

Apéndice IV: Ecuación de Debye-Hückel y la tabla de Kielland para coeficientes de actividad iónica

Expresión de Debye-Hückel para actividades iónicas individuales.

Valores del parámetro å del tamaño del ion para iones comunes encontrados en agua natural.

å × 10⁸	Ion
2.5	$\ce{NH^+_4}$
3.0	K⁺, Cl^–, $\ce{NO^-_3}$
3.5	OH^–, HS^–, $\ce{MnO^-_4}$ , F^–
4.0	SO₄^2–, PO₄^3–, HPO₄^2–
4.0–4.5	Na⁺, $\ce{HCO^-_3}$ , $\ce{H2PO^-_4}$ , $\ce{HSO^-_3}$
4.5	CO₃^2–, SO₃^2–
5	Sr²⁺, Ba²⁺, S^2–
6	Ca²⁺, Fe²⁺, Mn²⁺
8	Mg²⁺
9	H⁺, Al³⁺, Fe³⁺

Parámetros A y B a 1 Bar

Temperatura (°C)	A	B (× 10^–8)
0	0.4883	0.3241
5	0.4921	0.3249
10	0.4960	0.3258
15	0.5000	0.3262
20	0.5042	0.3273
25	0.5085	0.3281
30	0.5130	0.3290
35	0.5175	0.3297
40	0.5221	0.3305
50	0.5319	0.3321
60	0.5425	0.3338

Tabla de Kielland de coeficientes de actividad iónica a 25°C ordenados por el tamaño de los iones (basado en la ecuación de Debye-Hückel)

Carga	*tamaño, å**	Iones	l = 0.0005	0.001	0.0025	0.005	0.01	0.025	0.05	0.1
1	2.5	Rb⁺, Cs⁺, Ag⁺, $\ce{NH^+_4}$ , Tl⁺	0.975	0.964	0.945	0.924	0.898	0.85	0.80	0.75
	3	K⁺, Cl^–, Br^–, I^–, CN^–, $\ce{NO^-_3}$ , $\ce{NO^-_2}$ , OH^–, F^–, $\ce{ClO^-_4}$	0.975	0.964	0.945	0.925	0.899	0.85	0.805	0.755
	4	Na⁺, $\ce{IO^-_3}$ , $\ce{HCO^-_3}$ , $\ce{HSO^-_3}$ , $\ce{H2PO^-_4}$ , $\ce{CLO^-_2}$ , $\ce{C2H3O^-_2}$	0.975	0.964	0.947	0.928	0.902	0.86	0.82	0.775
	6	Li⁺, $\ce{C6H3O^-_2}$	0.975	0.965	0.948	0.929	0.907	0.87	0.835	0.80
	9	H⁺	0.975	0.967	0.950	0.933	0.914	0.88	0.86	0.83
2	4.5	Pb²⁺, Hg₂²⁺, SO₄^2–, CrO₄^2–, CO₃^2–, SO₃^2–, C₂O₄^2–, S₂O₃^2–, H citrate^2–	0.903	0.867	0.805	0.742	0.665	0.55	0.455	0.37
	5	Sr²⁺, Ba²⁺, Cd²⁺, Hg²⁺, S^2–, WO₄^2–	0.903	0.868	0.805	0.744	0.67	0.555	0.465	0.38
	6	Ca²⁺, Cu²⁺, Zn²⁺, Sn²⁺, Mn²⁺, Fe²⁺, Ni²⁺, Co²⁺, phthalate^2–	0.905	0.870	0.809	0.749	0.675	0.57	0.485	0.405
	8	Mg²⁺, Be²⁺	0.906	0.872	0.813	0.755	0.69	0.595	0.52	0.45
3	4	PO₄^3–, Fe(CN)₆^3–, Cr(NH₃)₆³⁺	0.796	0.725	0.612	0.505	0.395	0.25	0.16	0.095
3	9	Al³⁺, Fe³⁺, Cr³⁺, Sc³⁺, In³⁺, and rare earths	0.802	0.738	0.632	0.54	0.445	0.325	0.245	0.18

*Nótese que estos tamaños son valores redondeados para el tamaño efectivo en la solución de agua y no son el tamaño de iones simples, deshidratados. Para una discusión más detallada revise el artículo original de los cuales los valores han sido tomados de [J. Kielland, J. Amer. Chem. Soc., 59, 1675 (1973)].

Referencias

BERNER, R. A. 1971. Principles of Chemical Sedimentology. McGraw-Hill, New York, 240 pp.

KLOTZ, I. M. 1950. Chemical Thermodynamics. Prentice-Hall, Eaglewood Cliffs, N.J., 369 pp.

MANOV, G. G., R. G. BATES, W. J HAMER, and S. F. ACREE. 1943. Values of the constants in the Debye-Hückel equation for activity coefficients. J. Amer. Chem. Soc., 65, pp. 1765–1767.

Apéndice V: Función de Error Complementario (erfc)

$\text{erf}$ $(\beta) = \frac{2}{\sqrt{\mu}}\int_{0}^{\beta}e^{-e^2}d \epsilon$

$\text{erf}$ $(-\beta)$ = $\text{-erf}\beta$

$\text{erfc}$ $(-\beta) = 1 - \text{erf}$ $(\beta)$

$\pmb{\beta}$	erf ( $\pmb{\beta}$ )	erfc ( $\pmb{\beta}$ )
0	0	1.0
0.05	0.056372	0.943628
0.1	0.112463	0.887537
0.15	0.167996	0.832004
0.2	0.222703	0.777297
0.25	0.276326	0.723674
0.3	0.328627	0.671373
0.35	0.379382	0.620618
0.4	0.428392	0.571608
0.45	0.475482	0.524518
0.5	0.520500	0.479500
0.55	0.563323	0.436677
0.6	0.603856	0.396144
0.65	0.642029	0.357971
0.7	0.677801	0.322199
0.75	0.711156	0.288844
0.8	0.742101	0.257899
0.85	0.770668	0.229332
0.9	0.796908	0.203092
0.95	0.820891	0.179109
1.0	0.842701	0.157299
1.1	0.880205	0.119795
1.2	0.910314	0.089686
1.3	0.934008	0.065992
1.4	0.952285	0.047715
1.5	0.966105	0.033895
1.6	0.976348	0.023652
1.7	0.983790	0.016210
1.8	0.989091	0.010909
1.9	0.992790	0.007210
2	0.995322	0.004678
2.1	0.997021	0.002979
2.2	0.998137	0.001863
2.3	0.998857	0.001143
2.4	0.999311	0.000689
2.5	0.999593	0.000407
2.6	0.999764	0.000236
2.7	0.999866	0.000134
2.8	0.999925	0.000075
2.9	0.999959	0.000041
3	0.999978	0.000022

Apéndice VI: Desarrollo de la ecuación de diferencias finitas para el flujo en estado estacionario en un medio isotópico homogéneo

La ecuación de diferencias parciales que describe el flujo en estado estacionario en una región de flujo isotrópica, bidimensional y homogénea, es la ecuación de Laplace (Sección 2.11).

$\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial y^2}$ (A6.1)

Para encontrar la ecuación de diferencias finitas para un nodo interior en la cuadrícula nodal utilizado para discretizar la región de flujo, se debe reemplazar la derivada parcial de segundo orden en la Ec. (A6.1) por diferencias. Consideremos primero el primer término de la ecuación. Recuerde que la definición de derivada parcial con respecto a x de una función de dos variables h(x, z) es

$\frac{\partial h}{\partial x} = \substack {lim \\ \Delta x \rightarrow 0} \frac{h(x + \Delta x, z) - h(x, z)}{\Delta x}$ (A6.2)

En una computadora digital es imposible tomar el límite como $\Delta x \rightarrow 0$ , pero es posible aproximar el límite asignando a Δx algún valor pequeño arbitrario; de hecho, se puede hacer diseñando una red nodal con un espaciado de malla de Δx.

Para cualquier valor de z, ejemplo z₀ se puede expandir h(x, z₀) en una expansión de Taylor sobre el punto (x₀, z₀), como sigue:

$h(x, z_0) = h(x_0, z_0) + (x - x_0)\frac{\partial h}{\partial x}(x_0, z_0) + \frac{(x - x_0)^2}{2}\frac{\partial^2h}{\partial x^2}(x_0, z_0) + ...$ (A6.3)

Si tomamos x = x₀ + Δx (esto es conocido como diferencia hacia adelante) y abandonamos todos los términos de orden más grande que la unidad, se puede aproximar ∂h/∂x por

$\frac{\partial h}{\partial x}(x_0, z_0) = \frac{h(x_0 + \Delta x, z_0) - h(x_0, z_0)}{\Delta x}$ (A6.4)

Los términos abandonados de la expansión de Taylor representan el error de truncamiento en la aproximación de diferencias finitas.

Se puede obtener una expresión similar a la Ec. (A6.4) sustituyendo la diferencia hacia atrás, x = x₀ – Δx en la Ec. (A6.3). Esto produce:

$\frac{\partial h}{\partial x}(x_0, z_0) = \frac{h(x_0, z_0) - h(x_0 - \Delta x, z_0)}{\Delta x}$ (A6.5)

Para obtener la aproximación para ∂₂h/∂x₂, escribimos la ecuación diferencial en términos de ∂h/∂x usando una expresión de diferencia hacia adelante:

$\frac{\partial^2h}{\partial x^2}(x_0, z_0) = \frac{\frac{\partial h}{\partial x}(x_0 + \Delta x, z_0) - \frac{\partial h}{\partial x}(x_0, z_0)}{\Delta x}$ (A6.6)

y sustituyendo la expresión de la diferencia hacia atrás de Ec. (A6.5) en la Ec. (A6.6) para obtener:

$\frac{\partial^2h}{\partial x^2}(x_0, z_0) = \frac{h(x_0 + \Delta x, z_0) - 2h(x_0, z_0) + h(x_0 + \Delta x, z_0)}{(\Delta x)^2}$ (A6.7)

De forma similar, se puede desarrollar una expresión de diferencia para ∂₂h/∂x₂ como

$\frac{\partial^2h}{\partial x^2}(x_0, z_0) = \frac{h(x_0, z_0 + \Delta z) - 2h(x_0, z_0) + h(x_0, z_0 - \Delta z)}{(\Delta z)^2}$ (A6.8)

Para una malla cuadrada, Δx = Δz; añadiendo las Ecs. (A6.7) y (A6.8) para formar la ecuación de Laplace se obtiene

$\frac{1}{(\Delta z)^2}[h(x_0 + \Delta x, z_0) + h(x_0 - \Delta x, z_0) + h(x_0, z_0 + \Delta z) + h(x_0, z_0 - \Delta z) + 4h(x_0, z_0) = 0]$ (A6.9)

Si dejamos que (x₀, z₀) sea el punto nodal (i, j), la Ec. (A6.9) se puede reorganizar para producir

$h_{ij} = \frac{1}{4}(h_{i+1,j} + h_{i-1,j} + h_{i,j+1} + h_{i,j-1})$ (A6.10)

Lo cual es idéntico a la Ec. (5.24).

Apéndice VII: Solución analítica de Tóth para el flujo regional de agua subterránea

Tóth (1962, 1963) ha presentado dos soluciones analíticas para el problema de contorno representando un flujo en estado estacionario, en un campo de flujo isotrópico, vertical, bidimensional, saturado, homogéneo limitado al tope por un nivel freático y en los otros tres lados por límites impermeables (la celda sombreada de la Figura 6.1, reproduciendo en un sistema de coordenadas xz en la figura A7.1).

**Figura A7.1** Región de flujo para la solución analítica de Tóth.

Él, primero consideró el caso en el que la configuración del nivel freático es una línea recta de pendiente constante. Para este caso, la región de flujo en la Figura A7.1 es la región ABCEA. Ya que no es posible obtener una solución analítica en una región trapezoidal, Tóth aproximó la región actual de flujo por la región sombreada ABCDA. Él proyectó los valores de carga hidráulica que existen a lo largo del nivel freático AE sobre el límite superior AD de la región de solución. La aproximación es satisfactoria para una pequeña α.

La ecuación de flujo de Laplace es la ecuación:

$\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial z^2} = 0$ (7.1)

Para una región limitada por x = s y z = z₀, las condiciones de contorno son

$\frac{\partial h}{\partial x}(0, z) = \frac{\partial h}{\partial x}(s, z) = 0$ en AB y CD

$\frac{\partial h}{\partial x}(x, 0) = 0$ en BC (7.2)

$h(x, z_0) = z_0 + cx$ en AD

Donde c = tan α.

La solución analítica, obtenida por separación de variables, es

$h(x, z) = z_0 + \frac{cs}{2} - \frac{4cs}{\pi^2}\sum_{m=0}^{\infty} \frac{\cos[(2m + 1)\pi x/s] \cos h[(2m + 1)\pi z/s]}{(2m +1)^2 \cos h[(2m +1)\pi z_0/s]}$ (7.3)

Esta ecuación satisface tanto la ecuación de flujo (A7.1) como la de las condiciones de frontera (A7.2). Cuando se traza y se contornea, conduce a la red equipotencial que se muestra en la Figura A7.1. El flujo es desde el límite de recarga DF a través del campo al límite del caudal AF.

Tóth también consideró el caso donde la configuración del nivel freático es especificada como una curva seno sobreimpuesta sobre la línea AE (para representar una topografía ondulada). La condición de límite final luego se convierte en

$h(x, z_0) = z_0 + cx + \alpha \sin bx$ en AD (A7.4)

Donde $c = \tan \alpha$ , $a = a'/\cos \alpha$ , y $b = b'/\cos \alpha$ , siendo a’ la amplitud de la curva seno y b’ la frecuencia.

La solución analítica para este caso toma la forma de

$h(x, z) = z_0 + \frac{cs}{2} + \frac{\alpha}{sb}(1 - \cos bs) + 2 \displaystyle\sum_{m=1}^{\infty} \left[\frac{ab(1 - \cos bs \cos m\pi)}{b^2 - \frac{m^2\pi^2}{s^2}} + \frac{cs^2}{m^2\pi^2}(\cos m\pi - 1)\right] \times \frac{\cos (m\pi x/s) \cos h (m\pi z/s)}{s \cdot \cos h (m\pi z_0/s)}$ (A7.5)

La red equipotencial descrita por esta función es similar a aquella de la Figura 6.2 (b).

Apéndice VIII: Solución numérica del problema de contorno representado en una infiltración unidimensional sobre el sistema de recarga del agua subterránea

Consideremos un campo de flujo unidimensional, vertical, homogéneo con su límite superior en la superficie del suelo y su límite inferior bajo el nivel freático. Como fue señalado en la Sección 6.4, la ecuación de flujo por tal sistema es

$\frac{\partial}{\partial z}\left[K(\psi)(\frac{\partial\psi}{\partial z}+ 1)\right] = C(\psi)\frac{\partial\psi}{\partial t}$ (A8.1)

y las condiciones de contorno son

$\frac{\partial\psi}{\partial t} = \frac{R}{K(\psi)} - 1$ (A8.2)

en el tope, y

$\frac{\partial\psi}{\partial z} = \frac{Q}{K_0} - 1$ (A8.3)

en la base. La solución está en términos de carga de presión, $\psi(z, t)$ . La posición del nivel freático en cualquier tiempo está dada por el valor de z en el cual $\psi = 0$ . El aporte inicial incluye las condiciones iniciales, $\psi(z, 0)$ , la tasa de precipitación R, la tasa de recarga del agua subterránea Q y las curvas características no saturadas K(ψ) and C(ψ). Para $\psi \geq \psi_a$ , K = K₀ and C = 0, siendo $\psi_a$ la carga de presión de entrada del aire.

El esquema numérico de diferencia finita usada por Rubin & Steinhardt (1963), Liakopoulos (1965b) y Freeze (1969b) es el esquema implícito de Richtmyer (1957). Con este método el plano (z, t) es representado por una malla rectangular de puntos j = 1, 2, …, L a lo largo del eje z, y n = 1, 2, … a lo largo de del eje t. La distancia entre los nódulos verticales es Δz, y entre los pasos de tiempo, Δt. La solución en cualquier punto dado de la malla (j, n) es $\psi_j^n$ . Bajo esta notación, la forma de la Ec. (A8.1) en diferencias finitas para un nodo interior de (j = 2 a j = L – 1), luego de un adecuado reordenamiento, es

$C(\psi_{\text{III}})\left(\frac{\psi_j^n - \psi_j^{n-1}}{\Delta t}\right) = \frac{1}{\Delta z}\left\{\left[K(\psi_{\text{I}})(1 + \frac{1}{2\Delta z})[\psi_{j+1}^{n-1} + \psi_{j+1}^{n} - \psi_{j}^{n-1} - \psi_j^n]\right]\right\}$ (A8.4)

Donde los valores de $\psi_{\text{I}}$ , $\psi_{\text{II}}$ y $\psi_{\text{III}}$ los cuales determinan los valores de K y C a ser aplicados para un nodo particular en un paso de tiempo particular, están determinados por una extrapolación de los pasos de tiempo previos según lo descrito por Rubin & Steinhardt (1963).

Ecuaciones similares de diferencias finitas, las cuales incorporan las condiciones de contorno (A8.2) y (A8.3), pueden ser escritas por los nodos en los límites superiores e inferiores (j = 1, j = L).

Para cada paso de tiempo, las aproximaciones de diferencias finitas constituyen un sistema de L ecuaciones lineales algebraicas en L desconocidas. La forma general es

$-A_j\psi_{j+1}^n + B_j\psi_j^n - C_j\psi_{j-1}^n = D_j$ (A8.5)

Donde los coeficientes A, B, C y D varían con el nodo (j = 1, j = 2 a L – 1, j = L), con el paso del tiempo (n = 1, n = 2, n > 2) y con la saturación $(\psi \geq \psi_a, \psi < \psi_a)$ . Las variables de las cuales dependen A, B, C y D son los valores de ψ desde el paso del tiempo previo, los valores de frontera R y Q, y las relaciones funcionales K(ψ) y C(ψ).

Los valores de $\psi_j^n$ son calculados desde la relación de recurrencia

$\psi_j^n = E_j\psi_{j+1}^n + F_j$ (A8.6)

donde

$E_j = \frac{A_j}{B_j - C_j E_{j-1}}$

$E_j = \frac{D_j + C_j F_{j-1}}{B_j - C_j E_{j-1}}$

Los coeficientes E y F son calculados desde j = 1 hasta j = L y los ψ son calculados de nuevo desde j = L hasta j = 1 usando la Ec. (A8.6). Para j = 1, E₁ = A₁/B₁ y F₁ = D₁/B₁. Complicaciones surgen si el nodo superior j = L es saturado.

Si K(ψ) y C(ψ) son histeréticas, son usualmente generadas en programas computacionales en la forma de tablas de valores representando la humectación principal, el secado principal y las principales curvas de escaneo. El programa debe localizar la curva correcta en la base según si el nodo en cuestión está mojado o seco, y escaneando la historia pasada del nodo para determinar el valor de ψ al cual un cambio desde mojado a seco o viceversa ha ocurrido.

Apéndice IX: Desarrollo de ecuaciones de diferencias finitas de flujo transitorio en un acuífero confinado heterogéneo, anisotrópico, horizontal

La ecuación diferencial parcial que describe el flujo transitorio a lo largo de un medio saturado anisotrópico (Sección 2.11) es

$\frac{\partial}{\partial x} \left(K_x\frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial y} \left(K_y\frac{\partial h}{\partial y}\right) + \frac{\partial}{\partial z} \left(K_z\frac{\partial h}{\partial z}\right) = S_s\frac{\partial h}{\partial t}$ (A9.1)

Para un acuífero confinado horizontal de espesor, b, la forma bidimensional de la Ec. (A9.1) se reduce a

$\frac{\partial}{\partial x} \left(T_x\frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial y} \left(T_y\frac{\partial h}{\partial y}\right) = S_s\frac{\partial h}{\partial t}$ (A9.2)

Donde T_x y T_y son los componentes principales del tensor de transmisividad definido por T = Kb, y S es el coeficiente de almacenamiento definido por S = S_sb. Para encontrar la ecuación en diferencias finitas para un nodo interior en la cuadricula nodal usada para discretizar la región del flujo, se deben reemplazar las derivadas parciales en la Ec. (A9.2) por diferencias. Usando la definición desarrollada en el apéndice VI y la notación ijk de la Sección 8.8 y la Figura 8.26 (c), se puede escribir una expresión en diferencias finitas:

$\frac{\partial}{\partial x} \left(T_x\frac{\partial h}{\partial x}\right) \simeq \frac{1}{\Delta x}\left[\left(T_x\frac{\partial h}{\partial x}\right)_{i+1/2,j}^k - \left(T_x\frac{\partial h}{\partial x}\right)_{i-1/2,j}^k\right]$ (A9.3)

Donde el subíndice $(i + \frac{1}{2}, j)$ significa que la cantidad entre paréntesis es evaluada en el punto medio entre los nodos (i, j) e (i + 1, j), y el subíndice k significa que la cantidad entre paréntesis es evaluada para el paso del tiempo k. Se pueden aproximar aún más los términos a mano derecha de la Ec. (A9.3) por

$\left(T_x\frac{\partial h}{\partial x}\right)_{i+1/2,j}^k = (T_x)_{i+1/2,j}(h_{i+1,j}^k - h_{i,j}^k)\frac{1}{\Delta x}$ (A9.4A)

$\left(T_x\frac{\partial h}{\partial x}\right)_{i-1/2,j}^k = (T_x)_{i-1/2,j}(h_{i,j}^k - h_{i-1,j}^k)\frac{1}{\Delta x}$ (A9.4B)

Si se evalúa $(T_x)_{i+1/2,j}$ y $(T_x)_{i-1/2,j}$ por simples promedios de la forma

$(T_x)_{i+1/2,j} \simeq \frac{(T_x)_{i+1,j} + (T_x)_{i,j}}{2}$ (A9.5)

Luego esta expresión puede ser sustituida en las Ecs. (A9.4), y las Ecs. (A9.4) pueden por su parte ser sustituidas en la Ec. (A9.3) para obtener

$\frac{\partial}{\partial x} \left(T_x\frac{\partial h}{\partial x}\right) \simeq \frac{1}{2\Delta x^2} \left\{\left[(T_x)_{i+1/2,j} + (T_x)_{i,j}\right]h_{i+1,j}^k - \left[(T_x)_{i+1/2,j} + 2(T_x)_{i,j} + (T_x)_{i-1,j}\right]h_{i,j}^k + \left[(T_x)_{i,j} + (T_x)_{i-1,j}\right]h_{i-1,j}^k \right\}$ (A9.6)

Similarmente:

$\frac{\partial}{\partial y} \left(T_y\frac{\partial h}{\partial y}\right) \simeq \frac{1}{2\Delta y^2} \left\{\left[(T_y)_{i+1/2,j} + (T_y)_{i,j}\right]h_{i+1,j}^k - \left[(T_y)_{i+1/2,j} + 2(T_y)_{i,j} + (T_y)_{i-1,j}\right]h_{i,j}^k + \left[(T_y)_{i,j} + (T_y)_{i-1,j}\right]h_{i-1,j}^k \right\}$ (A9.7)

Finalmente, se puede aproximar la parte derecha de la Ec. (A9.2) por

$S\frac{\partial h}{\partial t} \simeq S_{i,j}\left(\frac{h_{i,j}^k - h_{i,j}^{k-1}}{\Delta t}\right)$ (A9.8)

Sustituyendo las Ecs. (A9.6), (A9.7) y (A9.8) por los tres términos en las Ecs. (A9.6) y agrupando los términos llegamos a una ecuación diferencial finita general para un nodo interno en un acuífero anisotrópico heterogéneo:

$Ah_{i,j}^k = Bh_{i,j-1}^k + Ch_{i+1,j}^k + Dh_{i,j+1}^k + Eh_{i-1,j}^k + F$ (A9.9)

donde

$A = \frac{1}{2\Delta x}^2[(T_x)_{i+1,j} + 2(T_x)_{i,j} + (T_x)_{i-1,j}] + \frac{1}{2\Delta y^2}[(T_y)_{i,j+1} + 2(T_y)_{i,j} + (T_y)_{i,j-1}] + \frac{S_{i,j}}{\Delta t}$

$B = \frac{1}{2\Delta y^2}[(T_y)_{i,j} + (T_y)_{i,j-1}]$

$C = \frac{1}{2\Delta x^2}[(T_x)_{i+1,j} + (T_x)_{i,j}]$

$D = \frac{1}{2\Delta y^2}[(T_y)_{i,j+1} + (T_y)_{i,j}]$

$E = \frac{1}{2\Delta x^2}[(T_x)_{i,j} + (T_x)_{i-1,j}]$

$F = \frac{S_{i,j}}{\Delta t} \cdot h_{i,j}^{k-1}$

Si el acuífero es homogéneo e isotrópico, entonces T_x = T_y = T para todas las (i, j) y S_{i, j} = S para todas las (i, j). Bajo estas condiciones, y para una malla cuadrada nodal con Δx = Δy, los coeficientes de la Ec. (A9.9) se vuelven

$A = \frac{4T}{\Delta x^2} + \frac {S_{i,j}}{\Delta t}$

$B = C = D = E = \frac{T}{\Delta x^2}$

$F = \frac{S}{\Delta t} \cdot h_{i,j}^{k-1}$

Si se divide por T/Δx², estos coeficientes son vistos como equivalentes a aquellos desarrollados de una forma menos rigurosa en la Sección 8.8 y presentados en conexión con la Ec. (8.57).

Trescott et al. (1976) han sugerido que hay cierta ventaja al utilizar la media armónica en lugar de la media aritmética en la Ec. (A9.5). Esta aproximación cambia los coeficientes en la ecuación finita diferencial pero no altera los conceptos que subyacen al desarrollo.

La Ec. (A9.9) está escrita en términos de valores de la carga hidráulica en cinco nodos al paso de tiempo k y un nodo al paso de tiempo k – 1. Eso es conocido como una aproximación de diferencia hacia atrás. Remson et al. (1965) indican que existen algunas ventajas computacionales al usar una aproximación diferencial-central, conocida como esquema Crank-Nicholson, que utiliza valores principales en cinco nodos en el paso de tiempo k y cinco nodos en el paso de tiempo k – 1. El proceso implícito de dirección alternativa (PIDA) utilizado por Pinder & Bredehoeft (1968) involucra dos ecuaciones de diferencias finitas, una en el plano xt y otra en el plano yt. Cada una usa valores principales en los tres nodos en el tiempo k y tres nodos en el tiempo k – 1.

Apéndice X: Derivación de la ecuación de advección-dispersión para transporte de solutos en un medio poroso

La ecuación derivada que se muestra aquí es un resumen de la ley de la conservación de la masa. Dicha derivación está basada en aquellas de Ogata (1970) y Bear (1972). Se asumirá que el medio poroso es homogéneo e isotrópico, que el medio es saturado, que el flujo es estacionario y que la ley de Darcy se aplica. Bajo la suposición de Darcy, el flujo es descrito por la velocidad lineal promedio, que lleva la sustancia disuelta por advección. Si este fuera el único mecanismo de transporte operativo, los solutos no reactivos que están siendo transportados por el flujo se moverían como un tapón. En realidad, existe un proceso adicional de mezcla, que se denomina dispersión hidrodinámica (Sección 2.13), que es causada por variaciones microscópicas en la velocidad dentro de cada canal de poro y de un canal a otro. Si deseamos describir el proceso de transporte en una escala macroscópica usando parámetros macroscópicos, y aun así tomar en cuenta el efecto de la mezcla microscópica, es necesario introducir un segundo mecanismo de transporte, en adición a la advección, para contar por la dispersión hidrodinámica.

Para establecer el resumen matemático de la conversión de las masas, será considerado el volumen elemental pequeño de un flujo de soluto hacia adentro y hacia afuera de un medio poroso (Figura A10.1). En unas coordenadas cartesianas el caudal específico v tiene componentes (v_x, v_y, v_z) y la velocidad lineal promedio $\vec{v} = v/n$ tiene componentes $(\vec{v}_x, \vec{v}_y, \vec{v}_z)$ . La tasa de transporte advectivo es equivalente a $\vec{v}$ . La concentración del soluto C es definida como la masa del soluto por unidad de volumen de solución. Por su parte, la masa del soluto por unidad de volumen de un medio poroso es nC. Para un medio homogéneo, la porosidad n es constante, y $\partial(nC)/\partial x = n\partial C/\partial x$ .

**Figura A10.1** Balance de masas en un elemento cubico en el espacio.

La masa del soluto transportada en dirección x por los dos mecanismos de transporte de soluto, puede ser representada como

Transporte por advección = $\vec{v}_xnC dA$ (A10.1)

Transporte por dispersión = $nD_x\frac{\partial C}{\partial x} dA$ (A10.2)

Donde D_x es el coeficiente de dispersión en la dirección x y dA es el área de sección transversal elemental del elemento cúbico. La dispersión del coeficiente D_x está relacionada con la dispersividad α_x y el coeficiente de difusión D* por la Ec. (9.4):

$D_x = \alpha_x\vec{v}_x + D$ (A10.3)

La forma del componente dispersivo incorporado en la Ec. (A10.2) es análoga a la primera ley de Fick.

Si F_x representa el total de masa de soluto por unidad de área transversal transportada en la dirección x por unidad de tiempo, entonces

$F_x = \vec{v}_xnC - nD_x\frac{\partial C}{\partial x}$ (A10.4)

El signo negativo antes del término dispersivo indica que el contaminante se mueve hacia la zona de concentración más baja. De forma similar, las expresiones en las otras dos direcciones son expresadas como

$F_y = \vec{v}_ynC - nD_y\frac{\partial C}{\partial y}$ (A10.5)

$F_z = \vec{v}_znC - nD_z\frac{\partial C}{\partial z}$ (A10.6)

La cantidad total de soluto entrando al elemento cúbico (Figura A10.1) es

$F_x dz \hspace{1mm} dy + F_y dz \hspace{1mm} dx + F_z dx \hspace{1mm} dy$

La cantidad total que deja el cubo es

$\left(F_x + \frac{\partial F_x}{\partial x}dx\right)dz \hspace{1mm} dy \left(F_y + \frac{\partial F_y}{\partial y}dy\right)dz \hspace{1mm} dx + \left(F_z + \frac{\partial F_z}{\partial z}dx \hspace{1mm} dy\right)$

donde los términos parciales indican el cambio espacial de la masa del soluto en la dirección especificada. La diferencia en la cantidad entrando y saliendo del elemento es, por tanto

$\left(\frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} + \frac{\partial F_z}{\partial z}\right) dx \hspace{1mm} dy \hspace{1mm} dz$

A causa de que la distancia disuelta es asumida como no reactiva, la diferencia entre el flujo dentro de un elemento y el flujo fuera del elemento equivale a la cantidad de sustancia disuelta acumulada en el elemento. La tasa de cambio de masa en el elemento es

$-n\frac{\partial C}{\partial t} dx \hspace{1mm} dy \hspace{1mm} dz$

La expresión completa de conservación de la masa, por tanto, se convierte en

$\frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} + \frac{\partial F_z}{\partial z} = n\frac{\partial C}{\partial t}$ (A10.7)

Sustituyendo las expresiones (A10.4), (A10.5), y (A10.6) en (A10.7) y cancelando la n de ambos campos, se obtiene

$\bigg[\frac{\partial}{\partial x}\left(D_x\frac{\partial C}{\partial x}\right) + \frac{\partial}{\partial y}\left(D_y\frac{\partial C}{\partial y}\right) + \frac{\partial}{\partial z}\left(D_z\frac{\partial C}{\partial z}\right)\bigg] -$ $\bigg[\frac{\partial}{\partial x}(\vec{v}_xC) + \frac{\partial}{\partial y}(\vec{v}_yC) + \frac{\partial}{\partial z}(\vec{v}_zC)\bigg] = \frac{\partial C}{\partial t}$ (A10.8)

En un medio homogéneo en el cual la velocidad es estacionaria y uniforme (es decir, no varía a lo largo del tiempo o el espacio), los coeficientes de dispersión D_x, D_y, y D_z no varían a lo largo del espacio, y la Ec. (A10.8) se vuelve

$\bigg[D_x\frac{\partial^2C}{\partial x^2} + D_y\frac{\partial^2C}{\partial y^2} + D_z\frac{\partial^2C}{\partial z^2} \bigg] - \bigg[ \vec{v}_x\frac{\partial C}{\partial x} + \vec{y}_x\frac{\partial C}{\partial y} + \vec{v}_z\frac{\partial C}{\partial z} \bigg] = \frac{\partial C}{\partial t}$ (A10.9)

En una dimensión,

$D_x\frac{\partial^2C}{\partial x^2} - \vec{v}_x\frac{\partial C}{\partial x} = \frac{\partial C}{\partial t}$ (A10.10)

En algunas aplicaciones, la dirección unidimensional es tomada como una coordenada curvilínea en la dirección del flujo a lo largo de la línea de flujo. La ecuación de transporte se vuelve, entonces

$D_l\frac{\partial^2C}{\partial l^2} - \vec{v}_l\frac{\partial C}{\partial l} = \frac{\partial C}{\partial t}$ (A10.11)

Donde l es la dirección de la coordenada a lo largo de la línea de flujo, D_l el coeficiente longitudinal de la dispersión, y $\vec{v}_l$ , la velocidad lineal promedio a lo largo de la línea de flujo.

Para un problema bidimensional, es posible definir dos direcciones de coordenadas curvilíneas, S_l y S_t, donde S_l es dirigida a lo largo de la línea de flujo y S_t es ortogonal con respecto a esta. La ecuación de transporte luego se vuelve

$D_l\frac{\partial^2C}{\partial S_l^2} + D_t\frac{\partial^2C}{\partial S_t^2} - \vec{v}_l\frac{\partial C}{\partial S_l} = \frac{\partial C}{\partial t}$ (A10.12)

Donde D_l y D_t son los coeficientes de dispersión en las direcciones longitudinales y transversales, respectivamente

Si v_l varía a lo largo de la línea de flujo y D_l y D_t varían a lo largo del espacio, la Ec. (A10.13) se vuelve

$\frac{\partial}{\partial S_l}\left(D_l\frac{\partial C}{\partial S_l}\right) + \frac{\partial}{\partial S_t}\left(D_l\frac{\partial C}{\partial S_t}\right) - \frac{\partial}{\partial S_l}(\vec{v}_lC) = \frac{\partial C}{\partial t}$ (A10.13)

Las ecuaciones (A10.8) a (A10.13) representan seis formas de la ecuación de dispersión-advección para el transporte de solutos en un medio poroso saturado. La ecuación (A10.11) es idéntica a la Ec. (9.3) en la sección 9.2. La solución a cualquiera de estas dos ecuaciones proporcionará la concentración del soluto C como una función del espacio y el tiempo. Tomará la forma de C(x, y, z, t) para las Ecs. (A10.8) y (A10.9); C(l, t) para la Ec. (A10.11); y C(S_l, S_t, t) para las Ecs. (A10.12) y (A10.13). Existen muchas soluciones analíticas muy bien conocidas a la forma simple de transporte en las ecuaciones. La ecuación (9.5) en la Sección 9.2 es una solución analítica a la Ec. (A10.11), y la Ec. (9.6) es una solución analítica a la Ec. (A10.9). En la mayor parte de situaciones de campo, análisis de dos o tres dimensiones son normalmente requeridos. Además, las velocidades son raramente uniformes y las dispersividades normalmente varían en el espacio. Para estas condiciones, métodos numéricos adaptados para computadores digitales deben ser usados para obtener soluciones.

El coeficiente de dispersión en un medio tridimensional, homogéneo, isotrópico y poroso, como aquel expresado en la Ec. (A10.12), es un tensor simétrico de segundo orden con nueve componentes. D_l y D_t son términos diagonales en la forma bidimensional. Las propiedades direccionales de los coeficientes de dispersión están causadas por la naturaleza del proceso dispersivo en unas direcciones longitudinales y transversales. Si el medio por sí mismo es anisotrópico, la descripción matemática del proceso dispersivo tomaría mucha más complejidad. Para los sistemas anisotrópicos no hay soluciones analíticas o numéricas disponibles. Solo para el medio isotrópico se tiene el coeficiente de dispersión representado exitosamente por métodos experimentales.

Es posible extender la ecuación de transporte para incluir los efectos del retardo de un transporte de soluto debido a adsorción, reacciones químicas, transformaciones biológicas, o decaimiento radioactivo. En este caso, el balance de masa llevado a cabo en el volumen elemental debe incluir un término de fuente-sumidero. Para retardo debido a la adsorción, la ecuación de transporte en un medio homogéneo en un sistema unidimensional a lo largo de la dirección de flujo toma la forma de

$D_l\frac{\partial^2C}{\partial l^2} - \vec{v}_l\frac{\partial c}{\partial l} + \frac{\rho b}{n}\frac{\partial S}{\partial t} = \frac{\partial C}{\partial t}$ (A10.14)

Donde ρb es la densidad aparente de masa del medio poroso, n la porosidad, y S la masa del constituyente químico adsorbido en una unidad de masa de la parte sólida del medio poroso. El primer término de la Ec. (A10.14) es el término de dispersión, el segundo es el término de advección y el tercero es el término de reacción. Esta forma del término de reacción es considerada en la Sección 9.2 en conexión con las Ecs. (9.9) a (9.13). Soluciones analíticas a la Ec. (A10.14) son presentadas por Codell et al. (in press.). La Figura 9.12 presenta algunas soluciones tipo llevadas a cabo por Pickens & Lennox (1976) usando métodos numéricos.