Flujo por redistribución de agua a través de perfiles de la zona no saturada bajo condiciones de uso de la tierra agrícola en la cuenca alta del río Páez, Cartago, Costa Rica

Flow by redistribution of water through profiles of the unsaturated zone under conditions of agricultural land use in the Páez upper watershed, Cartago, Costa Rica

Fluxo por redistribuição de água através de perfis da zona não saturada em condições de uso do solo agrícola na bacia superior do rio Páez, Cartago, Costa Rica

Pablo Ramírez-Granados1

Universidad Nacional, Costa Rica

Karolina Villagra-Mendoza2

Instituto Tecnológico de Costa Rica, Costa Rica

Federico Masís-Meléndez3

Instituto Tecnológico de Costa Rica, Costa Rica

Introducción

El movimiento del agua en el suelo tiene un papel primordial en la agricultura ya que permite entender algunos fenómenos como la disponibilidad de agua en el suelo, el transporte de contaminantes como agroquímicos y el flujo de nutrientes como el nitrógeno y el fósforo (Ramadas et al. 2016).

La dinámica del flujo de agua en el suelo está basada en dos procesos principales, la infiltración y la redistribución. La infiltración es un proceso por el cual el agua entra en el suelo, generalmente, a través de la frontera entre la atmósfera y el suelo, lo que produce un flujo hacia abajo que cambia los gradientes de contenido de humedad y presión de agua-poro (Sonaje, 2013).

Una vez finalizado el proceso de la infiltración, inicia un movimiento de agua hacia abajo a través del perfil del suelo llamado Redistribución de Agua del Suelo (RAS) (Kirkham, 2005). Tal proceso puede llegar a persistir por un tiempo conforme la humedad fluya dentro del perfil del suelo (Kendy et al. 2003, Chen et al. 2019).

En el RAS las principales variables que controlan el movimiento del flujo no saturado son el potencial gravitacional y el potencial mátrico (ψ) (Biswas et al. 1966; Gardner et al. 1970), los cuales permiten el flujo de agua a través de gradientes de carga hidráulica entre los horizontes. La RAS en la ZNS tiende a formar un frente de humedecimiento producto de la intrusión de aire desde la superficie (Raats, 1973) y la tortuosidad de medio (Wang et al. 2003), el cual es un proceso complicado debido a la inestabilidad que se da en el flujo (Jury et al. 2003) y a la histéresis (Staple, 1966). Ese último es el elemento más relevante en el flujo de agua (Rubin, 1967) además que es causado principalmente por el efecto de botella, el ángulo de contacto, el aire atrapado y la contracción e hinchazón del suelo (Horgan y Ball 2005; Zhai et al. 2021). Otros autores como Philip y Knight (1991) agregan que la relativa escasez de los estudios físico-matemáticos de la RAS se debe a las complicadas condiciones iniciales y a lo infructífero de la aplicación de las ecuaciones de la infiltración en la redistribución.

El movimiento por RAS en suelos no saturados puede ser explicado a través de los mecanismos de flujo que se pueden dar en él. Por ejemplo, el flujo constante y el flujo intermitente (Bouma, 1981). Dentro de estos flujos se incluyen otros como el flujo de pistón que ocurre lentamente a través de una gran parte del espacio poroso (Wiekenkamp et al. 2020) y el flujo preferencial, el cual se refiere al movimiento vertical del agua a lo largo de los macroporos o a través de una matriz no saturada (Van Stiphout et al. 1987) el cual puede ser mediado por los macroporos o por la matriz (Eguchi y Hasegawa 2008).

Para analizar el flujo por RAS en la ZNS se han aplicado diversos métodos. Charbeneau (1989) usó modelos físicos simples basados en continuidad y en la ley de Darcy, Philip (1992) utilizó soluciones para ecuaciones de difusión-convección no lineal, Wang y Zhang (2011) analizaron el flujo de agua en los macroporos usando trazadores y colorantes, Mavinbela y Rensburg (2011) utilizaron mediciones experimentales con monolitos, Ramadas et al. (2016) usaron soluciones analíticas para calcular la evolución de perfiles de contenido volumétrico y los flujos de agua, Lai et al. (2017) utilizaron el modelo de infiltración con RAS de Green-Ampt,, Chen et al. (2019) usaron un modelo conceptual bajo redistribución de la capacidad de campo; y Fuentes et al. (2020) aplicaron la resolución numérica de la ecuación vertical de Richards.

El presente trabajo tiene como objetivo analizar la RAS en perfiles de la ZNS en suelos dedicados a las actividades agrícolas en una parte de la zona norte de Cartago, dentro de la cuenca alta del río Páez. Esto debido a que en las zonas agrícolas existe un enorme consumo de agua y una menor eficiencia (Ortiz et al. 1999), además de la temática ligada a la aplicación de agroquímicos (Ramírez et al. 2014). La hipótesis de este trabajo es que a diversas profundidades existen aportes específicos de agua producto de los mecanismos de flujo por RAS en los horizontes del perfil.

Materiales y métodos

Descripción del área de estudio

Para la medición de la RAS, se instalaron tres monolitos (M1, M2 y M3), distribuidos en la parte alta de la cuenca del río Páez en Cartago, específicamente, en los distritos de Potrero Cerrado y Santa Rosa del cantón de Oreamuno (ver Figura 1). La cuenca tiene una precipitación media anual de 1600 mm, y una temperatura media anual de 13,8 °C. Se caracteriza por estar dentro de la zona de mayor productividad de hortalizas de Costa Rica con alrededor de un 80% (Ramírez et al. 2008) además que es una importante zona de producción de agua para consumo humano debido a los manantiales que se encuentran en la subcuencas de los ríos Reventado, Tiribí, Tatiscú, Páez, Birrís y Pacayas (Arellano et al. 2012). Sin embargo, existen importantes problemas originados por las actividades agrícolas como el alto uso de agroquímicos (Ramírez et al. 2014), el consumo de agua debido a los cultivos hortícolas (Narváez, 2013), los sistemas de riego y los requerimientos hídricos de los cultivos (Quesada, 2017) y la erosión por mecanización (Mehuys et al. 2009; Gómez et al. 2017) que pueden llegar a afectar la dinámica del agua en los suelos.

Las zonas donde se ubican los monolitos se caracterizan por estar en áreas de uso agrícola, donde se alternan los pastos y los cultivos. Dentro de los cultivos más importantes están la papa, el brócoli, la zanahoria, la coliflor, el repollo, la remolacha y el culantro.

Los materiales volcánicos sobre los cuales se desarrollan los suelos donde se encuentran los monolitos fueron descritos a partir de algunas perforaciones realizadas en el sector de Prusia, Sanatorio Durán, parque Prusia y la Pastora (Arellano et al. 2012), todos muy cercanos a la cuenca alta.

Los suelos son andisoles pertenecientes a los subórdenes de los orthents y udands en la parte más alta; udands, en la parte media; y ustands en la parte baja según el mapa de suelos de Mata y Sandoval (2016).

Figura 1: Distribución de los tres monolitos instalados de acuerdo con su ubicación distrital en la cuenca alta del río Páez,
cantón de Oreamuno, Cartago.

Fuente: Elaboración propia, 2022.

El Cuadro 1 muestra la ubicación geográfica y la altitud de cada monolito, así como el suborden de suelo y las dimensiones utilizadas para su construcción. Los tres monolitos fueron instalados en parcelas dedicadas a la producción agrícola y a la cobertura de pastos en reposo.

Cuadro 1. Ubicación y dimensiones de los monolitos construidos en la cuenca alta del río Páez

Fuente: Elaboración propia, 2022.

Construcción de los monolitos

La construcción de los monolitos se basó en una variante del método del perfil instantáneo o drenaje interno (Hillel et al. 1972; Kranz et al. 1998; Krisdani et al. 2009; Hari et al. 2010; Abbasi et al. 2012; Pirastru et al. 2017) con el fin de medir las variaciones del contenido volumétrico del agua (θ) en distintas profundidades de forma continua (Vachaud y Dane, 2002; Reynolds, 2008 B) para determinar los aportes de agua en el perfil en diferentes tiempos.

Para obtener θ se instalaron 9 reflectómetros de dominio de tiempo (TDR) modelo CS650 con 30 cm de guía (Campbell Scientific, Logan, UT) en una de las caras de cada monolito (Figura 2). La ubicación de cada sensor respondió a un arreglo rectangular en tres niveles, a 25, 50 y 75 cm de profundidad para el monolito 1 (M1) y el monolito 2 (M2); y hasta 50 cm para el monolito 3 (M3) debido a la aparición de material rocoso. En cada nivel se colocaron tres sensores a 50 cm de distancia horizontal entre ellos. Un registrador de datos, CR1000x (Campbell Scientific, Logan, UT), registró la lectura de cada sensor. El tamaño de cada monolito se estableció de acuerdo con un volumen elemental representativo (VER) de la ZNS (Figura 2).

Figura 2: Esquema de la instalación de los sensores θ en cada uno de los monolitos. La imagen ejemplifica el arreglo de sensores en el M1.

Fuente: Elaboración propia, 2022.

Cada monolito fue cubierto con una película de plástico transparente para mantener la estabilidad del monolito y evitar ingresos de agua lateral asumiendo los lados impermeables. Finalmente, se realizó una segunda cubierta con plástico negro durante el tiempo del ensayo para evitar la evaporación y la entrada de agua en el horizonte superior.

Las lecturas de los sensores se obtuvieron durante un período continuo de 12,5 días entre los meses de agosto y noviembre del año 2019. Las variaciones en θ fueron registradas y analizadas de acuerdo con el procedimiento descrito por Vachaud y Dane (2002) y Reynolds (2008 B). La cantidad de agua almacenada (S) en cada horizonte entre las profundidades 0 y z en el tiempo t fue calculada con la siguiente ecuación:

(1)

donde S = cantidad de agua almacenada (L), θ = contenido volumétrico de agua (L3/L3), z = profundidad de medición (L), t = tiempo (T).

El mecanismo de flujo fue identificado a partir de la interpolación de θ usando el método de kriging (Rubin, 2003), el cual es aplicado al cálculo de la superficie continúa formada por θ:

(2)

donde Ẑ (x0) = variable interpolada, 𝜆i = pesos, xi = localización geográfica, Z (xi) = valores muestreados, m = media.

Para la identificación de los tipos de flujos se visualizaron las isolíneas de θ en tres tiempos específicos, 180, 7200 y 18300 minutos en cada uno de los 3 monolitos.

Propiedades físicas e hidráulicas de los horizontes

Para cada horizonte se obtuvieron muestras inalteradas utilizando cilindros con un volumen aproximado de 100 cm3 y se recolectaron muestras disturbadas de 2 kg, tamizadas a través de una malla de 2 mm y secadas al aire. La textura se obtuvo a partir del método del hidrómetro de Bouyoucos (Kroetsch y Wang, 2008) y las curvas granulométricas fueron determinadas por tamizado en húmedo (Kroetsch y Wang, 2008). La materia orgánica (MO) fue obtenida por el método de redox dicromato (Skjemstad y Baldock, 2008). La densidad aparente (ρa) fue obtenida con el método del núcleo (Henríquez y Cabalceta, 1999; Hao et al. 2008) y la densidad real (ρr) con el método del matraz aforado (Henríquez & Cabalceta, 1999). La porosidad (φt) se calculó con la siguiente ecuación según Hao et al. (2008):

(3)

siendo φt = porosidad total, ρa = densidad aparente [M/L3], ρr = densidad real [M/L3].

La conductividad hidráulica saturada (Ksat) fue obtenida mediante carga hidráulica constante (Henríquez y Cabalceta, 1999; Reynolds, 2008 A) y luego calculada por la ecuación según Reynolds (2008 A):

(4)

donde V, [L3], es el volumen de agua durante el intervalo de tiempo Δt [T], L [L] es la longitud del núcleo, ΔH [L] es la diferencia en elevación entre el nivel de agua en el cilindro adicional y la salida del cilindro de muestreo y dc [L] es el diámetro del cilindro.

Las curvas de retención de humedad del suelo (ψ-θ) fueron obtenidas por medio de la extracción de presión (Reynolds y Topp 2008), la cual fue ajustada usando el modelo de van Genuchten (van Genuchten 1980) con el programa RECT (van Genuchten et al. 1980) utilizando la ecuación siguiente:

(5)

donde θs = humedad de saturación [L3/L3], θr = humedad volumétrica residual [L3/L3], ψ = potencial mátrico [L] y α [1/L], n y m = constantes.

Resultados y discusión

Efecto de las propiedades físicas en el movimiento del agua en el suelo

Los suelos están compuestos de partículas de tamaño variable, derivadas de la meteorización de las rocas y la transformación de los minerales primarios. Estos materiales conforman los horizontes de un perfil y son el medio por el cual fluye el agua luego de ser infiltrada. Para poder analizar las propiedades de los materiales, se realizaron dos análisis, un análisis textural para visualizar la distribución de las partículas totales; y, otro granulométrico, para observar las variaciones en la composición de la fracción más gruesa, la cual es dominante en estos suelos volcánicos. A continuación, se presentan las descripciones de los perfiles para cada monolito.

Monolito 1

● Ap (1a): gris oscuro, grano suelo, arenoso, agregados débiles subangulares, raíces muy finas, friable a muy friable, no adhesivo y arenoso.

● A2 (1b): capa de ceniza gris oscura, textura muy arenosa, muy friable, grano suelto, con poca o ninguna estructura y no adherente.

● A3 (1c): suelo ligeramente rojizo, friable, grano suelto con algunos agregados subangulares, de poca consistencia, textura al tacto arenosa y no adherente.

● A4 (1d): gris claro, grano suelto, formación de algunos agregados subangulares, friable textura arenosa y no adherente.

● A5 (1e): capa de ceniza, textura arenosa, gris claro, sin consistencia en seco y muy poca en húmedo al muestreo, no adherente, muy friable, grano suelo, sin estructura ni agregados.

● A6 (1f): gris claro a amarillento, formación de agregados subangulares, adherencia de baja a muy baja, friable y textura arenosa.

Monolito 2

● Ap (2a): Marrón oscuro, friable, grano suelto con formación de agregados subangulares, textura al tacto arenosa y poca adherencia.

● A2 (2b): Marrón claro, friable, grano suelto con algunos agregados subangulares, moderada consistencia, textura arenosa al tacto y adherencia de moderada a baja.

Monolito 3

● Ap (3a): Marrón claro, poco friable, grano suelto con algunos agregados subangulares, moderada consistencia a fuerte, textura arenosa al tacto, adherencia baja y presencia de material rocoso.

● A2 (3b): Marrón oscuro, poco friable, grano suelto con algunos agregados subangulares, moderada consistencia a fuerte, textura arenosa al tacto, adherencia baja y presencia de material rocoso.

● C1 (3c): Rojo oscuro, estructura masiva a moderada, con agregados subangulares, fuerte consistencia, textura arenosa al tacto y adherencia de moderada a baja.

● C2 (3d): Rojo oscuro, gravoso, desarrollo estructura fuerte, no plástico, consistencia extremadamente firme, textura arenosa al tacto y adherencia de moderada a baja.

La disposición de los horizontes de cada uno de los monolitos y su espesor se presenta en la Figura 3:

Figura 3: Disposición de los horizontes en cada uno de los perfiles para cada monolito.

Fuente: Elaboración propia, 2022.

El Cuadro 2 muestra la textura y el contenido de materia orgánica (MO) en los horizontes identificados para cada monolito. Se observa que, a lo largo de los perfiles de los tres monolitos, la textura es muy uniforme, la cual se clasificó como franco-arenosa (FAO, 2009). La MO tiende a ser mayor en el primer horizonte para todos los monolitos, sin embargo, conforme se profundiza en el perfil, se empieza a disminuir y luego a aumentar nuevamente.

Cuadro 2: Fracciones texturales y MO de los horizontes en cada monolito

Fuente: Trabajo de campo y laboratorio, 2020.

Las texturas concuerdan con las reportadas para la zona (Vega, 1999; Mata y Sandoval, 2016) y su origen se da a partir de depósitos volcánicos, especialmente cenizas volcánicas (Yatno y Suharta, 2011). Esta uniformidad textural está relacionada también con los materiales tobáceos subyacentes (Arellano et al. 2012) con los materiales volcánicos descritos en campo en sitios cercanos (Ramírez, 2007; Arellano et al. 2012; Zúñiga y Ramírez, 2015). La razón del incremento de la MO en profundidad en M1 y M3 está relacionada a los restos vegetales anteriores al cambio de vegetación (Kramer y Gleixner, 2008) que se han dado en esta zona, además del tipo de labranza de volteo (Espinoza, 2010). La MO tiene influencia significativa en la retención de agua del suelo (Ankenbauer y Loheide, 2016) debido a su afinidad con el agua, Yang et al. (2014) mencionan que en ψ menos negativos, la MO afecta la retención de agua del suelo principalmente alterando los parámetros estructurales, mientras a ψ más negativos la MO incrementa la capacidad absorbente del suelo para retener más agua.

El Cuadro 3 muestra la composición detallada de la fracción de arenas para cada monolito.

Cuadro 3: Detalle de la composición de la fracción arenosa (gruesa, media y fina) en porcentaje, para los monolitos

Fuente: Trabajo de campo y laboratorio, 2020.

En M1 hay variaciones en la fracción arenosa en todos los horizontes, mientras en el M2 2 y el M3 hay dominio de la fracción de arena gruesa en todos los horizontes, lo cual es esperable para suelos desarrollados bajo capas vulcanoclásticas (Pérez, 2000), estos materiales volcánicos son altamente permeables (Tejedor et al. 2003), debido a la presencia de macroporos (>75 μm). Los suelos arenosos muestran varios rangos de φ y están relacionados con la presencia de numerosos poros que resultan del arreglo de las partículas de arenas (Braund et al. 2004), de la actividad de la fauna y el desarrollo de las raíces (Braund et al. 2005) y los cambios en poros estructurales y biológicos que pueden ser modificados por compactación (Cuolon y Bruand, 1989). La relación de los poros con las fracciones arenosas radica en que estos pueden ayudar a acelerar la infiltración de agua (Zhou et al. 2013). Arya y Paris (1981) mencionan la relación entre el diámetro de poro y el agua retenida por presión capilar de acuerdo con la distribución granulométrica, la cual es un aspecto importante en el movimiento de agua a través de los horizontes.

El Cuadro 3 presenta ρa y ρr, así como la φ total para cada uno de los horizontes identificados en los monolitos M1, M2 y M3. En M1 se observa un comportamiento fluctuante de la ρa conforme se profundiza en el perfil del suelo. Por otro lado, en M2 y M3 la ρa tiende a disminuir conforme aumenta la profundidad del suelo. Consecuentemente, la φ total tiende a disminuir en M1 y a aumentar en M2 y M3 conforme se avanza en los horizontes del suelo.

Cuadro 4: Densidades y porosidades de los horizontes de los monolitos M1, M2 y M3.

Fuente: Trabajo de campo y laboratorio, 2020.

Los cambios en profundidad de la ρa obtenidos fueron reportados en otros perfiles (Alvarado y Forsythe, 2005) y están condicionados por el uso agrícola intenso que luego afecta a las propiedades hidráulicas del perfil (Assouline 2006; Dec et al. 2008). En suelos de composición arenosa como es el caso de los horizontes de los monolitos de la cuenca alta del río Páez, tanto la ρa como la φ están afectados por la uniformidad del empaquetamiento, la cantidad de compactación y la distribución del tamaño de grano (Tindall y Kunkel, 1999).

Los valores de ρr reportados para los horizontes de los monolitos están acordes con los valores para los suelos andisoles (Shoji et al. 1993) así como los cambios en profundidad (Vega, 1999). En M1, las diferencias en la profundidad den las ρr podrían asociarse a las diferentes composiciones mineralógicas de los materiales que originan los suelos, ya que en esta zona está el área de contacto entre las formaciones Sapper y Reventado Miembro Superior (Krushensky, 1972). En M2, los materiales son derivados de la Formación Birrís por lo que se mantienen sus ρr en profundidad, a la vez que en M3, las variaciones en las densidades están relacionadas a las que se encuentran en la Formación Cervantes. Esta última presenta gran parte de su superficie como roca desnuda y, por lo tanto, variaciones significativas en el desarrollo del suelo. En el caso de φ, el promedio de los valores de los horizontes de cada monolito es alto y su relación con la profundidad es evidente en M2 y M3, donde hay relaciones altas entre ambas variables y ninguna relación con el M1. La diferencia en M1 está relacionada a los horizontes de cenizas, de la misma forma que se vio en las anteriores relaciones. La φ puede variar en profundidad entre límites muy amplios debido a los distintos porcentajes de las fracciones de poros (Aranda et al. 1980) que se pueden dar especialmente en la fracción arenosa y a los cambios por el grado de fragmentación, compactación e inversión de la zona de labranza por el equipo de labranza, la compactación por tráfico y la coalescencia de la zona de labranza durante períodos de humedecimiento (Kay y VandenBygaart, 2002).

Efecto de las propiedades hidráulicas en el movimiento del agua

El Cuadro 5 detalla la Ksat para de M1, M2 y M3. En general, se observa un comportamiento irregular en el perfil de M1, con variaciones promedio de 87% entre horizontes, mientras que la Ksat aumenta entre los horizontes en el orden de 460% para M2 y 115% para M3.

Cuadro 5: Conductividad hidráulica saturada de los monolitos

Fuente: Trabajo de campo y laboratorio, 2020.

La Ksat es altamente dependiente de θ del suelo (Ng y Shi, 1998) y sus cambios en profundidad son causados por los efectos de la preparación de la tierra para el cultivo (Green et al. 2003) y la heterogeneidad de las propiedades del suelo entre los horizontes, lo que causa que algunas regiones sean más o menos favorables para el flujo (Sobieraj et al. 2004) así como los efectos de la estructura del suelo (Ghassemi y Pak, 2011) y la φ (Braund et al. 2005). Los valores de la Ksat se clasifican como moderados de acuerdo con la clasificación de Hazelton y Murphy (2016).

Las curvas ψ-θ de los tres monolitos en todos sus horizontes se muestran en la Figura 4. En M1, M2 y M3, específicamente en el horizonte 1, el θs se mantiene hasta un ψ más negativo (-33 kPa) en comparación con los horizontes más profundos. Esto incide en la Ksat, la cual llega a mantenerse a un valor máximo, con un ψ más negativos que en los perfiles más profundos, lo que mejora las condiciones del movimiento del agua en el suelo. El mayor de θs en la superficie de M2 y M3 con respecto a M1 puede ser debido a un mayor contenido de MO.

Figura 4: Curvas ψ-θ para los horizontes de los monolitos a) M1, b) M2 y c) M3 ajustadas con la ecuación de van Genuchten usando el modelo de conductividad de Mualem en el software RECT (van Genuchten et al. 1980).

Fuente: Elaboración propia, 2022.

En las curvas ψ-θ en los monolitos (Figura 4), la distribución del tamaño de grano está relacionada con la retención de humedad (Yang et al. 2004). La estructura de cada uno de los horizontes, especialmente la heterogeneidad, (Abbasi et al. 2012) así como el efecto de la MO (Salcedo et al. 2007), juegan un rol importante. El comportamiento de las curvas ψ-θ para los tres monolitos concerniente a elementos como los minerales con alta superficie especifica como los amorfos y los óxidos de Fe y Al (Salcedo et al. 2007) asociados a los suelos andisoles (Alvarado et al. 2014) y el laboreo intensivo y continuo (Ceballos et al. 2010; Gómez y Estrada, 2020) que provoca cambios estructurales en el suelo relacionados con la distribución de poros de este.

Variación del θ y dinámica del agua

La Figuras 5, la Figura 6 y la Figura 7 muestran el registro continuo del θ a tres diferentes profundidades para M1, durante 12,5 días.

En a (Figura 5), a una profundidad de 25 cm, el θ inicial fue de 0,122 cm3/cm3, el cual disminuyó a 0,103 cm3/cm3 en un lapso de 12,74 días, lo que dio una RAS de 4,75 mm. A 50 cm, el θ máximo fue de 0,150 cm3/cm3, el cual era el θ inicial al comienzo de la RAS en esa profundidad. Este θ disminuyó hasta un valor de 0,139 cm3/cm3 al final del tiempo de medición. Durante el lapso de 12,74 días, θ aumentó dos veces; a 6,92 y 10,90 días, el θ pasó de 0,141 a 0,145 cm3/cm3; así, la RAS en cada uno de los dos incrementos fue de 2 mm. A 75 cm, el θ inicial fue de 0,360 cm3/cm3, a 6,50 días, este pasó de 0,356 a 0,361 cm3/cm3; luego, a 10,77 días, pasó de 0,357 a 0,361 cm3/cm3. La RAS fue de 3,75 mm y 3,0 mm. a 75 cm para cada uno de los incrementos.

Figura 5: Comportamiento del θ en el perfil 1 para el monolito 1 en cada z medida en campo.

Fuente: Trabajo de campo, 2019.

En b (Figura 6), a 25 cm, el θ inicial era 0,083 cm3/cm3, el cual disminuyó hasta 0,068 cm3/cm3 en 12,74 días, lo que generó una RAS de 3,75 mm. A 50 cm de profundidad, al inicio de la medición, el θ era de 0,171 cm3/cm3, luego, disminuyó hasta 0,155 cm3/cm3 y solo se observó un incremento del θ a 6,58 días al pasar de 0,161 a 0,165 cm3/cm3, lo cual dio una RAS de 2,0 mm, la cual disminuyó hasta 0,161 cm3/cm3. A 75 cm, el valor del θ inicial fue de 0,365 cm3/cm3, el cual fue reduciéndose en el tiempo de medición y tuvo dos incrementos por RAS a 6,56 y 10,83 días; el primero, desde 0,360 a 0,367 cm3/cm3, lo cual significó un incremento de 5,25 mm; y el segundo, a 10,83 días, de 0,362 a 0,367 cm3/cm3, lo que representó un aporte por RAS de 3,75 mm.

Figura 6: Comportamiento del θ en el perfil 2 para el monolito 1 en cada z medida en campo.

Fuente: Trabajo de campo, 2019.

En c (Figura 7), el θ a 25 cm inicialmente tuvo un valor máximo de 0,103 cm3/cm3 con un valor final de 0,087 cm3/cm3 al concluir la medición, lo que generó un aporte de 4,0 mm en 12,74 días. A 50 cm se observaron dos incrementos debido a los aportes desde 25 cm: a 6,60 días el incremento pasó de 0,166 a 0,172 cm3/cm3, lo que generó una RAS de 3,0 mm; luego a 11,04 días, θ pasó de 0,166 cm3/cm3 a 0,171 cm3/cm3, lo que se traduce en una RAS de 2,5 mm. A 75 cm se observaron 3 aportes por RAS, a 4,54 días hubo un primer incremento de 0,361 a 0,366 cm3/cm3, lo cual dio una RAS de 3,75 mm; luego, un segundo aporte a los 6,60 días de 0,362 a 0,372 cm3/cm3 que generó una RAS de 7,5 mm; y un último incremento a 9,98 días de 0,364 a 0,371 cm3/cm3 que significó una RAS de 5,25 mm.

Figura 7: Comportamiento del θ en el perfil 3 para el monolito 1 en cada z medida en campo.

Fuente: Trabajo de campo, 2019.

En el M2 se presentan las variaciones en θ en cada perfil en 12,5 días. En a (Figura 8), el θ a 25 cm varió desde 0,466 a 0,367 cm3/cm3 en un lapso de 12,81 días, el cual se tradujo en una entrada de agua de 24,75 mm. A 50 cm, el θ mostró el mismo comportamiento de 25 cm. El θ cambió desde 0,489 a 0,374 cm3/cm3, lo cual significó una RAS de 57,5 mm. A 75 cm, el θ inicial fue de 0,491 cm3/cm3, luego, a 9,92 días, θ varió desde 0,421 a 0,428 cm3/cm3, lo que implicó una RAS de 5,25 mm.

Figura 8: Comportamiento del θ en el perfil 1 para el monolito 2 en cada z medida en campo.

Fuente: Trabajo de campo, 2019.

En b (Figura 9), el θ a 25 cm tuvo un valor inicial de 0,444 cm3/cm3, el cual disminuyó hasta 0,372 cm3/cm3 en un periodo de 12,81días que significó una entrada de agua de 18 mm. En 50 cm, el valor inicial de θ fue de 0,456 cm3/cm3, el cual se redujo hasta un valor de 0,397 cm3/cm3 en 10,13 días, luego hubo un incremento de θ a 0,400 cm3/cm3 que en RAS significó 1,5 mm. A 75 cm, θ fue de 0,472 cm3/cm3, luego, a 9,81 días, se dio otro incremento desde 0,418 a 0,426 cm3/cm3 que dio una RAS de 6 mm.

Figura 9: Comportamiento del θ en el perfil 2 para el monolito 2 en cada una z medida en campo.

Fuente: Trabajo de campo, 2019.

En c (Figura 10), el θ en 25 cm fue de 0,457 cm3/cm3, el cual decreció hasta 0,409 cm3/cm3 durante el tiempo de medición, lo que en RAS significó un aporte de 12 mm. A 50 cm, el θ, al momento de la medición inicial, fue de 0,458 cm3/cm3; a 9,90 días, hubo un aporte desde 0,406 a 0,419 cm3/cm3 que produjo una RAS de 6,5 mm. A 75 cm, el θ fue de 0,467 cm3/cm3 que fue disminuyendo en 9,81 días hasta llegar a 0,434 cm3/cm3, el cual recibió un aporte por RAS de 9,75 mm.

Figura 10: Comportamiento del θ en el perfil 3 para el monolito 2 en cada z medida en campo.

Fuente: Trabajo de campo, 2019.

En el M3 se presentan las diferencias en θ en cada perfil para un tiempo de 12,5 días. En a (Figura 11), la variación en θ a 25 cm, durante 13,77 días, fue de 0,421 hasta 0,380 cm3/cm3 al final del tiempo de medición, cuyo aporte de agua fue de 10,25 mm. A 37,5 cm, el θ fue de 0,473 cm3/cm3, el cual mostró en su drenaje durante la medición tres aportes por RAS; a 0,73 días, el θ paso de 0,463 a 0,502 cm3/cm3, el cual dio en 14,625 mm; luego, a 5,63 días, el incremento pasó de 0,452 a 0,476 cm3/cm3, lo que generó 9 mm y una RAS de 3,75 mm debido al incremento desde 0,447 a 0,457 cm3/cm3 en el tiempo de 10,77 días. A 50 cm, a 0,70 días hubo un primer aporte de 0,472 a 0,492 cm3/cm3, lo que dio una RAS de 10 mm, luego, se presentó un segundo incremento a 2,04 días de 0,47 a 0,474 cm3/cm3 que resultó una RAS de 2 mm. Asimismo, se dio un tercer incremento desde los 2,8 días, de 0,47 a 0,475 cm3/cm3 que generó una RAS de 2,5 mm. A los 6,81 días se dio otra RAS de 6 mm, lo cual sucedió por el incremento desde 0,466 a 0,478 cm3/cm3. Además, a los 10,85 días, hubo un último incremento desde 0,462 a 0,466 cm3/cm3 para una RAS de 2 mm.

Figura 11: Comportamiento del θ en el perfil 1 para el monolito 3 en cada z medida en campo.

Fuente: Trabajo de campo, 2019.

En b (Figura 12), a 25 cm el θ fue de 0,502 cm3/cm3, el cual disminuyó durante 13,77 días hasta 0,443 cm3/cm3, cuyo aporte de agua se devino en 14,75 mm. A 37,5, el θ fue de 0,473 cm3/cm3 al que durante su drenaje se le observaron 4 períodos de RAS desde la capa superior. El primer incremento se dio a 0,73 días desde 0,465 a 0,490 cm3/cm3, este generó una RAS de 9,375 mm; el segundo ocurrió incremento a los 2,19 días desde 0,464 a 0,478 cm3/cm3 que significó 5,25 mm de RAS; el tercero a los 5,60 días, el θ cambio de 0,4595 a 0,4815 cm3/cm3 que dio una RAS de 8,25 mm; y el último incremento fue a los 10,73 días desde 0,456 a 0,4655 cm3/cm3, lo que implicó una RAS de 3,5625 mm. A 50 cm, el θ en el suelo fue de 0,516 cm3/cm3, el cual, conforme fue disminuyendo, experimentó 5 aportes por RAS. A los 0,75 días hubo un primer incremento desde 0,511 a 0,549 cm3/cm3, este significó 19 mm de RAS; luego a 2,08 días, el θ ascendió desde 0,506 hasta 0,541 cm3/cm3, lo cual implicó una RAS de 17,5 mm; después a los 5,96 días se dio un tercer incremento desde 0,503 a 0,512 cm3/cm3 que en RAS significó 4,5 mm; a los 9,42 días, el θ tuvo un incremento de 0,496 a 0,502 cm3/cm3 que significó un valor de RAS de 3 mm; y por último, a los 10,71 días, la capa tuvo un aporte desde 0,495 a 0,501 cm3/cm3 que implicó una RAS de 3 mm, igual al incremento anterior.

Figura 12: Comportamiento del θ en el perfil 2 para el monolito 3 en cada z medida en campo.

Fuente: Trabajo de campo, 2019.

En c (Figura 13), el θ a 25 cm fue de 0,420 cm3/cm3 que descendió hasta 0,389 cm3/cm3, lo que generó una entrada de agua de 7,75 mm en 13,77 días. En 37,5 cm, el θ fue de 0,482 cm3/cm3, en cuya profundidad se observaron 5 aportes por RAS en los siguientes tiempos: 0,69; 1,98; 2,90; 5,75 y 10,46 días. Para el primer aporte, el θ pasó de 0,474 a 0,494 cm3/cm3, lo cual fue 7,5 mm de RAS; el segundo aporte fue de 0,470 a 0,493 cm3/cm3 que generó un incremento de 8,625 mm; el tercero fue de 7,5 mm, el cual resultó en el cambio de 0,474 a 0,494 cm3/cm3; el cuarto incremento fue desde 0,487 a 0,467 cm3/cm3 que resultó en una RAS de 7,5 mm; mientras aporte de 0,465 a 0,473 cm3/cm3 que concluyó en una RAS de 3 mm. A 50 cm, el θ inicial fue de 0,474 cm3/cm3, este empezó a drenar y experimentó un incremento a 0,73 días desde 0,467 a 0,489 cm3/cm3, que en RAS, significó 11 mm; luego a los 2,02 días, hubo otro incremento desde 0,466 a 0,480 cm3/cm3 que implicó un aporte de 7 mm; a los 2,90 días un nuevo incremento por RAS permitió un cambio desde 0,469 a 0,489 cm3/cm3, cuyo RAS fue de 8 mm; a los 5,11 días, un aporte desde 0,463 a 0,471 cm3/cm3 que significó un aumento de 4 mm; y por último, a los 11,10 días, hubo un cambio desde 0,458 a 0,462 cm3/cm3, que dio como resultado 2 mm de RAS.

Figura 13: Comportamiento del θ en el perfil 3 para el monolito 3 en cada z medida en campo.

Fuente: Trabajo de campo, 2019.

En cuanto a los tipos de flujo de agua que se dan en cada uno de los monolitos, la interpolación θ (z,t), mediante kriging se permitió visualizar el comportamiento de θ (z,t) y el tipo de flujo.

Mecanismos de redistribución en los monolitos

En M1 (Figura 14), el comportamiento del flujo es de tipo pistón donde el agua infiltrada va redistribuyéndose en el perfil de cada monolito de forma que es observable un desplazamiento del frente de humedad con un comportamiento casi uniforme, con dispersión y mezcla despreciable.

Figura 14: Visualización del tipo de flujo de agua en el monolito 1 a partir de la interpolación de θ (z,t) de los sensores.

Fuente: Elaboración propia, 2022.

En M2 (Figura 15), el flujo de agua infiltrada presenta un flujo por RAS en el perfil de tipo preferencial orientado hacia la derecha en la cara del monolito. En este caso, es observable una dirección del flujo de humedad orientado lateralmente.

Figura 15: Visualización del tipo de flujo de agua en el monolito 2 a partir de la interpolación de θ (z,t) de los sensores.

Fuente: Elaboración propia, 2022.

En M3 (figura 16), se pueden visualizar dos tipos de flujos, inicialmente un flujo que presenta un comportamiento de pistón en los primeros 40 cm; luego, otro que comienza a seguir un flujo preferencial hacia el centro de la base del monolito.

Figura 16: Visualización del tipo de flujo de agua en el monolito 3 a partir de la interpolación de θ (x,t) de los sensores.

Fuente: Elaboración propia, 2022.

Los horizontes de los monolitos, a pesar de tener homogeneidad como capas dada su composición textural, pueden presentar anisotropía horizontal debido a condiciones de humedad inicial (Philip, 1991), lo cual origina comportamientos diferentes por pequeñas diferencias entre secciones húmedas y secas (Zhuang et al. 2017) en cada horizonte, así como cambios temporales (Schwen et al. 2014). Durante el drenaje, el θ gradualmente disminuye con el tiempo en todos los monolitos debido a las texturas arenosas que presentan (Hillel y Talpaz 1977; Ortiz et al. 1999).

El flujo por RAS está relacionada a la tortuosidad del flujo a través de la configuración de los poros y los pasajes de flujo (Ghassemi y Pak 2011) de acuerdo con la disposición del medio poroso (Koponen et al. 1997) en su textura y estructura

Las variaciones en las curvas de descenso de θ están relacionadas a pequeñas variaciones laterales en la succión capilar (Kung, 1990). Los cambios en θ son graduales y muy pequeños en cada uno de los sensores debido al flujo, con disminuciones lentas en el tiempo. La diferenciación de los flujos de pistón y preferencial es difícil de separar por las pequeñas variaciones en θ en el suelo (Zhang et al. 2019 A).

En los monolitos, los frentes de humedad identificados presentaron dos comportamientos diferentes. En M1 y M3, el flujo de pistón observado obedece a un desplazamiento uniforme del frente de humedecimiento a través de las profundidades medidas, hecho relacionado a que ambos presentan una mayor cantidad de horizontes con respecto a M2, lo que hace que el flujo sea constante y vaya remplazando el θ en los poros (Bouma, 1981) sin que se manifieste un flujo preferencial específico. Yao y Hendrickx (1996) afirman que los frentes de humedecimiento se estabilizan bajo tasas de infiltración bajas. Así en el caso del M1, las pequeñas variaciones en θ identificadas durante el proceso de medición y el número de horizontes con granulometrías distintas contribuyen a la estabilidad del flujo por RAS en los perfiles de los monolitos. En M2, el flujo de agua acorde con las variaciones en θ es de tipo preferencial, identificado por una mayor convergencia lateral (Jarvis, 2007), donde los mayores valores están orientados hacia una dirección específica, seguidos por la absorción lateral (Bouma, 1981). Gao et al. (2014) mencionan que el uso de la tierra y el θ inicial podrían incrementar la velocidad del frente de humedecimiento y, con ello, el flujo preferencial. El M2 presentaba pastos con ganado, lo cual incrementa el desarrollo de macroporos debido a las raíces y cambios en las características y la geometría de la red porosa estructural (Mossadeghi et al. 2016), lo que lleva a que se desarrolle flujos no uniformes (Gerke et al. 2010). En el caso de los árboles, el incremento del flujo preferencial por las raíces ha sido mencionado por otros autores (Zhang et al. 2019 B), en el caso del M2, si bien otros autores (Kan et al. 2019) indican que en pastos el flujo de pistón es más dominante, los cambios constantes en el uso de la tierra pueden llegar a presentar un terreno descubierto y, por lo tanto, cambios en la dinámica del flujo de agua en profundidad con el desarrollo de la vegetación.

Con respecto al M3, su comportamiento dual con flujos de pistón al inicio y, preferencial después, ha sido evidenciado previamente (Zhang et al. 2019 A; Xiang et al. 2019). El M3, a diferencia de los otros, presentó al inicio un comportamiento de flujo de pistón y, conforme aumentó la cantidad de material piroclástico en profundidad, se incrementaron las fracturas y por lo tanto el desarrollo de pasos preferenciales para el flujo vertical (Šimůnek et al. 2003).

Conclusiones

El presente trabajo tuvo como objetivo presentar un primer estudio al proceso de RAS en el suelo a partir de la experimentación en campo con sensores de θ. Durante su ejecución, la hipótesis de trabajo fue confirmada, donde se pudo contabilizar los aportes de agua en profundidad producto del flujo por RAS en los horizontes de perfiles de la ZNS.

El flujo por RAS es un proceso complejo y depende muchísimo de las propiedades físicas e hidráulicas de los suelos que, a su vez, también presentan una enorme variabilidad espacial de acuerdo con las condiciones locales que se tengan.

Durante el proceso de RAS en los perfiles de la ZNS, se pudieron encontrar los mecanismos que tienen lugar dentro de él, como lo son los flujos de pistón y preferencial.

El mecanismo de flujo varía de un monolito a otro y pudo ser identificado a partir de la interpolación de los θ de acuerdo con su posición en cada uno de los perfiles.

La determinación del flujo de agua por RAS será de utilidad en el estudio de los procesos de flujo de agua que se dan en las zonas agrícolas, entre ellos, se tienen el flujo de agroquímicos y pesticidas, así como la recarga de agua. El entendimiento de este proceso podría responder por qué podrían migrar los contaminantes agrícolas hacia las capas acuíferas y el desarrollo y recuperación de la vegetación.

Referencias

Abbasi, F., Javauaxm M., Vanclooster, M. y Feyen, J. 2012. Estimating hysteresis in the soil water retention curve from monolith experiments. Geoderma 189-190: 480-490.

Alvarado, A. y Forsythe, W. (2005). Variación de la densidad aparente en órdenes de suelos de Costa Rica. Agronomía Costarricense 29(1), 85-94.

Alvarado, A., Mata, R. y Chinchilla, M. (2014). Arcillas identificadas en suelos de Costa Rica a nivel generalizado durante el período 1931-2014: I. Historia, metodología de análisis y mineralogía de arcillas en suelos derivados de cenizas volcánicas. Agronomía Costarricense 38(1), 75-106.

Ankenbauer, K. y Loheide, S. (2016). The effects of soil organic matter on soil water retention and plant water use in a meadow if the Sierra Nevada, CA. Hydrological Processes (31), 891-901.Aranda, M., Arrúe, J., Ugarte; J. y Moreno, F. (1980). Propiedades físicas y químicas en relación con la porosidad de los suelos: Influencia del ciclo natural de humectación-desecación. IV Hinchamiento y contracción. Anales de Edafología y Agrobiología 39 (3-4), 523-537

Aranda, M; Arrúe, J; Ugarte; J; Moreno, F. (1980). Propiedades físicas y químicas en relación con la porosidad de los suelos: Influencia del ciclo natural de humectación-desecación. IV Hinchamiento y contracción. Anales de Edafología y Agrobiología 39 (3-4): 523-537.

Arellano, F., Vásquez, M., Suárez, J. y Rafael, E. (2012). Identificación de zonas de protección y de captura de las nacientes, subcuencas de los ríos Reventado, Tiribí, Tatiscú, Páez, Birrís y Pacayas, Provincia de Cartago, Costa Rica. Licitación Abreviada 2012LA-000068-87900 Informe Final. COMCURE-MINAET-Dirección de Agua-Hidrogeotecnia Ltda.

Arya, L. y Paris, J. (1981). A physicoempirical model to predict the soil moisture characteristic from particle-size distribution and bulk density data. Soil Science Socienty of America Journal (45), 1023-1030.

Assouline, S. (2006). Modeling the relationship between soil bulk density and the water retention curve. Vadose Zone Journal (5), 554-563.

Biswas, T., Nielsen, D. y Biggar, J. (1966). Redistribution of soil water after infiltration. Water Resources Research 2(3), 513-524.

Bouma, J. (1981). Soil morphology and preferential flow along macropores. Agricultural Water Management (3), 235-250.

Braund, A., Hartmann, C., Ratana, S., Sindhusen, P., Poss, R. y Hardy, M. (2004). Composition, fabric, and porosity of an arenic haplustalf of Northeast Thailand: Relation and penetration resistance. Soil Science Society of America Journal (68), 185-193.

Braund, A., Hartmann, C. y Lesturgez, G. (2005). Physical properties of tropical sandy soils: A large range of behaviours. Management of Tropical Sandy Soil for Sustainable Agriculture. International Union of Soil Science.

Ceballos, D., Hernández, O. y Vélez, J. (2010). Efecto de la labranza sobre las propiedades físicas en un andisol del Departmento de Nariño. Revista de Agronomía 23(1), 40-48.

Charbeneau, R. (1989). Liquid moisture redistribution: Hydrologic simulation and spatial variability. En Morel-Seytoux (ed). (1989). Unsaturated Flow in Hydrologic Modeling: Theory and Practice. Kluwer Academic Publishers, pp.127-161.

Chen, S., Huo, Z., Xu, X. y Huang, G. (2019). A conceptual agricultural water productivity model considering under field capacity soil water redistribution applicable for arid and semi-arid areas with deep groundwater. Agricultural Water Management (213) 309-323.

Cuolon, E. y Bruand, A. (1989). Effects of compaction on the pore space geometry in sandy soils. Soil & Tillage Research (15), 137-152.

Dec, D., Dörner, J., Becker, O. y Horn, R. (2008). Effect of bulk density on hydraulic properties of homogenized and structured soils. Revista de Ciencia del Suelo y Nutrición Vegetal 8(1), 1-13.

Eguchi, S. y Hasegawa, S. (2008). Determination and characterization of preferential water flow in unsaturated subsoil of andisol. Soil Science Society of America Journal 72(2), 320-330.

Espinoza, Y. (2010). Efecto de la labranza sobre la materia orgánica y tamaño de agregado en un suelo cultivado con maíz en condiciones tropical. Bioagro 22(3), 177.184.

FAO. (2009). Guía para la Descripción de Suelos. 4ta edición. FAO.

Fuentes, S., Trejo, J., Quevedo, A., Fuentes, C. y Chávez, C. (2020). Modeling soil water redistribution under gravity irrigation with Richards Equation. Mathematics 8(9), 1581.

Gao, Z., Xu, X., Zhao, J., Yu, M., Zhang, S. y Zhao, C. (2014). Field study on preferential flow under different land uses in the Loess Hilly region. Nature Environment and Pollution Technology 13(3), 483-490.

Gardner, W., Hillel, D. y Benyamini, Y. (1970). Post-irrigation movement of soil water: 1. Redistribution. Water Resource Research 6(3), 851-861.

Gerke, H., Germann, P. y Nieber, J. (2010). Preferential and unstable flow: From the pore to the catchment scale. Vadose Zone Journal 9(2), 207-212.

Ghassemi, A. y Pak, A. (2011). Pore scale study of permeability and tortuosity for flow through particulate media using Lattice Boltzmann method. International Journal for Numerical and Analytical Methods in Geomechanics 35(8), 886-901.

Gómez, N., Solórzano, M. y Villagra, K. (2017). Cuantificación de la erosión hídrica en función de diferentes técnicas de mecanización para minimizar la contaminación del agua por sedimentos en la parte alta de la cuenca del río Reventazón. Informe Final. Tecnológico de Costa Rica.

Gómez, N. y Estrada, R. (2020). Conservación de suelos mediante la modificación de la frecuencia de labranza: Un caso en Costa Rica. Revista de Ciencias Ambientales 54(1), 123-139.

Green, T., Ahuja, L. y Benjamin, J. (2003). Advances and challenges in predicting agricultural management effects on soil hydraulic properties. Geoderma (116), 3-27.

Hao, X., Ball, B., Culley, J., Carter, M. y Parkin, G. (2008). Soil Density and Porosity. En Carter, M. y Gregorich, E. (eds). (2008). Soil Sampling and Methods of Analysis pp. 743-760. 2da edición. Canadian Society of Soil Science.

Hari, K., Ojha, C., Chandramouli, P. y Madramootoo, C. (2010). Estimation of unsaturated hydraulic parameters from infiltration and internal drainage experiments. Journal of Irrigation and Drainage Engineering.

Hazelton, P. y Murphy, B. 2016. Interpreting Soil Test Results: What do all the numbers mean? 3era edición. CSIRO Publishing.

Henríquez, C. y Cabalceta, G. (1999). Guía Práctica para el Estudio Introductorio de los Suelos con un Enfoque Agrícola. Asociación Costarricense de la Ciencia del Suelo.

Hillel, D., Krentos, V. y Stylianou, Y. (1972). Procedure and test of an internal drainage method for measuring soil hydraulic characteristics in situ. Soil Science 114(5), 395-400.

Hillel, D. y Talpaz, H. (1977). Simulation of soil water dynamics in layered soils. Soil Science 123 (1), 54-62.

Horgan, G. y Ball, B. (2005). Modelling the effect of water distribution and hysteresis on air-filled pore space. European Journal of Soil Science (56), 647-634.

Jarvis, N. (2007). A review of non-equilibrium water flow and solute transport in soil macropores: principles, controlling factors and consequences for water quality. European Journal of Soil Science (58), 523-546.

Jury, W., Wang, Z. y Tuli, A. (2003). Conceptual model of unstable flow in unsaturated soil during redistribution. Vadose Zone Journal (2), 61-67.

Kan, X., Cheng, J., Hu, X., Zhu, F. y Li, M. (2019). Effects of grass and forest and the infiltration amount on preferential flow in karst regions of China. Water 11 (8): 1634.

Kay, B. y VandenBygaart, A. (2002). Conservation tillage and depth stratification of porosity and soil organic matter. Soil & Tillage Research (66), 107-118.

Kendy, E., Gérard, P., Walter, M.,Zhang, Y., Liu, C. y Steenhuis, T. (2003). A soil-water-balance approach to quantify groundwater recharge from irrigated cropland in the North China Plain. Hydrological Processes 17, 2011-2031.

Kirkham, M. (2005). Principles of Soil and Plant Water Relations. Elsevier.

Koponen, A., Kataja, M. y Timonen, J. (1997). Permeability and effective porosity of porous media. Physical Review E 56(3), 3319-3325.

Kramer, G. y Gleixner, G. (2008). Soil organic matter in soil depth profiles: Distinct carbon preferences of microbial groups during carbon transformation. Soil Biology and Biochemistry (40), 425-433.

Kranz, W., Kanwar, R. y Pederson, C. (1998). Collection and monitoring of one-meter cubic soil monoliths for leaching studies. Transactions of the American Society of Agricultural Engineers 41(2): 333-344.

Krisdani, H., Rahardjo, H. y Leong, E. (2009). Use of instantaneous profile and statistical methods to determine permeability functions of unsaturated soils. Canadian Geotechnical Journal (46), 869-874.

Kroetsch, D. y Wang, C. (2008). Particle Size Distribution. En Carter, M. y Gregorich, E. (eds). (2008). Soil Sampling and Methods of Analysis pp. 713-725. 2da edición. Boca Raton, Canadian Society of Soil Science.

Krushensky, R. (1972). Geology of Istaru Quadrangle, Costa Rica. Bulletin 1358. United States Geological Survey.

Kung, K. (1990). Preferential flow in a sandy vadose zone: 1. Field observation. Geoderma (46), 51-58.

Lai, W., Ogden, F., Steinke, R. y Talbot, C. (2017). An efficient and guaranteed stable numerical method for continuous modeling of infiltration and redistribution with a shallow dynamic water table. Water Resources Research (51), 1514-1528.

Mata, R. y Sandoval, D. (2016). Mapa Digital de Suelos de Costa Rica. Centro de Investigaciones Agronómicas. Universidad de Costa Rica.

Mehuys, G., Tiessem, K., Villatoro, M., Sancho, F. y Lobb, D. (2009). Erosión por labranza con arado de disco en suelos volcánicos de ladera en Costa Rica. Agronomía Costarricense 33(2), 205-215.

Mavinbela, S. y Rensburg, L. (2011). In-situ evaluation of internal drainage in layered soils (Tukulu, Sepane and Swartland). Hydrological Earth System Science 8: 9797-9841.

Mossadeghi, M., Arvidsson, J., Keller, T., Koestel, J., Lamandé, M., Larsbo, M. y Jarvis, N. (2016). Effects of subsoil compaction on hydraulic properties and preferential flow in a Swedish clay soil. Soil & Tillage Research (156), 91-98.

Narváez, I. (2013). Percepción sobre la tendencia de caudales, precipitación, temperatura y cambio de uso del suelo con relación al uso y manejo del agua en la zona norte de Cartago, Costa Rica. (Trabajo de Graduación). Centro Agronómico Tropical de Investigación y Enseñanza.

Ng, C. y Shi, Q. (1998). A numerical investigation of the stability of unsaturated soi slopes subjected to transient seepage. Computers and Geotechnics 22(1). 1-28.

Ortiz, J., Gavrilov, I., Palacios, O. y Acosta, R. (1999). Pérdidas de agua de riego por percolación profunda durante el proceso de infiltración. Terra Latinoamericana 17(2), 115-124.

Pérez, F. (2000). The influence of surface volcaniclastic layers form Haleakala (Maui, Hawaii) on soil water conservation. Catena (38) 301-332.

Philip, J. (1991). Horizontal redistribution with capillary hysteresis. Water Resource Research 27(7): 1459-1469.

Philip, J. y Knight, J. (1991). Redistribution of soil water from plane, line, and point sources. Irrigation Science (12), 169-180.

Philip, J. (1992). Exact solutions for redistribution by nonlinear convection-diffusion. Journal of Australian Mathematical Society (33) 362-383.

Pirastru, M., Marrosu, R., Di Prima, S., Keesstra, S., Giadrossich,F. y Niedda, M. (2017). Lateral saturated hydraulic conductivity of soil horizons evaluated in large-volume soil monoliths. Water 9(11), 862.

Quesada, J. (2017). Determinación de los requerimientos de agua para diferentes cultivos de la zona norte de Cartago. (Trabajo Final de Graduación). Tecnológico de Costa Rica.

Ramadas, M., Ojha, C. y Govindaraju, R. (2016). Analytical models of infiltration and redistribution for unsaturated flow in soils with vertically non-uniform saturated hydraulic conductivity. ISH Journal of Hydraulic Engineering, 22:2, 181-192.

Ramírez, P. (2007). Caracterización de la dinámica de flujo mediante la aplicación de un modelo numérico hidrogeológico: Caso de la cuenca del río Birrís, Cartago, Costa Rica. Revista Geológica de América Central 34-35: 83-97.

Ramírez, L., McHugh, A y Alvarado, A. (2008). Evolución histórica y caracterización socioeconómica de la cuenca media del río Reventado, Cartago, Costa Rica. Agronomía Costarricense 32(2), 53-72.

Ramírez, F., Fournier, M., Ruepert, C. y Hidalgo, C. (2014). Uso de agroquímicos en el cultivo de papa en Pacayas, Cartago, Costa Rica. Agronomía Mesoamericana 25(2), 337-345.

Raats, P. (1973). Unstable wetting fronts in uniform and nonuniform soils. Soil Science Society of America Proceedings (37), 681-685.

Reynolds, W. (2008 A). Saturated Hydraulic Properties: Laboratory Methods. En Carter, M. y Gregorich, E. (eds). (2008). Soil Sampling and Methods of Analysis pp. 1013-1024. 2da edición. Boca Raton, Canadian Society of Soil Science.

Reynolds, W. (2008 B). Unsaturated Hydraulic Properties: Instantaneous Profile. En Carter, M. y Gregorich, E. (eds). (2008). Soil Sampling and Methods of Analysis pp. 1129-1138. 2da edición. Boca Raton, Canadian Society of Soil Science.

Reynolds, W. y Topp, G. (2008). Soil Water Desorption and Imbibition: Tension and Pressure Techniques. En Carter, M. y Gregorich, E. (eds). (2008). Soil Sampling and Methods of Analysis pp. 981-997. 2da edición. Boca Raton, Canadian Society of Soil Science.

Rubin, J. (1967). Numerical method for analyzing hysteresis-affected, post-infiltration redistribution of soil moisture. Soil Science Society of America Proceedings 31: 13-20.

Rubin, Y. (2003). Applied Stochastic Hydrogeology. Oxford University Press.

Salcedo, E., Galvis, A., Hernández, T., Rodríguez, R., Zamora, F., Bugarin, R. y Carrillo, R. (2007). La humedad aprovechable y su relación con la materia orgánica y superficie específica de suelo. Terra Latinoamérica 25(4), 419-425.

Schwen, A., Zimmermann, M. y Bodner, G. (2014). Vertical variations of soil hydraulic properties within two soil profiles and its relevance for soil water simulations. Journal of Hydrology (516), 169-181.

Shoji, S., Nanzyo, M. y Dahlgren, R. (1993). Volcanic Ash Soils: Genesis, Properties and Utilization. Developments in Soil Science 21. Elsevier.

Šimůnek, J., Jarvis, N., van Genuchten, M. y Gärdenäs, A. (2003). Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. Journal of Hydrology (272), 14-35.

Skjemstad, J. y Baldock, J. (2008). Total and Organic Carbon. En Carter, M. y Gregorich, E. (eds). (2008). Soil Sampling and Methods of Analysis pp. 225-237. 2da edición. Boca Raton, Canadian Society of Soil Science.

Sobieraj, J., Elsenbeer, H. y Cameron, G. (2004). Scale dependency in spatial patterns of saturated hydraulic conductivity. Catena (55), 49-77.

Sonaje, N. (2013). Modeling of infiltration process – A review. Indian Journal of Applied Research 3(9) 226-230.

Staple, W. (1966). Infiltration and redistribution of water in vertical columns of loam soil. Soil Science Society of America Proceedings (30), 553-558.

Tejedor, M., Jiménez, C. y Díaz, F. (2003). Volcanic materials as mulches for water conservation. Geoderma (117), 283-295.

Tindall, J. y Kunkel, J. (1999). Unsaturated Zone Hydrology for Scientist and Engineers. Prentice Hall.

Van Genuchten, M. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society American Journal (44), 892-898.

Van Genuchten, M., Leij; F. y Yates, S. (1980). The RECT code for quantifying the hydraulic functions of unsaturated soils. U.S. Salinity Laboratory. U.S. Department of Agriculture. Agriculture Research Service.

Van Stiphout, T., van Lanen, H., Boersma, O. y Bouma, J. (1987). The effect of bypass flow and internal catchment of rain on the water regime in a clay loam grassland soil. Journal of Hydrology (95). 1-11.

Vachaud, G. y Dave, J. (2002). Instantaneous Profile. En Dane, J. y Topp, G. (eds). (2002). Methods of Soil Analysis, Part 4 Physical Methods pp. 937-945. Soil Science Society of America.

Vega, O. (1999). Estudio semidetallado de suelos de la zona de Tierra Blanca, Potrero Cerrado y Llano Grande, Cartago. Ministerio de Agricultura y Ganadería.

Wang, Z., Wu, L., Harter, T., Lu, J. y Jury, W. (2003). A field study of unstable preferential flow during soil water redistribution. Water Resources Research 34(4), 1075.

Wang, K. y Zhang, R. (2011). Heterogeneous soil water flow and macropores described with combined tracers of dye and iodine. Journal of Hydrology (397) 105-117.

Wiekenkamp, I., Huisman, J., Bogema, H. y Vereecken, H. (2020). Effects of deforestation on water flow in the vadose zone. Water 12 (1), 2-16.

Xiang, W., Si, B., Biswas, A y Li, Z. (2019). Quantifying dual recharge mechanisms in deep unsaturated zone of Chinese Loess Plateau using stable isotopes. Geoderma (337), 773-781.

Yang, H., Rahardjo, H., Leong, E. y Fredlund, D. (2004). Factors affecting drying and wetting soil-water characteristic curves of sandy soils. Canadian Geotechnical Journal (41), 908-920.

Yang, F., Zhang, G., Yang, J., Li, D., Zhao, Y., Liu, F., Yang, R. y Yang, F. (2014). Organic matter controls of soil water retention in an alpine grassland and its significance for hydrological processes. Journal of Hydrology (519), 3086-3093.

Yao, T. y Hendrickx, J. (1996). Stability of wetting fronts in dry homogeneous soils under low infiltration rates. Soil Science Society of America Journal (60), 20-28.

Yatno, E. y Suharta, N. (2011). Andisols derived from pyroclastic liparite tuff: Their properties and their management strategy for agricultural development. Jurnal Tanah Dan Iklim (33), 49-64.

Zhai, Q., Rahardjo, H., Satyanaga, A., Zhu, Y., Dai, G. y Zhao, X. (2021). Estimation of wetting hydraulic conductivity function for unsaturated sandy soil. Engineering Geology 285: 106034.

Zhang, Z., Si, B., Li, H. y Li, M. (2019 A). Quantify piston and preferential water flow in deep soil using Cl- and soil water profiles in deforested apple orchards on the Loess Plateau, China. Water 11(10), 2183.

Zhang, J., Lei, T., Qu, L., Zhang, M., Chen, P., Gao, X., Chen, C. y Yuan, L. (2019 B). Method to quantitatively partition the temporal preferential flow and matrix infiltration in the forest soil. Geoderma 347(10), 150-159.

Zhou, B., Li, Y., Wang, Q., Jiang, Y. y Li, S. (2013). Preferential water and solute transport through sandy soil containing artificial macropores. Environmental Earth Science (70), 2371-2379.

Zhuang, L., Hassanizadeh, S., Kleingeld, P. y van Genuchten, M. (2017). Revisiting the horizontal redistribution of water in soils: Experiments and numerical modeling. Water Resources Research 53: 75.

Zúñiga, H. y Ramírez, P. (2015). Geología de la parte norte del Cantón Central de Cartago, Costa Rica. Informe de consultoría. Hidrogeotecnia Ltda.

---