1. Introducción
El crecimiento demográfico acelerado del Valle de Toluca, evidente desde la década de los 60’s, se ha atribuido a la cercanía con la Ciudad de México, así como al asentamiento de una de las zonas industriales más importantes del país (Sánchez y Orozco, 2015; GEM, 2019). Al margen de la expansión urbana y el aumento en la demanda hídrica, surgió la necesidad de explorar recursos hídricos subterráneos para suministrar agua a una población que anteriormente se abastecía de lagunas (Albores, 1995). Las necesidades de agua para fines domésticos, agrícolas e industriales del Valle de Toluca son cubiertas mediante 1,109 pozos de bombeo. Sumado a esto, se localizan 198 pozos pertenecientes al Sistema Lerma, cuyo fin es extraer agua para la Ciudad de México (Figura 1).
La rápida explotación de agua subterránea desde 1950 ha tenido efectos como la desecación de las lagunas Chignauapan, Tultepec y San Bartolo (Esteller y Díaz, 2002; Esteller y Andreu, 2005; Rudolph et al., 2006; Calderhead et al., 2010a; Hancox et al., 2010; Salas et al., 2011), así como subsidencia regional con tasas de hasta 9 cm/año en zonas de depósitos aluviales y lacustres (Calderhead et al., 2010a, 2010b, 2011, 2012; Castellazzi et al., 2017). De manera complementaria, se han manifestado fracturas en la Zona Metropolitana del Valle de Toluca (ZMVT), afectando edificios, tuberías y vías de comunicación (Murillo, 2008; CICEM, 2008a; CICEM, 2008b; Castellazzi et al., 2017). A fin de analizar tales discontinuidades, se han realizado estudios enfocados a la cartografía y medición por técnicas InSAR (Figueroa, 2004; Almazán, 2017; Castellazzi et al., 2017; López, 2019). Por su parte, el gobierno del estado reporta las fracturas en los Atlas de Riesgos, y coloca carteles en las áreas afectadas a modo de prevención. No obstante, la profundización en el mecanismo que genera tales discontinuidades en el Valle de Toluca continúa siendo un área de investigación en desarrollo.
Diversos estudios apuntan a la disminución de niveles piezométricos, generación de conos de abatimiento y presencia de estratos rígidos con montículos enterrados como promotores de la compactación diferencial, fracturas y reactivación de fallas (Pacheco et al., 2006; Gambolati y Teatini, 2015; Ochoa et al., 2015; Carreón et al., 2016; Frigo et al., 2019; Howard y Zhou, 2019; Schuck et al., 2020; Li et al., 2021; Li et al., 2022). A este respecto, las fallas son planos de rotura en macizos rocosos con movimiento relativo entre los bloques asociadas a procesos tectónicos (González et al., 2002). En comparación, las fracturas consisten en grietas de tensión asociadas a una compactación diferencial que genera componentes en el plano vertical y horizontal (Burbey, 2002, Conway, 2016). De acuerdo con tales definiciones, resulta adecuado caracterizar a las discontinuidades del Valle de Toluca como fracturas.
Este estudio se propuso con el objetivo de evaluar la configuración geológica/geométrica de la cuenca del Valle de Toluca, la evolución piezométrica y su posible relación con la ubicación de fracturas de la ZMVT. En ese marco de investigación, se plantea una revisión y reinterpretación del modelo geológico vigente (Calderhead et al., 2011) con base en nueva información de pozos. Se analizaron las configuraciones y gradientes piezométricos por métodos geoestadísticos durante el periodo de 1975 a 2018. Este estudio constituye el primer acercamiento en la modelación de desplazamientos tridimensionales en el acuífero del Valle de Toluca.
2. Contexto geológico e hidrogeológico
La cuenca del Valle de Toluca con una extensión aproximada de 2,100 km2, se localiza al centro de México, en el marco de la Faja Volcánica Transmexicana (FVTM). Está delimitada al este por la Sierra de las Cruces (3,000 msnm), y al oeste por el volcán Nevado de Toluca (4,560 msnm) (Figura 3). Los materiales geológicos se agrupan en tres clases: rocas ígneas, materiales piroclásticos y depósitos sedimentarios (Figura 2). Las rocas ígneas se constituyen principalmente por andesitas y basaltos conformando el Nevado de Toluca, La Sierra de las Cruces, así como una serie de domos y conos cineríticos localizados en el oeste del valle. Le sobreyacen materiales piroclásticos representados por lahares, bloques, cenizas, flujos de pómez, tobas y brechas volcánicas que afloran al pie de las sierras, que corresponden a la formación Tarango (Macías et al., 1997). Los depósitos aluviales y lacustres que rellenan la cuenca, se interdigitan con materiales volcánicos.
En el ámbito estructural la zona está afectada por sistemas de fallas regionales que atraviesan las estructuras volcánicas mayores. El sistema de fallas San Antonio (SFSA) se localiza entre los volcanes San Antonio y Nevado de Toluca con una orientación general NE-SW. Al sur, el sistema de fallas Tenango (SFT) exhibe una orientación E-W y se caracteriza por fallas curvas y discontinuas, ruptura y formación de zonas de transferencia (García et al., 2000). En el margen este, el sistema de fallas de la Sierra las Cruces presenta tres orientaciones principales N-S, NE-SW y E-W (García et al., 2008).
En la ZMVT se han registrado fracturas en zonas circundantes al cerro La Teresona y las colonias Santa María Totoltepec y San Pedro Totoltepec (Figueroa, 2004; Murillo, 2008; GEM, 2013; Arroyo, 2016; Almazán, 2017; Castellazzi et al., 2017; GEM, 2019; López, 2019). Se distinguen las fracturas Miltepec (NE-SW), fracturas Coatepec (N-S), fracturas El Calvario (NE-SW y NW-SE). Se destaca la fractura Totoltepec con dirección NE-SW, que presenta un desplazamiento horizontal de 80 cm y vertical de hasta 1.40 m (Figura 3).
La organización hidrogeológica se asocia a un sistema multicapa complejo, constituido por dos acuíferos separados por un acuitardo. La unidad acuífera superior de tipo libre se caracteriza por ser un medio poroso compuesto por depósitos clásticos no consolidados e intercalados con flujos piroclásticos y horizontes de pómez. Esta unidad engloba a las formaciones Chalma, Tarango, así como los flujos andesíticos/basálticos y los depósitos aluviales (Figura 2b), alcanzando un espesor aproximado de 200 m (CONAGUA, 2002; Cervantes y Armienta, 2004; Calderhead et al ., 2010a; 2010b). En el acuitardo intermedio se engloban tobas lacustres de baja permeabilidad y arcillas lacustres. El acuífero inferior corresponde a un medio fracturado de tipo semiconfinado, compuesto fundamentalmente por andesitas y basaltos.
La Sierra de las Cruces representa la principal área de recarga debido al alto grado de fracturamiento por enfriamiento de las rocas aflorantes de la cima (Birkle et al., 1998). Asimismo, este flujo se integra con el proveniente del Nevado de Toluca para, a su vez, continuar con una dirección S-N similar al cauce del río Lerma hacia la cuenca de Ixtlahuaca/Atlacomulco.
3. Materiales y métodos
3.1. Integración de información geológica e hidrogeológica
Los datos utilizados en este trabajo comprenden un análisis de la geología superficial y modelos de elevación digital recopilados de fuentes gubernamentales de acceso libre. En tanto que, la estratigrafía y piezometría histórica se obtuvieron de estudios previamente realizados en el acuífero de Toluca (Tabla 1). De igual modo, se integra información reciente de niveles estáticos proporcionada por la Comisión Nacional del Agua (CONAGUA, 2018).
3.2. Conceptualización del medio geológico mediante modelación tridimensional
3.2.1. Clasificación de materiales geológicos
En el Valle de Toluca, las fracturas se han atribuido a la excesiva extracción de agua y a la variación en la compresibilidad del medio poroso entre macizos rocosos y materiales no consolidados (CICEM, 2008a). Para su análisis, resulta necesaria una agrupación de materiales con características hidrogeológicas y propiedades de deformación del medio poroso similares. Los lineamientos utilizados para la agrupación de litologías son: (1) Procedencia del material: origen volcánico o sedimentario, (2) Tamaño de grano de materiales sedimentarios: granos gruesos y finos, y (3) Contenido de arcilla. Se evaluaron 260 columnas litológicas correspondientes a sondeos profundos, sondeos de penetración estándar y pozos a cielo abierto.
3.2.2. Interpolación de horizontes
En ausencia de fallas, pliegues u otras estructuras complejas en el Valle de Toluca, resultan útiles los métodos numéricos explícitos cuyo objetivo es establecer la elevación de contactos geológicos en una determinada ubicación (Wellmam y Caumon, 2018). A este respecto, se aplicó la interpolación por distancia inversa ponderada (Ecuación 1). En el método se establecen pesos de ponderación con base en la distancia entre el punto de valor conocido y el punto de estimación, asumiendo que los puntos cercanos guardan mayor relación.
Donde n es el número de puntos usados, fi es el valor de la variable en el punto analizado y wi es el peso ponderado asignado a cada punto. Para su aplicación se asignaron números de identificación de contacto litológico a cada columna con base en la elevación topográfica de cada unidad y la concordancia con estratos adyacentes. En las zonas carentes de perfiles litológicos se añadieron 8 “Dummy points”, los cuales consisten en pozos hipotéticos y cuyo fin es facilitar la interpolación de estratos (Pawlowsky et al., 1993). La estratigrafía de los Dummy points se basa en la reportada en pozos cercanos. Se elaboraron 116 secciones transversales para establecer una condición inicial de geometría en la interpolación 3d.
El dominio del modelo se estableció con base en la geología del Valle de Toluca, tomando como referencia el contacto entre el acuífero granular y el acuífero fracturado (Figura 4a). Lo anterior se basa en el hecho que la unidad granular es donde reside la mayor cantidad de pozos de extracción, asimismo, en la región montañosa no se cuentan con información suficiente de estratigrafía, piezometría y propiedades hidráulicas.
3.2.3. Asignación de la conductividad hidráulica
En la asignación de la conductividad hidráulica de arcillas y gravas, se retomaron valores obtenidos por métodos de tamaño de grano, pruebas de bombeo, pruebas de recuperación (slug test) y ensayos de permeabilidad Lefranc, cuya información fue sintetizada en diversos estudios desarrollados en el Valle de Toluca. Asimismo, se contemplaron valores propuestos en la literatura para suelos y depósitos volcánicos (Juárez y Rodríguez, 2005; Baumann et al., 2019), así como los propuestos por Calderhead et al. (2011).
3.3. Geoestadística aplicada a la evolución del nivel estático
El nivel de agua subterránea es considerado una variable regionalizada, siendo los métodos geoestadísticos los más utilizados para su interpolación (Kitanidis, 1997; Chilès y Delfiner, 1999; Desbarats et al., 2002; Peterson y Barnett, 2004; Lu et al., 2009; Nikroo et al., 2010). La ventaja de la implementación de tal enfoque, en general, es el análisis de la continuidad espacial de la variable utilizando funciones teóricas que parten del variograma (Ecuación 2).
donde: n es el número de pares experimentales separados por un vector h, Z(xi) es el valor de la variable de interés en el punto xi con i=1…n. Se utilizaron los niveles estáticos medidos en 42 multipiezómetros a una profundidad de sonda de 150 m durante el periodo de 1975 a 2018 (CONAGUA, 2018). Las técnicas de exploración de datos aplicadas comprenden histogramas de frecuencia y variogramas direccionales a 0°, 45°, 90° y 135°, cuyo fin es determinar tendencias de anisotropía.
Las elevaciones de nivel estático fueron simuladas aplicando varios variogramas teóricos, cuyos resultados se evaluaron por validación cruzada con el método Leave one out. En esta validación, se omite cada uno de los datos del set y su valor se estima por kriging ordinario puntual haciendo uso de la función teórica propuesta (Oliver y Webster, 2014). La selección del modelo teórico se basó en los resultados de los índices de ajuste (Ecuación 3) que, por lo general, evalúan qué tanto se subestima o sobreestiman los datos. Los índices utilizados son el error medio (ME), la raíz cuadrada del error cuadrático medio (RMSE) y la tasa de desviación media cuadrada (MSDR), cuyos valores esperados son ME= 0, RMSE= 0, MSDR= 1.
donde, n es el número de datos, Pi es el valor simulado, Oi es el valor observado y σk2 es la varianza de kriging. Una vez establecida la función matemática con los menores errores, se realizó la interpolación de nivel por kriging ordinario (Ecuación 4) a través del software Surfer v. 13 (Golden Software Inc., 2015).
donde
3.4. Interpretación de fracturas en el contexto geológico y piezométrico
La construcción de secciones transversales en las zonas de principales fracturas permitió delimitar la geología subsuperficial y, por ende, evaluar la configuración de las rocas y los estratos sedimentarios. La integración del modelo geológico tridimensional, evolución piezométrica y ubicación de conos de abatimiento por extracción de agua permite formular una conceptualización del mecanismo que genera las fracturas.
4. Resultados y discusión
4.1. Modelo geológico tridimensional
Debido a que los macizos rocosos se encuentran mayormente a grandes profundidades, las variaciones piezométricas son pocos significativas. Se agruparon los basaltos, andesitas y riolitas en el grupo de rocas volcánicas. La clasificación de depósitos piroclásticos parte de su baja densidad y alta porosidad, que, a su vez, se caracteriza por lahares, bloques, cenizas, flujos de pómez, tobas blandas poco cementadas y brechas volcánicas. En el caso de los suelos, la categorización considera el tamaño de grano como punto de partida, distinguiéndose gravas, arenas, limos y arcillas. Las gravas y arenas se agruparon en una sola clase. Los estratos con contenido de finos, de acuerdo con la clasificación SUCS (Unified Soil Classification System) reportada en los perfiles estratigráficos se dividieron en granos gruesos embebidos en matriz fina y arcillas puras. De esta manera, es posible establecer una diferencia en el grado de consolidación de los estratos y potencial factor de subsidencia diferencial. Los 6 tipos de materiales aplicados a la reclasificación de todos los pozos son 1) Arcillas, 2) Gravas y arenas, 3) Granos gruesos en matriz fina, 4) Conglomerado, 5) Depósitos piroclásticos, 6) Sólidos volcánicos (Figura 4b).
La interpolación de estratos partió de las secciones transversales construidas (Figura 4c), obteniendo un modelo geológico (Figura 4d) conformado por 14 capas (Tabla 2). De manera general, se identifica un estrato de rocas volcánicas como basaltos y andesitas, sobreyacidos por depósitos piroclásticos. Sobre esta secuencia, se identifican conglomerados, gravas y arenas, conforme la distancia a la superficie es menor, la granulometría disminuye alcanzando lechos arcillosos. Se identifican dos domos volcánicos que subyacen al centro de la ciudad de Toluca y a la colonia de San pedro Totoltepec respectivamente. Este último coincide con la zona estable ante subsidencia caracterizada por Arroyo (2016) y denominada Punto de Control en el estudio de subsidencia vertical realizado por Calderhead et al. (2010a).
Material | Clave | Descripción | Espesor promedio (m) | Rango de conductividad hidráulica K (m/s) |
---|---|---|---|---|
Capa 14 | C 14 | Arcillas | 11.36 | 1.0E-11 - 1.0E-07 |
Capa 13 | GS 13 | Gravas y arenas | 4.56 | 1.0E-04 - 1.0E-01 |
Capa 12 | GF 12 | Granos gruesos en matriz fina | 8.96 | 1.0E-07 - 1.0E-04 |
Capa 11 | FP 11 | Flujos piroclásticos | 5.68 | 1.0E-06 - 1.0E-02 |
Capa 10 | AB 10 | Andesitas y basaltos fracturados | 0.90 | 8.0E-09 - 3.0E-04 |
Capa 9 | C 9 | Arcillas | 0.92 | 1.0E-11 - 1.0E-07 |
Capa 8 | GS 8 | Gravas y arenas | 33.04 | 1.0E-04 - 1.0E-01 |
Capa 7 | CO 7 | Conglomerado | 11.61 | 1.0E-05 - 1.0E-03 |
Capa 6 | GF 6 | Granos gruesos en matriz fina | 31.96 | 1.0E-07 - 1.0E-04 |
Capa 5 | GS 5 | Gravas y arenas | 28.72 | 1.0E-04 - 1.0E-01 |
Capa 4 | CO 4 | Conglomerado | 3.44 | 1.0E-05 - 1.0E-03 |
Capa 3 | FP 3 | Flujos piroclásticos | 27.35 | 1.0E-06 - 1.0E-02 |
Capa 2 | GS 2 | Gravas y arenas | 8.58 | 1.0E-04 - 1.0E-01 |
Capa 1 | AB 1 | Andesitas y basaltos fracturados | 138.81 | 8.0E-09 - 3.0E-04 |
4.2. Evolución piezométrica
4.2.1. Ajuste del variograma teórico
Debido a los diferenciales topográficos de hasta 200 metros en la cuenca, en los histogramas de frecuencia se identificaron mediciones de 2,750 - 2,820 msnm que se alejan de la media (2,510 - 2,620 msnm) y, en consecuencia, aumentaban la desviación estándar del set. En tal sentido, los piezómetros de las serranías fueron descartados, y con las estaciones restantes, se elaboraron variogramas direccionales a 0°, 45°, 90° y 135° (siendo 0° la dirección horizontal y 90° la dirección vertical). Las mayores semivarianzas se presentaron a 0° y 45°, mientras que, las configuraciones más estables se observaron a 90° y 135°, destacando la característica anisotrópica de la carga hidráulica, siendo los ejes principales de anisotropía determinados por la orientación del flujo de agua subterránea (Kitanidis, 1997). A este respecto, las funciones teóricas que se ajustaron mejor son las que incluyen un coeficiente de anisotropía comprendido entre 1.804 y 2 para un ángulo comprendido entre 118° y 135° (Tabla 3).
Parámetro | 1975 | 1984 | 1994 | 2004 | 2015 | 2018 |
---|---|---|---|---|---|---|
Datos | 32 | 30 | 28 | 30 | 28 | 27 |
Función de variograma | Lineal | Lineal | Exp | Lineal | Lineal | Lineal |
Pendiente | 0.03 | 0.035 | NA | 0.03 | ||
Meseta | NA | NA | 925.2 | NA | NA | NA |
Rango | NA | NA | 16000 | NA | NA | NA |
Factor de anisotropía | 1.804 | 2 | 2 | 1.804 | 1.804 | 1.804 |
Ángulo de anisotropía | 135 | 135 | 120.4 | 135 | 135 | 135 |
ME (m) | -0.6930 | -0.6344 | 0.1425 | -0.3599 | -0.1944 | -0.1676 |
RMSE (m) | 8.8030 | 9.3824 | 8.5307 | 8.3495 | 9.0941 | 8.7082 |
MSDR | 0.9748 | 0.9711 | 0.9646 | 0.9685 | 0.9647 | 0.9633 |
4.2.2. Evaluación del comportamiento evolutivo del nivel estático
Las superficies del nivel estático muestran un flujo del agua subterránea con dirección preferencial SE-NW, siendo las áreas de recarga la Sierra de la Cruces y el volcán Nevado de Toluca. De esta manera, las isopiezas de los años 1975 y 1984 indican un patrón regional del flujo subterráneo proveniente de las áreas de recarga hacia el cauce del río Lerma, que simula al enrutamiento de agua superficial (Figuras 5a y 5b). El análisis del año 1994 presenta un incipiente cono de abatimiento, así como la deformación de las líneas equipotenciales (Figura 5c). En las décadas posteriores, se observa una disminución de la distancia relativa entre isolíneas y, por lo general, un retroceso hacia las zonas altas, lo que sugiere un aumento del gradiente hidráulico en el centro de Toluca (Figura 5d, 5e y 5f). La existencia de un cono de abatimiento en el análisis del año 2015 denota una deformación local de las equipotenciales, comportamiento que incluso se mantiene en el año 2018.
Las estaciones con los residuales máximos y mínimos fueron analizadas para cada año de krigeado. En general, de 1975 al 2004, las estimaciones más inexactas se localizan en las zonas altas cercanas al Nevado de Toluca y, por lo general, a la salida de la cuenca. En contraste, las mejores predicciones se obtuvieron en las áreas de planicie aledañas al cauce del río Lerma. En los últimos análisis, sin embargo, los mayores residuales se ubican en el corredor industrial (Figura 5e) y el centro de Toluca (Figura 5f). Esta configuración de residuales indicaría que la tendencia del flujo subterráneo ha sido influenciada por perturbaciones locales, particularmente la extracción de agua subterránea. Por consiguiente, la aplicación de variogramas teóricos regionales para el periodo entre 2015 y el 2018 es poco recomendada. Los mayores descensos piezométricos acumulados se ubican en las zonas cercanas a las fracturas, reafirmando la relación existente con el bombeo de agua subterránea (Figura 6). Particularmente, en la fractura Totoltepec, el foco de extracción localizado al sur de esta coincide con los estratos de sedimentos no consolidados y el bloque que se compacta.
4.3. Integración de información geológica e hidrogeológica
Con base en las secciones trasversales derivades del modelo geológico 3d, el volcán andesítico que subyace al centro de Toluca funciona como una zona mecánicamente estable en comparación con los sedimentos de grano fino que se depositan en sus alrededores (Figura 7a). Lo anterior explicaría la ausencia de fracturas en el centro de la ciudad, así como la existencia de éstas con múltiples orientaciones en el área circundante. Curiosamente, la fractura Totoltepec se manifiesta en zonas cercanas a un segundo domo de rocas volcánicas. De acuerdo con los registros litológicos de pozos cercanos, se registra una secuencia de basalto andesítico en el pozo 15W, la cual, queda cubierta por un estrato de gravas y arenas (Figura 7b). En este contexto, la fractura Totoltepec, ubicada al sur de dicha estructura, coincide con una secuencia de depósitos piroclásticos, conglomerados, granos en matriz arcillosa y depósitos de grano grueso (en orden de localización del depósito). De tal columna se destacan los estratos con alto contenido arcilloso por su susceptibilidad a procesos de consolidación (Gambolati y Teatini, 2015). La zona de transición roca-sedimento podría suponer un área de esfuerzos cortantes en la cresta de los macizos rocosos generando subsidencia diferencial y fracturas (Figura 8).
5. Conclusiones
Este estudio se propuso con el objetivo de evaluar la configuración geológica/geométrica de la cuenca del Valle de Toluca, la evolución piezométrica y su posible relación con la ubicación de fracturas de la ZMVT. La reinterpretación e interpolación de los registros estratigráficos denota una geometría de estratos altamente irregular. La ubicación de fracturas muestra correlación con la presencia de domos volcánicos y descensos piezométricos. Las discontinuidades en los alrededores del centro de Toluca concuerdan con el domo andesítico subyacente reportado en estudios anteriores (Calderhead et al., 2010a). Por otro lado, el bloque que baja en la fractura Totoltepec coincide con una secuencia de materiales piroclásticos y sedimentarios, donde destacan estratos de granos gruesos en matriz fina de hasta 70 m de espesor. Con base en las clasificaciones derivadas de Sondeos de Penetración Estándar, así como la estratigrafía reportada en los multipiezómetros y pozos del Sistema Lerma, las matrices finas son conformadas mayormente por arcillas que por limos. Lo cual indicaría que el paquete de granos gruesos en matriz fina presentaría un grado de consolidación ante diferenciales piezométricos.
Los resultados de la evolución piezométrica muestran la disminución en los niveles estáticos de manera general en la cuenca, así como un cambio en las direcciones de flujo a partir del análisis del año 1994. Del análisis geoestadístico aplicado a los niveles estáticos históricos, las principales inexactitudes obtenidas en el krigeado se ubican en las zonas de mayor abatimiento. Lo anterior tiene sentido ya que tales perturbaciones locales son difícilmente modeladas por variogramas regionales. A fin de contemplar el escarpado relieve de la cuenca, como una alternativa se sugiere la aplicación de cokriging con una deriva topográfica. Tomando juntos estos resultados, la zona de transición roca-sedimento no consolidado, la presencia de materiales susceptibles a consolidación, además de los descensos piezométricos por extracción de agua, concuerdan con la ubicación de las fracturas reportadas. El conocimiento de los factores que intervienen en el fracturamiento, sin duda, constituye uno de los primeros pasos en la delimitación de zonas de riesgo.