Por: Jorge Zafra, Principal Hydrologist Minera Las Bambas; Luis Suescún, Senior Hidrogeologist Modeler, y Robin Dufour, director general de DHI Perú.

Resumen

La minería a tajo abierto es una actividad económica relevante en varios países del mundo, incluyendo Perú. Sin embargo, puede tener un impacto significativo en el medio ambiente y en la calidad del agua subterránea, especialmente si no se toman medidas adecuadas de manejo y monitoreo de la hidrogeología local. Por esta razón, los estudios hidrogeológicos son esenciales para comprender la fase operacional de la minería a tajo abierto y los efectos que con la ejecución se pueden presentar en los recursos hídricos superficiales y subterráneos, y para diseñar estrategias de mitigación, optimización y operación de los proyectos mineros.

En este contexto, se han desarrollado diferentes metodologías y herramientas para modelar la hidrogeología de las labores mineras a tajo abierto, y mejorar la precisión y eficiencia de los modelos numéricos que se utilizan para realizar estas estimaciones de impactos y avances operacionales. Entre estas metodologías se encuentra la modelación inversa, una herramienta clave en la hidrogeología que permite abordar la incertidumbre en los modelos numéricos y conceptuales de agua subterránea. Sin embargo, los modeladores a menudo ejecutan y analizan separadamente la calibración estacionaria y transitoria debido a restricciones operacionales en los códigos numéricos y al desconocimiento de metodologías de trabajos eficientes en los proyectos mineros. 

Por esta razón, se han explorado nuevas formas de optimizar el pre y post procesamiento, así como acoplar la corrida de modelos estacionarios y transitorios en la fase de calibración, con el fin de ahorrar tiempo de modelamiento y contar con mayor tiempo para tareas de evaluación de los resultados, análisis de incertidumbre, etc.

Para abordar este desafío, se analizó una importante mina de cobre a cielo abierto del tajo Ferrobamba de la mina Las Bambas en Perú, que utilizó un modelo geológico en Leapfrog y numérico 3D en Feflow para representar la geología e hidrogeología compleja de la región andina en su condición de preminado y actual. Del modelo hidrogeológico conceptual, trabajos de campo y pruebas hidráulicas ejecutadas permitieron identificar zonas de reducción de conductividad hidráulica a profundidad, zonas fracturadas y kársticas, que llevaron a una parametrización de más de 1,500 zonificaciones. Se emplearon métodos de programación de IFM en Python para minimizar los tiempos de corrida de las instalaciones mineras y técnicas de factorización, puntos piloto y multiplicadores, en un acoplamiento de calibración estacionaria-transitorio en Pest para minimizar los tiempos de cálculo y apoyar la evaluación de impactos ambientales requeridos por el EIA.

Palabras clave: Modelo hidrogeológico, Mina a tajo abierto.

Introducción

Para el inicio del análisis y con el objetivo de identificar la extensión de la zona de estudio se procedió a delimitar el área de cuencas hidrogeológicas que forman parte del sistema a evaluar.

Las cuencas hidrogeológicas se definieron con base en la superficie freática, unidades hidrogeológicas y direcciones de flujo. En la mayoría de los casos las cuencas hidrogeológicas coinciden con cuencas de aguas superficiales. Sin embargo, en algunos casos, divergen, debido a los sistemas kársticos que pueden tener flujo entre cuencas. Estos se ubican en la cuenca del manantial Challhuapuquio, la cuenca Chinchina ubicada al sur de Sulfobamba, Pumamarca, situada al norte de Ferrobamba, la cuenca del río Tambo y las cuencas kársticas dentro del río Ferrobamba.

Se desarrollaron conceptos de balance hídrico para la roca fracturada, la cubierta del cuaternario y los sistemas de caliza kárstica y cómo se relacionan entre sí. Asimismo, la distribución de la precipitación se estimó en base a datos de descarga de manantiales y arroyos. En general, se estima alrededor del 30% a 45% de pérdidas por evapotranspiración, entre 10% a 40% de escorrentía, y alrededor del 30% de interflujo para suelos cuaternarios y sistemas de rocas fracturadas. El sistema fracturado que alcanza el nivel freático está entre 10% a 20% de la precipitación total para los suelos cuaternarios y sistemas de rocas fracturadas y solo alrededor del 12% del agua que se infiltra en el sistema kárstico es retenida por el sistema fracturado y microkárstico de movimiento más lento.

Las tasas de drenaje se estimaron en modelos de hoja de cálculo basados en datos y experiencia para sistemas cuaternarios, de roca fracturada y kársticos como porcentaje de recarga para los meses siguientes al evento de recarga. El drenaje rápido ocurre como interflujo y flujo a través de los conductos kársticos, la matriz y el flujo de fractura son un drenaje constante y los sistemas mixtos tienen una tasa de drenaje intermedia.

Se desarrolló un balance en hoja de cálculo con todos estos insumos hidrológicos que permite estimar el promedio mensual de recarga y descarga. En general, el modelo de balance hídrico fue capaz de simular la estación seca completa y, a veces, la húmeda temprana.

Área de estudio

La unidad minera Las Bambas se encuentra ubicada en la Sierra sur del Perú, en la región Apurímac, entre los distritos de Challhuahuacho, Tambobamba y Coyllurqui, provincia de Cotabambas, y el distrito del Progreso, provincia de Grau, a una altitud que varía entre los 3,800 y 4,600 msnm., aproximadamente 75 km al suroeste de la ciudad de Cusco. El área de la unidad minera Las Bambas se encuentra dentro de la intercuenca del Alto Apurímac que drena sus aguas hacia la vertiente del Atlántico. La mayor parte de la mina se ubica dentro de la subcuenca del río Ferrobamba y la cabecera de las subcuencas de los ríos Pamputa, Punanqui y Pumamarca (Figura 1).

Área de mina

Las Bambas es una mina de cobre que produce el 2% del suministro mundial del metal rojo y uno de los principales productores en el mundo (Surran, 2022). En la unidad minera se ubican tres yacimientos de los cuales uno de ellos está minado a cielo abierto y dos de ellos se planea que se integren junto con la operación en los siguientes años con la misma técnica de extracción mineral. Los depósitos de minerales de Las Bambas son tipo pórfido Skarn de Cu-Mo-(Ag-Au). Se encuentran en la franja de yacimientos de tipo pórfido y Skarn del SE del Perú, asociada al Batolito Andahuaylas - Yauri del Eoceno Oligoceno (Perelló et al., 2003) emplazado en rocas sedimentarias del Mesozoico, siendo la más importante la Formación Ferrobamba (Cretáceo inferior-superior).

Actualmente, el área de mina comprende el tajo Ferrobamba, el depósito de desmonte de material estéril (Huancarane), el depósito de relaves (Figura 2), los campamentos y demás instalaciones auxiliares, las cuales constituyen y tendrán conexión con una extensión de la vía férrea existente hacia el puerto de Matarani para su exportación (Percy, 2018).

Tajo Ferrobamba

El tajo Ferrobamba se ubica al sur del valle del río del mismo nombre, en la subcuenca 7, con valores de elevación entre los 3,800 y 4,400 msnm. (Piteau, 2010). Su operación comenzó en 2015 excavando en la cumbre de la Formación Ferrobamba sobre el lado norte del río. El tajo tiene una forma ligeramente ovalada a circular con un diámetro de unos 1,750 m a 2,000 m. Las elevaciones del fondo del tajo al final de operación se estiman que estén entre los 3,465 a 3,495 msnm., aproximadamente. La excavación y drenaje del tajo requirió de un canal de coronación del río Ferrobamba alrededor del tajo y el drenaje de la estribación de caliza Ferrobamba que fue minada (Golder, 2010).

Depósito de relaves

El depósito de relaves de la mina Las Bambas (TSF por su nombre en inglés: Tailing Storage Facility) se localiza en la cabecera de la cuenca del río Ferrobamba, a 8 km al noroeste de la población de Challhuahuacho, en la región Apurímac. De acuerdo con el informe de Modelación Numérica de la Expansión del TSF de mina Las Bambas (DHI 2021b) los potenciales receptores de interés para la infiltración en esta zona del depósito son la caliza de la Formación Ferrobamba al sur del depósito de relaves, el tajo Ferrobamba ubicado a 3.5 km aguas abajo del depósito y la cuenca limítrofe del río Tambo a 2.5 km al sur del depósito.

Depósito de desmonte Huancarane

Ubicado al noroeste del tajo Ferrobamba, está diseñado para contener alturas de desecho de roca de hasta 350 m, elevaciones superficiales de casi 4,220 msnm. y un volumen de relleno total de 424 Mm³ en un área aproximada de 4 km². En él se dispondrá del material extraído y excedente del tajo Ferrobamba, cuyo volumen se estima sea en aproximadamente 193 Mm³.

Las Bambas planea incluir una serie de componentes adicionales aumentando su producción entre 380 K y 400 K toneladas métricas/año de concentrado de cobre a mediano plazo (Surran, 2022), junto con la operación actual del tajo principal, cuyo plan de minado está proyectado hasta el 2043.

De la información geológica existente, datos observados de niveles piezométricos, etc., definición del plan de minado histórico y técnicas de implementación numéricas en Python con su administrador de interfaz Feflow - IFM (en inglés, Feflow Interface Manager) desarrollada por DHI, se realizó una implementación y calibración acoplada para determinar los efectos que tendrán los componentes actuales y futuros con la expansión en su próximo Estudio de Impacto Ambiental (EIA).

Metodología

Feflow Interface Manager

Feflow cubre una amplia gama de funciones para la simulación de flujo y transporte de medios porosos, accesible a través de una interfaz de usuario integral. No obstante, hay casos en los que se desea y resulta útil un control ampliado del usuario sobre los procesos internos del software. Los ejemplos típicos incluyen la importación de datos de fuentes definidas por el usuario, simulaciones combinadas con otro software, interdependencias extendidas entre parámetros, procesos físicos específicos y formatos de salida específicos del usuario. Muchos de estos casos pueden adaptarse al concepto de complemento Feflow y su interfaz de programación abierta controlada por IFM (Figura 3).

Si bien el uso de los complementos disponibles no requiere ninguna habilidad de programación, el IFM también brinda soporte para la codificación de complementos propios por parte de modeladores avanzados con cierta experiencia en programación en C/C++ y Python.

La columna vertebral de los IFM son las funciones API (en inglés, Application Programming Interface y están disponibles en https://dhi.github.io/ifm/), las cuales proporcionan un conjunto de métodos, funciones y protocolos que permiten a los modeladores:

1. Realizar operaciones simples y complejas dentro y fuera de la interfaz.

2. Optimizar el flujo de trabajo diario como modelador de aguas subterráneas.

3. Ahorrar tiempo y optimizar los tiempos de ejecución de una simulación.

4. Cambiar o leer los conjuntos de datos de procesamiento previo y posterior: antes, durante y después de una simulación, lo cual son procesos tediosos y repetitivos y consumen tiempo para el modelador en sus tareas diarias.

5. La automatización de los flujos de trabajo puede aumentar su eficiencia en los modelos numéricos de procesamiento previo y posterior.

6. Realizar operaciones dentro y fuera de la interfaz como la creación, modificación y ejecución de simulaciones. Esto puede incluir la capacidad de definir propiedades del medio poroso, condiciones de contorno, fuentes y sumideros etc., así como extraer resultados y analizar los datos generados por las simulaciones.

7. Acoplamiento con otras tecnologías (por ejemplo, aguas superficiales, etc.).

8. Acoplamiento de simulaciones en estado estacionario y transitorio y problemas de modelado inverso o calibración automática.

Modelo geológico numérico

FloSolutions (2021) elaboró un modelo hidrogeológico conceptual de la zona de estudio para la unidad minera. A partir de este, se elaboró la construcción del modelo geológico regional empleando técnicas de modelamiento implícito empleando el software Leapfrog Geo, aplicando controles al modelo consistente con la secuencia de unidades geológicas y sus orientaciones y relación sedimentaria: erosión o depositación o estructural – intrusión (Figura 5) y teniendo en cuenta la integración de la siguiente información:

1. Contexto geológico e interpretación basándose en el trabajo de recopilación geológico, realizado en la fase previa al inicio del modelamiento geológico (Figura 4).

2. Uso de mapas topográficos digitales Aster GDEM provenientes de la información abierta de la NASA.

3. Elaboración de secciones regionales a partir de mapas geológicos desarrollados en los trabajos de recopilación de información.

4. Revisión y unificación de criterios geológicos a partir de los modelos geológicos locales de los trabajos ubicados en las áreas de operación minera actuales y futuras.

5. Validación y revisión extensiva del modelo a partir de geólogos especialistas en geología regional y geólogos locales de Las Bambas.

6. Definición de sistemas de fallas principales (Falla Tambo, Falla Ferrobamba, Falla Pamputa y Falla Record), más de 100 sistemas de fallas representados con elementos discretos en el modelo numérico.

7. Definición de las principales unidades hidrogeológicas, tales como:

ν Intrusivos: todas las unidades intrusivas sin ninguna diferenciación usando contactos verticales: se distribuyen en gran parte de la unidad minera, condicionando el comportamiento hidráulico de la mayoría de los componentes mineros. Están compuestas principalmente de pórfidos de cuarzo- monzonita.

ν Skarn/mármol: afloran en superficie en sectores localizados en la zona de Chalcobamba y Sulfobamba a partir de la interacción entre unidades intrusivas y sedimentarias (calizas).

ν Formación Ferrobamba: afloran en gran parte de la unidad minera, siendo de principal importancia en la zona de Ferrobamba. Se compone de sucesiones de calizas micríticas.

ν Formación Mara: se compone de lutitas, y afloran a las afueras de la zona de influencia directa de los tajos.

ν Formación Soraya: se compone de cuarzoarenitas, y se distribuyen a las afueras de la zona de influencia directa de los tajos, principalmente hacia el sur de la unidad minera.

ν Formación Labra: se compone principalmente de areniscas y se distribuyen en el sur de la unidad minera.

ν Formación Cachíos: aflora a las afueras de la influencia directa de los tajos, en el sur y norte de la unidad minera. Se compone principalmente de lutitas.

ν Depósitos aluviales: estos se componen de acumulaciones de grava, arena, limo y arcilla con clastos subangulosos a angulosos de diferente composición. Afloran en toda la extensión de la unidad minera.

ν Depósitos Morrénicos: se componen de bloques angulosos rellenados con arcillas, limos y arenas. Afloran en toda la extensión de la unidad minera.

ν Volcánicos: son flujos de toba aflorando a lo largo del río Challhuahuacho.

ν Paleocanal: subyace a los volcánicos dentro del río Challhuahuacho y los valles bajos del Ferrobamba. Se compone de depósitos aluviales y volcánicos en sus unidades superiores.

Modelo hidrogeológico numérico

Dominio y geometría

El dominio del modelo numérico se definió con dimensiones aproximadas de 26 km (este a oeste) por 22 km (norte a sur) con un área superficial de aproximadamente 367 km². Una vez definida el área de estudio, se creó una malla estructurada de tipo triangular en el software Feflow, incluyendo refinamientos en zonas particulares. La malla del modelo contiene 681,938 nodos, 1’272,880 elementos y consta de 16 capas, las cuales se extienden hasta una profundidad de ~2 km bajo la superficie del terreno (bst) (Figura 5).

Propiedades hidráulicas

FloSolutions (2021) interpretó mediante análisis estadísticos aproximadamente 1,500 pruebas Lugeon, pruebas de inyección / extracción instantánea y pruebas de bombeo realizadas y recolectadas a diferentes profundidades, con el fin de caracterizar la conductividad hidráulica como el almacenamiento especifico. Lo anterior, definió un criterio conceptual que restringió la parametrización del modelo numérico en cinco zonas de profundidad en el Equivalent Porous Media (EPM) anisotrópico definido en el modelo numérico para cada unidad hidrogeológica, fallas y diques.

FloSolutions (2021) identificó que las formaciones ordenadas de mayor a menor permeabilidad, según los datos y conceptos hidrogeológicos, son las siguientes:

Ferrobamba > Paleocanal = Aluvial > 

Glaciofluvial > Skarn > Complejo intrusivo mineralizado > Intrusivo (indiferenciado) > Soraya > Mara

En la Tabla 1 se observan los valores conceptuales de conductividad hidráulica, Ss y Sy para cada unidad hidrogeológica del modelo numérico. Para las fallas y, de acuerdo con el concepto, estas fueron subdivididas y parametrizadas en siete zonas, las cuales tienen parámetros hidráulicos dependientes (en un factor de 5 o 10) de acuerdo con las unidades geológicas que cruzan (Tabla 1), dando así una subdivisión de más de 100 zonas y más de 500 subdivisiones/sectorizaciones para parametrizar el modelo numérico (Figura 6).

Uno de los desafíos al momento de configurar un modelo numérico por parte de un modelador es parametrizarlo y crear selecciones de elementos que serán usadas en la etapa de calibración. Este proceso manualmente tomaría días o semanas dada la complejidad del modelo numérico, pero con este proceso automático la parametrización es en segundos o minutos dada la dimensionalidad del modelo.

Al tener más de 500 subdivisiones y 10 parámetros hidráulicos (i.e., Kxyz, porosidad, Ss, parámetros de Van Genuchten) para la zona de estudio se procedió a implementar un script en Python empleando IFM. Para esto el usuario debe proporcionar como datos de entrada los parámetros y nombres de selecciones futuras en una tabla que tiene asociado un código geológico de una distribución elemental creada por el usuario en Feflow previamente y el modelo numérico. A continuación, se muestra la estructura general simplificada del código:

Table= open(".\Table.dat", "r") # tabla de entrada *.dat

path=".\model.fem" # modelo numérico *.fem

Geo= "Geology" # Nombre distribución elemental doc=ifm.loadDocument(atributos) # carga el archivo *.fem #nombre distribución elemental de usuario

doc.getElementalRefDistrIdByName(atributos)

loop with geology

# para crear selecciones de elementos

doc.createSelection(atributos)

loop with doc.getNumberOfElements():

# Obtiene los códigos de la distribución elemental

doc. doc.getElementalRefDistrValue(atributos)

# asigna elementos a la selección creada

doc.setSelectionItem(atributos)

# Asignación de los parámetros doc.setMatXConductivityValue3D(atributos)# Kx doc.setMatFlowCompressibility(atributos)# Ss (...)

doc.saveDocument()

doc.CloseDocument()

Anisotropía

A partir del modelo conceptual, se han identificado patrones estructurales y deposicionales en la Formación Ferrobamba, así como en las estructuras de falla, que han generado anisotropía, la que varía en orientación y magnitud tanto localmente como en toda el área regional, y está influenciada en parte por la presencia de capas estratigráficas, unidades plegadas, unidades alteradas y zonas de falla.

En general, la conductividad hidráulica en el área de estudio presenta variaciones espaciales, debido a las principales direcciones (x, y, z) que son influenciadas por las unidades plegadas y falladas. Por lo general, la relación de conductividad hidráulica es de 5 a 10, siendo menor en la componente vertical (K_v) en comparación con la componente horizontal (K_h). Sin embargo, en algunos casos, es posible que la conductividad hidráulica vertical (K_v) supere la conductividad hidráulica horizontal (K_h).

Para la zona de estudio, se han definido seis sectores de anisotropía según el concepto donde se encuentran intercalaciones de calizas margosas de la Formación Ferrobamba. Estas intercalaciones pueden aumentar aún más la anisotropía de los acuíferos, especialmente debido a la presencia de diques paralelos a ellas.

Para representar este comportamiento del flujo a través del EPM junto con el cambio de orientación del tensor de conductividad se aplicaron modificaciones de este usando ángulos eulerianos y siguiendo el orden de rotación con base al esquema de Euler “clásico” como una rotación intrínseca. Lo anterior, significa que hay una rotación alrededor del nuevo eje z con un ángulo Φ y una rotación siguiente alrededor del nuevo eje x con ángulo Θ. Estos cambios paramétricos fueron incorporados para la Formación Ferrobamba y las fallas, teniendo en cuenta para estas últimas el valor de referencia del dip y strike, de tal modo que el flujo a través de ellas sea preferentemente a lo largo de estas. Lo anterior, permitirá modelar la anisotropía de manera precisa y realista, considerando la orientación de los planos y líneas de interés en el EPM.

Condiciones de contorno

Se implementaron condiciones de borde asociadas a los componentes naturales del sistema (i.e., recarga por precipitación, descarga por flujo base en quebradas y manantiales) y asociadas a las componentes de infraestructura minera que potencialmente interactúan con el sistema de aguas subterráneas.

La recarga por precipitación proviene de la estimación de la recarga del modelo APLIS incluida en el modelo hidrogeológico conceptual y se aplicó en los elementos de la capa superior del modelo numérico como una condición de borde de flujo (i.e., ingreso/salida en la primera capa del modelo o “In/outflow on top/bottom” en Feflow).

Para la calibración en régimen transitorio, la recarga, además de estar espacialmente distribuida, incluye la variación estacional entre el ciclo anual húmedo y seco través de un factor de estacionalidad y multiplicadores con puntos piloto para ajustar las tasas de salida de las subcuencas de estudio. De acuerdo con el concepto, las áreas con más altas tasas de recarga ocurren en las calizas donde se desarrollan epikarst con un orden de 40% a 70% de la Precipitación Media Anual (PMA). Los cuaternarios poseen tasas de recarga en un orden de 10% a 20% de la PMA y los intrusivos están en el orden de 10% a 30% y las demás unidades con un 3%, respectivamente.

La condición de borde que representa la descarga regional como flujo base de los lagos se implementó en el modelo como condición de borde de carga fija con un nivel de referencia igual a la elevación de la capa 1 del modelo (“Hydraulic Head BC” en Feflow) y para las quebradas, manantiales y bofedales una condición de borde de drenaje libre (“Seepage Face” en Feflow), respectivamente, aplicando una restricción de flujo permitiendo únicamente la descarga del agua subterránea y no recarga.

Para los ríos Record, Ferrobamba, Pamputa, Challhuahuacho, Tambo, Chila, Pumamarca, Conjaca, Jajatuna, Anchapillay y Ceuñama se representaron utilizando la condición de borde de transferencia o de tercer tipo (Cauchy), la que tiene de referencia la carga de elevación igual a la elevación de la capa 1 y las tasas de transferencia de entrada y salida. Por simplicidad, se asumió que la separación de esta capa fue de 3 m para mantener un área constante a lo largo de todos los tramos de río.

Para el río Ferrobamba, dado que su trazado natural fue modificado y reemplazado por un canal de coronación que desvió el cauce por la operación del proyecto desde el 2015, fue necesario inactivar secuencialmente sus tramos de acuerdo con las entradas de las instalaciones mineras (i.e., TSF, depósito de desmonte Huancarane y tajo Ferrobamba) y cambiando el tipo de condición de borde de Cauchy a Neumann y Dirichlet, respectivamente (Figura 7).

En un proceso estándar de modelación no es posible hacer este cambio directamente en la interfaz o código numérico desde la configuración inicial que realiza el usuario, ya que comúnmente se agrega solamente una condición de borde a un nodo durante toda la simulación, por tanto, para esta representación particular de cambio de condición de borde, donde existe superposición, se optó por usar IFM para cambiar el tipo de condición de borde con base al crecimiento del componente minero. Para esto el usuario debe proporcionar como datos de entrada series temporales de condiciones de borde de los componentes mineros, distribuciones nodales temporales por cada componente minero previamente.

A continuación, se muestra la estructura general simplificada del código:

# callback para hacer cálculos antes de cada iteración

def preTimeStep(doc):

loop with doc.getNumberOfNodes()

# identifica configuración de 

crecimiento de la instalación con base a una

# distribución nodal de usuario

if doc.getNodalRefDistrValue

(atributos) > 0:

# contrasta el tiempo de 

simulación con el valor

if doc.getAbsoluteSimulationTime() >=

# remueve previa configuración de contorno

doc.getNodalRefDistrValue(atributos): doc.setBcFlowTypeAndValueAtCurrentTim e(atributos)

# aplica nueva condición de borde

doc.setBcFlowTypeAndValueAtCurrentTim e(atributos)

Para los componentes mineros el tajo Ferrobamba, fue implementado en condiciones transitorias empleando un esquema de desarrollo de taludes de los tajos. Durante la operación de mina del tajo Ferrobamba para la fase de calibración (2015 - 2021) se emplearon los planes de minado anual de la unidad minera y se generaron series temporales de elevación para simular la excavación progresiva.

La excavación se simuló con una combinación de la inactivación de los elementos de la malla y la implementación de condiciones de borde de carga constante, a la cual se le asignó una restricción de flujo para permitir únicamente la descarga de agua subterránea en las paredes de los tajos. Además, se asignó a cada condición de borde una serie de tiempo individual correspondiente a la elevación del LOM o topografía de cada tajo.

Para la recarga del depósito de desmonte se aplicó el concepto de infiltración para los nodos de los elementos de la capa superior. Esta condición de infiltración se representa con la condición de borde tipo 2 o Neumann, en la que se aplicó una serie temporal de tasa de infiltración obtenida a través de simulación 1D de (DHI 2021a), aplicando un factor de reducción de 0.2 (asumiendo una efectividad de 80% de los sistemas de subdrenaje ubicados bajo ellos).

Calibración

El modelo Feflow se calibró para objetivos de niveles y flujos base en subcuencas hidrogeológicas y flujos observados de la operación del tajo Ferrobamba en condición seca para las condiciones de régimen estacionario y transitorio. Se realizó la calibración utilizando el software de estimación de parámetros (PEST) desarrollado por Doherty (2018) disponible en Fepest. Inicialmente, se calibró el modelo estacionario y después de que los niveles simulados en condición preminado eran aceptables con el concepto, datos observados y estadísticos inferiores al 10% de NRMSE. Posteriormente, se continuaba con la calibración transitoria desde el inicio de la operación hasta septiembre de 2021, la cual incluye el desarrollo y la expansión de las instalaciones de la unidad minera usando Fepest inicialmente.

Este proceso repetitivo se realizó de forma paralelizada en Fepest con Beopest de calibración transitoria en una máquina virtual de 48 cores con procesador Intel(R) Xeon(R) Gold 6248R CPU @ 3.00 GHz y 2.99 GHz (dos procesadores) y 256 Gb de RAM. Sin embargo, esta forma de trabajo “clásica” en modelación de aguas subterráneas tomaba alrededor de 60 horas para el modelo de calibración estacionario y una semana para tres iteraciones en promedio para el modelo transitorio, dada la complejidad y miles de zonificaciones creadas para este EPM. Por tanto, se desarrolló una tecnología de acoplamiento de calibración estacionaria - transitoria empleando Python e IFM, lo que permitió minimizar los tiempos de corrida de calibración, optimizar los flujos de trabajo de horas de pre y procesamiento de resultados, minimizar el riesgo humano y aumentar el grado de satisfacción del modelador al reducir su “estrés” en el trabajo diario (Figura 8 y Figura 9).

Como datos de entrada, el usuario debe definir los archivos de control que requiere PEST en su archivo *.pst y hacer las modificaciones pertinentes en los archivos requeridos por PEST (i.e.,*.obs, *.par, *.ins, etc.) para la correcta ejecución por cada modelo (i.e., estacionario y transitorio). Posterior a esto, debe introducir algunos scripts adicionales en los archivos Feflow (*.fem), que permitirán exportar la carga hidráulica del modelo estacionario y actualizarla en el modelo transitorio junto con sus parámetros calibrados.

Seguidamente, se debe modificar los archivos de control de la ejecución en Fepest en formato Bach (*.bat) para que se ejecute inicialmente la simulación estacionaria, exporte los parámetros y carga hidráulica simulada, y posteriormente, los actualice en el modelo de calibración transitorio al finalizar.

Una vez que se ejecute correctamente este procedimiento, el usuario podrá incorporar códigos adicionales que permitan realizar el acople con un método de paralelización disponible en las herramientas de PEST desarrolladas por Doherty (2018) (i.e., Beopest, PEST- HP) y optimizar los tiempos de corrida mediante un acople paralelizado de la calibración. En la Figura 9 se puede observar de forma general la metodología o “paso a paso” de la calibración estacionaria y transitoria y en la Figura 10 una vista general del funcionamiento desde una ventana de comandos con un “worker” asignado para la calibración.

Resultados

Como datos de entrada paramétricos se incorporaron más de 1,500 zonificaciones para áreas de conductividad hidráulica, zonas de transferencia para los ríos, parámetros de conductos kársticos y parámetros de recarga representados mediante multiplicadores y la metodología de puntos piloto (Doherty 2018) en PEST, lo que significa que los valores del parámetro de recarga definido como valor inicial no será reemplazado por el valor del parámetro PEST antes de la simulación, sino que se multiplicará el valor del parámetro actual por el valor del parámetro PEST.

Para los datos de entrada observados, de un universo de 264 piezómetros disponibles para el proyecto, se realizó un proceso de verificación de calidad de la información (tanto de construcción como de datos de niveles) para seleccionar los niveles piezométricos que se emplearon en la calibración acoplada quedando un total de 189 puntos de observación para la calibración estacionaria, y 201 puntos de observación para la calibración transitoria.

El modelo de agua subterránea también consideró durante la calibración estacionaria y transitoria los caudales base promedio de temporada seca asociados con subcuencas hidrogeológicas seleccionadas de interés ambiental por parte del modelo conceptual. Cada subcuenca corresponde a una delimitación hidrogeológica definida en el modelo conceptual y no coincide necesariamente con las subcuencas hidrográficas. Las subcuencas fueron separadas entre subcuencas prioritarias (asociadas directamente con el proyecto) y subcuencas no prioritarias (no asociadas directamente).

El mismo enfoque se realizó con datos observados de manantiales con caudal mayor a 1 L/s, ya que la calibración de flujos base de manantiales con un modelo numérico regional debe ser tratada con cuidado, debido a que generalmente los caudales son producto de condiciones hidrogeológicas locales difíciles de representar en un modelo regional. Como datos de entrada adicionales, para llevar a cabo la calibración del modelo transitorio, se utilizaron datos observados pertinentes relacionados con los flujos de infiltración del tajo para el periodo 2015 - 2021 en la calibración transitoria.

Una vez construidos los archivos *.fps (archivos de Fepest) y *.pst con sus archivos de control de PEST, se procedió a realizar una reponderación de pesos por todas las observaciones corriendo una iteración y la herramienta de la suite de PEST pwtadj (Perturbation Weights Adjustments o ajustes de pesos de perturbación) para darle pesos a los datos observados según su importancia relativa para el ajuste del modelo.

Del acople de la calibración estacionaria - transitoria se empleó la misma máquina virtual para comparar los tiempos de corrida inicial, obteniendo así una reducción del tiempo de optimización total a dos días o aproximadamente un 80% del tiempo total empleado por separado de 10 días. Posterior a esto, se parametrizaron los modelos estacionario y transitorio final con los resultados contenidos en los archivos de salida de PEST (i.e., *.rec y *.par) y se procedió a realizar la evaluación de resultados.

Niveles piezométricos en condición estacionaria

De la calibración estacionaria se obtuvo una calibración preminada para los niveles piezométricos con un residual medio de -13 m, un error medio absoluto de 25 m y un error cuadrático medio normalizado (NRMSE, por su nombre en inglés Normalized Root Mean Square Error) de 5% (Tabla 3). Esta es una estadística de calibración general aceptable para un sistema de flujo de lecho rocoso en un área de diferencias topográficas significativas e inferior al 10% de lo recomendable para guías de modelación internacionales. En la Figura 11 se puede observar la distribución de residuales y el gráfico de dispersión (en inglés, scatter plot), el cual es muy útil para visualizar y analizar la relación entre los valores observados y simulados durante el proceso de calibración.

Flujos base en condición estacionaria

Para los flujos base de las subcuencas hidrogeológicas, se observa un ajuste razonable de los caudales en órdenes de magnitud. Siendo en términos globales un 90% similar para la zona de estudio donde se tenía un total de flujo base observado de 505 L/s respecto a 2,230 L/s simulados en temporada seca para condiciones estacionarias.

Flujos en manantiales en condición estacionaria

Para los flujos de los manantiales, en la Figura 12 se observa un ajuste razonable de los órdenes de magnitud de los flujos de salida de los manantiales priorizados presentes en la Formación Ferrobamba a excepción de los puntos MA-PA-10 y MA-FU-110, donde no se simula caudal. Respecto de los manantiales priorizados que no hacen parte de este grupo con características kársticas, el modelo no logra simular adecuadamente los caudales estimados, siendo en general cerca de un orden de magnitud más alto que los caudales simulados. Lo anterior, posiblemente está correlacionado con los parámetros de recarga estimados desde el concepto o la falta de caracterización hidráulica de las unidades presentes en estas zonas para representar la dinámica piezométrica local del sistema.

Balance hídrico en condición estacionaria

La Figura 13 muestra los componentes del balance hídrico simulado para la calibración estacionaria. Las entradas totales al sistema subterráneo son de 3,203 L/s distribuidos en 3,160 L/s de recarga por precipitación (similar a los 3,218 L/s de recarga por precipitación conceptual), 325 de lagos y 42 L/s de recarga de los ríos. Respecto a las descargas, el modelo simula una descarga de flujo base hacia los ríos de 1,097 L/s, quebradas en 1,572 L/s, valores que están dentro del intervalo promedio conceptual por subcuenca prioritaria observada. Lo anterior, evidencia que el balance hídrico de la calibración estacionaria es adecuado y razonable, dado el conocimiento actual del sistema regional de la zona de estudio.

Niveles piezométricos en condición transitoria

Los resultados de la calibración transitoria se evaluaron para el periodo 2012 (fecha más antigua con datos de nivel piezométrico) hasta 2021, ajustando los parámetros hidráulicos con el fin de lograr un ajuste adecuado de tendencias de niveles piezométricos y caudales del tajo Ferrobamba observados desde el inicio de operaciones en 2015, principalmente.

En términos generales, la calibración transitoria logra representar las condiciones de operación del proyecto para el periodo 2012 - 2021. Por otra parte, existen sectores en los que, dada la baja permeabilidad y recarga del sistema, no permiten representar la estacionalidad de los niveles piezométricos o incluso simular la dinámica de la carga hidráulica que presenta el tajo Ferrobamba ante los fenómenos de compartimentalización, sistema de fallas, diques, intrusiones y demás, que limitan el ajuste local de los niveles piezométricos a escala regional.

Al realizar un cálculo de estadísticos transitorios (Tabla 4) empleando el MAE, este presenta un buen ajuste y poca variación durante los primeros años (2013 a 2015), manteniéndose inferior a los 30 m con residuales medios de -7 m y se incrementa hasta un máximo de 74 m para 2020, con residuales medios -47 m. En general, el modelo Feflow es capaz de simular el comportamiento general del sistema de agua subterránea ante la incertidumbre y limitaciones de información existente para la unidad minera.

A modo de visualización en la Figura 14 se muestran doce hidrogramas de niveles piezométricos observados y simulados para la zona de estudio. En ella, se aprecia que en general los niveles observados tienen en algunos casos ascensos y descensos de niveles significativos, producto ya sea de la estacionalidad o de la propia operación de agua en el tajo, que no son simulados por el modelo numérico. Sin embargo, la tendencia general de niveles en este sector si es correctamente simulada, con una ligera tendencia al descenso del nivel piezométrico hacia el final de la simulación.

Flujos en el tajo Ferrobamba en condición transitoria

Finalmente, en la Figura 14, se presenta un gráfico de barras que representa el ajuste de los caudales observados y simulados durante el periodo de calibración (2015 - 2021). En general, se aprecia un buen ajuste entre los caudales observados y los caudales simulados durante los periodos de época seca. Sin embargo, es importante destacar que no se incluyó el periodo húmedo en el análisis de calibración, debido a la falta de confiabilidad en los datos observados de infiltración. Durante el periodo húmedo, se registraron tiempos de almacenamiento prolongados en piscinas o sumideros antes de ser bombeados a la superficie. Estos registros atípicos dificultan la precisión de los datos y su uso en el proceso de calibración. Por lo tanto, se decidió omitir el periodo húmedo en el análisis y centrarse únicamente en los datos confiables de la época seca para evaluar el ajuste de los caudales observados y simulados.

Conclusiones

1. Los estudios hidrogeológicos desempeñan un papel crucial en la minería a tajo abierto y en particular en los Andes peruanos. Estos permiten comprender los efectos de la actividad minera en los recursos hídricos y diseñar estrategias efectivas de mitigación. La gestión adecuada del agua subterránea se considera fundamental para minimizar los impactos ambientales y garantizar la sostenibilidad de los proyectos mineros.

2. La optimización de la calibración acoplada en los modelos hidrogeológicos es un área de interés en hidrogeología. Pocas experiencias se tienen documentadas y mucho menos el procedimiento metodológico expuesto en el presente artículo. 

3. Con la exploración de nuevas formas de combinar la calibración estacionaria y transitoria, se busca ahorrar tiempo en el proceso de modelado y brindar más oportunidades para evaluar los resultados y analizar la incertidumbre, un ejemplo de esto es que con este proyecto se redujo en aproximadamente un 80% los tiempos de ejecución de la calibración. 

4. La eficiencia en los proyectos mineros es un objetivo importante, y la aplicación de herramientas numéricas avanzadas puede contribuir a lograrlo para optimizar los recursos computacionales y la huella de carbono que emplean, aminorar el riesgo humado por el error en el pre y post procesamiento de resultados para modelos complejos y la optimización del presupuesto para este tipo de proyectos en consultoría.

5. Cada vez es más común la utilización de tecnologías de modelamiento y programación en la hidrogeología minera. Herramientas computacionales como son los softwares Leapfrog y Feflow que, junto con el lenguaje de programación Python, están siendo utilizadas para mejorar la precisión y eficiencia de los modelos hidrogeológicos. Estas tecnologías permiten una mejor estimación de los impactos ambientales y operacionales, así como una toma de decisiones más fundamentada en la gestión de los recursos hídricos.

6. El modelo Feflow fue calibrado en estado estacionario con niveles piezométricos con un NRMS del 5% y simula adecuadamente los flujos base observados a nivel conceptual y/o simulados. Los parámetros calibrados de conductividad hidráulica y almacenamiento están dentro de los intervalos conceptuales excepto para algunas unidades intrusivas de baja permeabilidad, en donde el intervalo conceptual presenta mayor incertidumbre. Para la calibración transitoria los niveles piezométricos y caudales en el tajo Ferrobamba denotaron un ajuste aceptable para el modelo regional (con NRMSE transitorio de niveles que varía entre 2% y 24%), logrando simular en particular y adecuadamente la tendencia de flujos pasivos de descarga del tajo Ferrobamba para los periodos secos.

7. Finalmente, se destaca que para cualquier proyecto hidrogeológico se requiere de un enfoque multidisciplinario. La colaboración entre geólogos, hidrogeólogos y otros expertos resulta esencial para comprender la complejidad geológica e hidrogeológica de los proyectos mineros. La integración de información geológica, estudios de campo y pruebas hidráulicas permite un desarrollo más completo de los modelos geológicos y facilita la toma de decisiones informadas en relación con la gestión del agua subterránea y el correcto análisis de la incertidumbre de los modelos numéricos hidrogeológicos.

Bibliografía

DHI. 2021a. Modelos 1D infiltraciones en Depósito de Desmontes – Bambas Proyecto MEIA 4.

DHI, G. 2021b. Modelación numérica de la expansión del Depósito de Relaves TSF1 de Mina Las Bambas, Perú.

Doherty, J. 2018. Model-independent Parameter Estimation User Manual Part I: PEST, Sensan and Global Optimisers. Watermark Numerical Computing.

FloSolutions. 2021. Estudio Hidrogeológico de Mina Las Bambas para el MEIA4.

Golder, A.I. 2010. Estudio de Impacto Ambiental Proyecto Minero Las Bambas, Anexo H4 Línea Base de Hidrogeología. Lima.

Percy, J. 2018. Influencia de los factores de riesgo que provocan accidentes de trabajo, aplicación de la seguridad basada en el comportamiento y evaluación de resultados, Unidad Minera las Bambas. Arequipa.

Perelló, J., Carlotto, V., Zarate, A., Ramos, P., Posso, H., Neyra, C., Caballero, A., Fuster, N. & Muhr, R. 2003. Porphyry-Style Alteration and Mineralization of the Middle Eocene to Early Oligocene Andahuaylas-Yauri Belt, Cuzco Region, Peru. Economic Geology, 98, 1575-1605, doi: 10.2113/gsecongeo.98.8.1575.

Piteau, A. 2010. Modelo Numérico del Flujo de Agua Subterránea y Evaluación de los efectos del flujo de Agua Subterránea. Lima.

Surran, C. 2022. Las Bambas copper mine expansion approved by Peru regulators (NYSEARCA:COPX). World Wide Web Address: https://seekingalpha.com/news/3817451-las-bambas-copper-mine-expansion-approved-by- peru-regulators.