Área de estudio y parámetros climáticos

El área de estudio se localiza al este de los Andes colombianos entre los departamentos de Boyacá y Arauca con una extensión de ~47565 ha (Figura 1). Geológicamente, la SRC se encuentra ubicada en la Cuenca Cordillera Oriental, la cual esta divida por el alto de Santander en dos subcuencas principales, Tablazo Magdalena al occidente y Cocuy al este (^{Etayo et al., 1969}; ^{Cooper et al., 1995}). La SRC localizada en la Subcuenca del Cocuy alcanza altitudes de hasta 5345 msnm y según ^{Rekowsky et al. (2018)} es el glaciar más extenso de Colombia. Las condiciones topográficas de la SRC permitieron una mayor acumulación de hielo principalmente en la vertiente occidental (^{IDEAM, 2001}; ^{Herrera y Ruiz, 2008}; ^{Rekowsky et al., 2018}).

Figura 1 Localización del área de estudio en imágenes satelitales adquiridas de Sentinel-2 (cuadrícula mundial en gránulos o tiles de 100×100 km² por unidad) y Landsat-8 (cuadrícula mundial de referencia WRS-2). 

En esta región, predominan paisajes montañosos con vastas áreas de roca expuesta dominantemente de edad cretácica, las cuales corresponden con depósitos marinos y transicionales formando un espeso cordón sedimentario de aproximadamente 7500 m, sobre los cuales descansan abundantes depósitos cuaternarios de tipo aluvial y glaciar (Figura 2) (^{Cooper et al., 1995}).

Figura 2 A. Mapa geológico (modificado de ^{Fabre et al., 1984}, ¹⁹⁸⁵). B. Estratos buzantes de edad cretácica. C. Masa glaciar de la Sierra Nevada del Cocuy (SRC). 

Esta zona se caracteriza por presentar un régimen de precipitación monomodal en la vertiente oriental y un régimen bimodal en la vertiente occidental (^{Rodríguez, 2015}). Dicho fenómeno es el resultado del doble movimiento de la ZCIT, conocida por ser una región de baja presión, la cual controla la dinámica de los vientos alisios provenientes del noreste y sureste sobre el océano y sobre las masas de tierra en áreas cercanas al Ecuador, generando mayor nubosidad y precipitación (^{Waliser y Jiang, 2015}; Yan, 2015). Este es el caso de la SRC en los Andes colombianos (^{Rodríguez, 2015}; ^{Rekowsky et al., 2018}).

Los valores de precipitación para la estación de la Sierra Nevada del Cocuy reportados por el Instituto de Hidrología, Meteorología y Estudios Ambientales de Colombia (IDEAM) indican un régimen bimodal y valores medios mensuales de precipitación máxima en 24 horas, entre 25 mm y 65 mm (Figura 3A). El análisis de este parámetro permite determinar los meses con condiciones apropiadas para la obtención de las imágenes, además de que posibilita tener un conocimiento base de la relación entre la precipitación y la cobertura de nieve. De acuerdo con ^{Gutzler y Preston (1997)} cuando la extensión espacial y la duración de la capa de nieve invernal son altas, la cantidad de precipitación en verano es menor y viceversa. Los valores registrados de temperatura media máxima en 24 horas oscilan entre 11,42 y 13,51°C para los diferentes meses del año (Figura 3B). Las variaciones de precipitación y temperatura dependen en gran medida de la época del año, localización y altura (^{Rodríguez, 2015}), presentando una relación inversa. Además, dichos valores pueden exhibir fluctuaciones dependiendo de las estaciones meteorológicas analizadas en consecuencia de la dinámica climática local.

Figura 3 Valores máximos de precipitación y temperatura en 24 horas (1980-2020), estación Sierra Nevada del Cocuy. A. Precipitación máxima en 24 horas media mensual. B. Temperatura máxima diaria mensual en 24 horas. 

Debido al cambio climático, durante las últimas 4 décadas se ha observado una disminución de las masas glaciares en los Andes colombianos (^{Guillen et al., 2004}; ^{Ceballos et al., 2006}; ^{Morris et al., 2006}; ^{Herrera y Ruiz, 2008}; ^{Bautista y Correa, 2018}; ^{Rekowsky et al., 2018}; ^{Quintero, 2018}), lo cual resulta de vital importancia, ya que en estos glaciares nacen numerosas fuentes hídricas. La SRC presenta una divisoria de aguas en dirección S-N con dos vertientes principales, la occidental hacia el río Magdalena y la oriental hacia el río Orinoco (^{IDEAM, 2019}).

La disposición topográfica de la SRC, así como los factores climáticos generan un aporte constante del recurso hídrico principalmente a la vertiente occidental hacia los ríos Lagunillas, Cóncavo, Corralitos, San Pablín y Cardenalillo, los cuales fluyen en dirección al río Nevado que posteriormente desemboca en el río Chicamocha (^{IDEAM, 2019}). En la vertiente oriental nace el río Mortiñal y las quebradas Patio Bolos y las Tapias que desembocan en el río Casanare y los ríos Ratoncito y Cubogón los cuales fluyen hacia el río Arauca (^{IDEAM, 2019}). Las aguas producto de la fusión glaciar de las cuales se benefician algunos pobladores rurales de las zonas más altas, tributan hacia los páramos y estos a su vez contribuyen a la provisión de agua de los principales municipios de la región (e.g., Güicán, Cocuy, Panqueba, etc.) (^{IDEAM, 2019}).

Conjunto de datos de sensores remotos

Fueron clasificadas 11 imágenes ópticas las cuales se tomaron preferencialmente en la época seca de la SRC para evitar la ocurrencia de nubes. Una escena del sensor TM del satélite Landsat-4 para el año de 1987, 3 escenas de los sensores TM y MSS del satélite Landsat-5 para los años de 1991, 1997 y 2009, 2 escenas del sensor ETM+ del satélite Landsat-7 para los años de 2000 y 2003, 3 escenas de los sensores OLI+TIRS del satélite Landsat-8 para los años de 2014, 2016 y 2017, y 2 escenas del sensor MSI del satélite Sentinel-2 de la constelación COPERNICUS para los años de 2019 y 2021. Estas imágenes fueron descargadas usando la herramienta GloVis del Servicio Geológico Americano (USGS), y la base de datos de Google Earth Engine.

La aplicación de técnicas de uso conjunto de datos provenientes de Landsat-8 y Sentinel-2, conceden escenas con mayor resolución temporal debido a que las áreas de estudio de cobertura glaciar son ampliamente caracterizadas por su alta nubosidad, obteniendo así imágenes de mejor calidad (^{Li y Roy, 2017}; ^{Singh et al., 2020}; ^{Andreassen et al., 2021}). Además, al ser datos de libre acceso ofrecen la posibilidad de ser un medio eficaz para llevar a cabo análisis de áreas y coberturas (^{Wang et al., 2018}). Una de las mayores ventajas del conjunto de datos Landsat 4, 5 y 7 son sus imágenes históricas, al cual se acredita el registro más extenso de imágenes satelitales (^{Wulder et al., 2019}), permitiendo realizar un análisis multitemporal de la cobertura glaciar en la SRC.

Adicionalmente, fueron utilizadas 5 escenas ópticas de alta resolución como datos de referencia: 3 escenas de ^{Google Earth (2020)}, para los años de 2009, 2014 y 2017 de entre 1,0 a 3,0 m de resolución espacial dependiendo del mosaico; y 2 escenas de “Planet’s high-resolution, analysis-ready mosaics of the world’s tropics” provistas por la Iniciativa Internacional sobre Clima y Bosques de Noruega (NICFI) a través de la plataforma ^{Planet (2020)} de 4,77 m de resolución espacial y que corresponden a los años de 2016 y 2019. Durante enero del 2021 se llevó a cabo una salida de campo que permitió tomar puntos de control con el fin de corroborar las clases propuestas en la clasificación de las imágenes (área glaciar, sueloroca, vegetación y agua), empleando un sistema de posicionamiento global GPS con una precisión de 3 m para una localización confiable en tiempo real. A lo largo de aproximadamente 93 km de caminos en el área de estudio, se recopilaron 500 puntos de control del suelo (Ground Control Point-GCP) para validar la clasificación de la imagen Sentinel-2 MSI de 2021. La fase de campo estuvo restricta a los caminos permitidos al público y a la disponibilidad de guías locales para el acceso al Parque Nacional Natural Sierra Nevada El Cocuy.

Procesamiento digital de las imágenes

Cada una de las escenas procesadas, Landsat-4 (1987), Landsat-5 (1991, 1997, 2009), Landsat-7 (2000, 2003), Landsat-8 (2014, 2016, 2017), y Sentinel-2 (2019, 2021) se convirtieron de valores de número digital (DN) a reflectancia del suelo (GR). La conversión de las escenas se realizó en el módulo de corrección atmosférica (ATCOR) del Software PCI Geomatics 2016. Posteriormente, las imágenes de 1987 y 1997 fueron ortorrectificadas. Todas las imágenes de los diferentes sensores se remuestrearon a un tamaño de píxel de 30*30 m y se recortaron en función del área de estudio. Además, se crearon máscaras de nubes para las imágenes de los años 2003 y 2009 (únicas escenas que presentaron nubes) mediante mapeo manual, dibujando polígonos alrededor de las nubes evidenciadas en la composición RGB 7-4-3 para Landsat 5 y 7 (Tabla 1). En esta composición RGB se incluye una banda del infrarrojo de onda corta en el canal rojo. Para la nieve se debe tener en cuenta su alta absorción en el infrarrojo de onda corta (^{Domine et al., 2006}), por lo cual en esta composición RGB 7-4-3, la cobertura de nieve se visualizará en color cian, mientras que las nubes se verán blancas. Utilizando estas máscaras se extrajeron los píxeles correspondientes a nubes, siendo excluidos durante la clasificación de las imágenes.

Tabla 1 Bandas usadas durante la clasificación para cada sensor, las unidades de longitud de onda se encuentran en μm. 

El índice NDSI propuesto por ^{Dozier (1989)}, relaciona una banda del infrarrojo de onda corta y una del espectro visible. Para este estudio el índice NDSI fue calculado en todas las imágenes analizadas (Ecuación 1) (Figura 4). Las bandas escogidas para el cálculo de este índice en cada sensor fueron seleccionadas en función de su longitud de onda (Tabla 1), siendo consideradas como verde (Green) las bandas con longitudes de onda entre 0,5 y 0,6 μm, y como infrarrojo de onda corta (SWIR) las bandas con longitudes de onda entre 1,55 y 1,75 μm.

(1)

Figura 4 Distribución espacial del Normalized Difference Snow Index (NDSI) para las escenas de 1987, 1991, 1997, 2000, 2003, 2009, 2014, 2016, 2017, 2019 y 2021. 

Clasificación supervisada orientada a píxel

Se definieron 4 clases principales dentro del área de estudio (área glaciar, suelo-roca, vegetación y agua). Para la extracción estadística de las clases espectrales, se recolectaron al menos 25 muestras de polígonos de entrenamiento para cada clase (~13341 píxeles para el área glaciar, ~48442 píxeles para suelo-roca, ~1325 píxeles para vegetación y ~1621 píxeles para agua). Las muestras de entrenamiento se definieron por clase, basados en los GCP colectados en campo, dibujando polígonos alrededor de los puntos mediante mapeo manual y considerando las variaciones correspondientes a cada año. Estos datos se complementaron con las imágenes de alta resolución. A pesar de la diferencia de hasta 40 años entre las imágenes y la adquisición de datos de campo, todas las descripciones de campo georreferenciadas y fotografiadas en el 2021 coincidieron en su mayoría con las clases observadas en las imágenes de años anteriores.

El dibujo de estos polígonos de muestreo se realizó a partir de una composición color RGB 5-4-3 para los sensores ETM, TM y ETM+, RGB 6-5-4 para el sensor OLI y RGB 11-8-4 para el sensor MSI. Posteriormente, se definieron las bandas que fueron usadas para la clasificación de cada escena, incluyendo valores de reflectancia de la región del visible y del infrarrojo cercano (Tabla 1), así como el índice NDSI (^{Metsämäki et al., 2015}). Se empleó el modelo de clasificación de máxima verosimilitud, el cual mostró resultados más acordes con los datos de referencia. Este algoritmo clasifica las imágenes utilizando un proceso de cálculo interactivo, garantizando la máxima diferenciación estadística posible basada en los datos espectrales (^{Rodrigues y Souza-Filho, 2011}).

Evaluación de la exactitud de la clasificación

Para evaluar la exactitud de la clasificación supervisada, se adoptó un modelo de matriz de confusión normalizada y estadística Kappa. Un total de 500 GCP fueron utilizados para evaluar la exactitud de la clasificación de la imagen de 2021 (dando un enfoque principal a la cobertura del área glaciar) (Figura 5). Para evaluar la exactitud de las clasificaciones de los años 2019, 2017, 2016, 2014 y 2009 se emplearon las imágenes de referencia de Google Earth y Planet con 500 puntos de validación localizados de manera aleatoria y estratificados de acuerdo al área de cobertura de cada clase con respecto a la escena total (~250 para el área glaciar, ~150 para el suelo-roca, ~30 vegetación y ~70 de agua). Debido a que las clasificaciones de las imágenes más antiguas no tenían disponibles datos de referencia (1987, 1991, 1997, 2000, y 2003), se usaron las mismas imágenes de Landsat 4, 5 y 7 en el proceso de evaluación de exactitud, generando 500 puntos de validación distribuidos aleatoriamente por toda la escena (Figura 5).

Figura 5 Clases definidas en el análisis de imágenes para la Sierra Nevada del Cocuy (SRC) y número de puntos para la validación de la exactitud. 

Para cada imagen de Landsat y Sentinel-2 fueron calculadas la exactitud de usuario y productor, exactitud global e índice Kappa individual con el fin de definir la concordancia de los datos (^{Congalton y Green, 2008}; ^{Santos et al., 2020}). El uso del índice Kappa en este estudio tiene el fin de enriquecer la evaluación de exactitud de las clasificaciones. Aunque este índice sea discutido (^{Pontius y Millones, 2011}), varias investigaciones actuales emplean el índice Kappa como referencia para la validación del mapeo de áreas con cobertura glaciar (e.g., ^{López-Moreno et al., 2020}; ^{Sood et al., 2021}; ^{Zhang et al., 2020}).

Datos Climáticos

Se evaluaron los valores medios anuales de precipitación y temperatura accedidos a través de la plataforma del Atlas climatológico de Colombia del IDEAM. Esta plataforma ofrece datos medidos reportados por una o dos estaciones meteorológicas por departamento, desde 1981 hasta el 2010, y la posibilidad de acceder a servicios de mapas interactivos. A partir de dichos valores se estableció de manera general la dinámica de los parámetros meteorológicos para el área de estudio y su repercusión en las tasas de retroceso glaciar para la SRC.

Clasificación general y evaluación de la exactitud de mapas multitemporales

La clasificación supervisada para las imágenes de Landsat y Sentinel-2, fue basada en 4 clases (área glaciar, suelo-roca, vegetación y agua), que contienen características distintivas, utilizando los beneficios de cada tipo de datos. Esto permite aprovechar la complementariedad de los datos espectrales para extraer más información y permitir una identificación más clara entre las diferentes clases (^{Ducrot et al., 2010}). En la Figura 6 se muestran las clases definidas en la zona de estudio y los valores correspondientes de área (ha) durante los años 1987, 1991, 1997, 2000, 2003, 2009, 2014, 2016, 2017, 2019 y 2021. El retroceso glaciar se hace más evidente sobre la vertiente occidental de la SRC (^{Herrera y Ruiz, 2008}), donde progresivamente la clase correspondiente al área glaciar cambia a la clase suelo-roca. Las otras clases no presentan una relación de cambio evidente entre los diferentes periodos de tiempo.

Figura 6 Imágenes clasificadas de Landsat y Sentinel-2 que muestran las variaciones temporales y espaciales para la Sierra Nevada del Cocuy (SRC) en los años 1987, 1991, 1997, 2000, 2003, 2009, 2014, 2016, 2017, 2019 y 2021. 

Las exactitudes globales de todas las escenas clasificadas son superiores al 86% (Tabla 2), mientras que los índices Kappa son superiores a 0,76. Generalmente, el índice Kappa para cada clase es superior a 0,7 excepto para la vegetación en las imágenes del 2009 y el 2017, y para la clase de suelo-roca en el 2009 (Figura 7A). Estos valores indican que el desacuerdo general es muy bajo y que las clases presentan un alto nivel de concordancia (^{Souza-Filho et al., 2018}). La exactitud del productor es normalmente superior al 80%, excepto para las clases de agua y vegetación en algunas de las imágenes (1987, 2000, 2009, 2016 y 2017) (Figura 7B). A pesar de esto, los valores para el área glaciar, foco de este estudio, estuvieron siempre por encima de 96,4%. En cuanto a la exactitud del usuario, esta se mantuvo por encima del 80% con excepción de las coberturas de suelo-roca, agua y vegetación para las imágenes del 2009, 2016 y 2017 (Figura 7C). No obstante, dicha exactitud para la clase de área glaciar presentó valores superiores a 96,5% (Figura 7C).

Figura 7 A. Índice Kappa. B. Exactitud del productor. C. Exactitud del usuario. Datos presentados por cada clase para las imágenes de 1987, 1991, 1997, 2000, 2003, 2009, 2014, 2016, 2017, 2019 y 2021. 

Análisis de cambios multitemporales

Los resultados de clasificación exhibieron un mayor valor de cobertura glaciar en la escena de 1987 con un área aproximada de 2899,6 ha (Tabla 3). De igual forma la clase de suelo-roca expuesta es dominante en todas las imágenes (>50% de cobertura), debido a las condiciones topográficas del área de estudio, substrato y clima. Desde 1987 hasta el 2021 la SRC exhibe variación en cuanto a las clases definidas anteriormente, lo cual puede relacionarse con el crecimiento temporal de vegetación durante las estaciones más húmedas y fusión glaciar, cambiando durante todos los intervalos de tiempo.

Tabla 3 Valores de cobertura glaciar Sierra Nevada del Cocuy (SRC) basado en las escenas analizadas en este estudio. 

La clase de área glaciar presenta valores mínimos para el año 2009 con un área aproximada de 1488,6 ha, no obstante, y respecto a la siguiente escena (2014) existe ganancia de área glaciar en este periodo de tiempo (Figura 6). Dicho valor es aproximadamente 300 ha menor al que se registró en la última escena analizada (2021) (Tabla 3); sin embargo, cabe resaltar que los valores altos de cobertura glaciar (~1800 ha) que presenta la imagen del 2021 están sujetos a la nevada registrada en la SRC en enero del mismo año. Esta nevada según reportes locales se derritió parcialmente y se acumuló únicamente en las zonas topográficamente más elevadas, por tal motivo puede que no genere datos representativos relacionados con el aumento de la masa glaciar para las imágenes analizadas posteriormente.

La Figura 8 muestra el área glaciar en la SRC para las 11 imágenes analizadas. Los valores exhiben normalmente un patrón decreciente, con cambios en las tasas de fusión glaciar en los diferentes intervalos de tiempo. Según el ^{IDEAM (2001)}, este retroceso toma lugar principalmente en los bordes del glaciar con una tasa aproximada de pérdida de 15 m/año. Pese a esta tendencia, no pueden descartarse los valores de ganancia de área glaciar durante intervalos específicos (e.g., 2009-2014, 2019-2021). Esto indica la posible incidencia de los valores más bajos de exactitud durante la clasificación del 2009 y que la variación en el área glaciar de la SRC corresponde también con otros factores climáticos, los cuales no son constantes a lo largo del tiempo como por ejemplo grandes anomalías en el fenómeno de El Niño y La Niña.

Figura 8 Cambios en el área glaciar de la Sierra Nevada del Cocuy (SRC) en el intervalo 1987-2021. 

En cuanto al cambio en las coberturas, el aumento del área glaciar en el período 2009-2014 se debe a un incremento en el régimen de lluvias asociado al fenómeno de La Niña. Se puede apreciar que el área glaciar se convirtió principalmente en la clase sueloroca, lo cual se confirma por los amplios depósitos de morrenas que se encuentran circundando actualmente la SRC (Figuras 2 y 5). De igual forma, se observa que la mayor pérdida de área glaciar se presenta sobre la vertiente occidental, además, las áreas glaciares de pequeña extensión prácticamente desaparecieron.

Condiciones meteorológicas

La zona de estudio presenta un ciclo anual bimodal en la precipitación (^{Poveda et al., 2006}; ^{Rekowsky et al., 2018}), que depende en gran medida del comportamiento de la ZCIT, en la que los vientos alisios del NE y SE, con sus masas de aire cálido y cargados de humedad chocan entre sí en la cordillera. Allí se condensan y precipitan hacia la vertiente oriental de la SRC, creando un efecto de sombra hacia la parte occidental de la cordillera, determinando así el régimen de lluvias (^{Muñoz et al., 2005}). Existe una clara diferenciación entre el régimen pluviométrico del costado oriental y occidental de la Cordillera Oriental. La región oriental de la SRC tiene un régimen de lluvias monomodal (Figura 9A) y bastante estacional principalmente hacia las sabanas, con un periodo de lluvias entre abril y noviembre y uno seco entre diciembre y marzo. Este comportamiento se debe a que esta zona recibe mayor influencia de los alisios del NE, por lo que en la época de invierno del hemisferio norte el aire viene frío y, por tanto, con baja humedad.

Figura 9 Valores de precipitación estaciones características de cada uno de los regímenes climatológicos, entre los años 2005-2009 de la SRC. A. Valores de precipitación mensual anual, estación Saravena - vertiente oriental. B. Valores de precipitación mensual anual, estación Chiscas - vertiente occidental. Modificado de ^{Muñoz et al. (2005)}. 

El régimen de lluvias hacia el occidente de la SRC es bimodal (Figura 9B), por el movimiento de la ZCIT. Los periodos húmedos ocurren entre abril-mayo y septiembre - octubre, cuando la ZCIT está pasando sobre la región, y los periodos secos entre diciembre - marzo y junio - julio, cuando está en su posición más meridional y septentrional respectivamente (^{Muñoz et al., 2005}). El flanco occidental del área cubre los municipios de El Cocuy, Chita, Chiscas y parte de Güicán. Este flanco se localiza en una sombra de lluvia o sotavento, por lo que en general la precipitación es menor en la vertiente occidental al no recibir la influencia directa de los alisos (^{Muñoz et al., 2005}), lo que explica por qué el retroceso glaciar se presenta principalmente en esta vertiente.

De igual forma, las altas montañas en los Andes colombianos se caracterizan por poseer un clima frío con fluctuaciones en los valores medios de temperatura anual, brillo solar y humedad relativa, esta última al parecer con valores más altos en cercanías a la SRC (Figura 10) (^{Silva, 2010}; ^{IDEAM, 2012}; ^{Rekowsky et al., 2018}). A pesar de que para Colombia los datos meteorológicos exhiben un aumento en los valores medios anuales de precipitación, las estaciones cercanas a la SRC revelan una tendencia de disminución en los valores de precipitación y un aumento significativo en la temperatura en los últimos 30 años (^{Rekowsky et al., 2018}). Así, se entiende que los procesos de fusión glaciar se relacionan con cambios climáticos en diferentes escalas temporales asociados con fenómenos de alta variabilidad climática como El Niño y La Niña.

Figura 10 Principales factores meteorológicos en el área de estudio. A. Precipitación media anual. B. Temperatura media anual. C. Brillo solar media diaria anual. D. Humedad relativa media anual. La caja roja representa el área de estudio en la Sierra Nevada del Cocuy (SRC) (modificado de ^{IDEAM (2010a)}). 

Problemas de exactitud para las clases en el análisis multitemporal

Las exactitudes obtenidas para las clasificaciones de las imágenes fueron satisfactorias, pero algunas consideraciones sobre el número y la distribución espacial de los puntos de control deben ser tomadas en cuenta. La recopilación de los GCP está restringida por la accesibilidad a la SRC y son relevantes únicamente para la escena del 2021. Dichos puntos se encuentran limitados a los bordes del glaciar y solo fueron tomados sobre la vertiente occidental, por tal motivo, la distribución espacial no es ideal estadísticamente. Debido a esto, es probable que algunos de los puntos presenten autocorrelación espacial haciendo que la concordancia media sea mayor (^{Guiot y de Vernal, 2011}; ^{Souza-Filho et al., 2018}).

En cuanto a las otras escenas, la validación se efectuó apenas con imágenes de alta resolución disponibles para los años 2009, 2014, 2016, 2017 y 2019, lo que genera un mayor grado de incertidumbre en las imágenes adquiridas previamente a estos datos; sin embargo, los valores generales de exactitud sugieren que la estrategia de validación fue consistente para todas las imágenes. La exactitud global y el índice Kappa fueron calculados con la finalidad de medir el grado de concordancia con los puntos de control. No se llevó en consideración información sobre la naturaleza de los desacuerdos como son el de asignación y el cuantitativo, los cuales son importantes para considerar diferencias de áreas de clase entre dos escenas (^{Souza-Filho et al., 2018}).

La exactitud del productor, así como el error de omisión, muestran qué tan bien fueron clasificados los puntos de control en las imágenes analizadas (^{Foody, 2002}; ^{Souza-Filho et al., 2018}). De acuerdo con estos valores, las exactitudes del productor fueron similares en valor y con una tendencia uniforme especialmente para la cobertura de área glaciar, sin embargo, existen valores significativamente altos de error de omisión asociados a las clases de agua y vegetación. Esto puede deberse a la cantidad de puntos tomados para la validación de estas clases o a respuestas espectrales similares, como por ejemplo la encontrada entre la vegetación de páramo de pequeño porte y el suelo-roca. Esta consideración no representa problemática para este estudio debido a que el enfoque principal es cuantificar el área glaciar y en comparación a las demás clases, esta presentó los valores de exactitud más altos y de error más bajos (<3%). Los valores de exactitud del usuario y el error de omisión indican qué tan genuinas son las clases definidas y la proporción de las diferentes clases que han sido asignadas respectivamente (Foody, 2002). Las medidas, a pesar de ser uniformes presentaron variación entre las clases y los años analizados (1987-2021), sugiriendo que en ocasiones ciertas clases incluían proporciones importantes de otras clases, lo cual ocurrió principalmente en la vegetación y el agua. Además, cabe resaltar que las imágenes usadas no presentan nubes sobre el área glaciar lo que corrobora la exactitud de los resultados obtenidos.

Variación multitemporal del área glaciar

En estudios previos, el mapeo del área glaciar en la SRC fue realizado empleando el NDSI el cual presenta una distinción representativa para la delimitación de áreas glaciares en comparación con las imágenes de composición colorida (^{Bautista y Correa, 2018}; ^{Rekowsky et al., 2018}); sin embargo, en la mayoría de los casos este procedimiento responde a la creación de máscaras de nieve, y cálculo de las áreas correspondientes, lo cual limita los resultados a procesos de interpretación visual. La metodología usada en este estudio ofrece la posibilidad de incluir durante la clasificación una mayor cantidad de información espectral de las bandas de los sensores seleccionados, y no restringir la información del modelo de clasificación, únicamente a las bandas utilizadas en el índice NDSI. Las bandas espectrales permiten hacer una caracterización apropiada de las coberturas (^{Varade et al., 2017}), facilitando la distinción entre las clases propuestas. La diferenciación entre la clase suelo-roca y la clase de área glaciar, se basa en un pico de reflectancia para la nieve en la longitud de onda del visible (^{Painter et al., 1998}), y una alta absorción en el infrarrojo cercano (^{Domine et al., 2006}); mientras, en la clase suelo-roca existe un pico de reflectancia en el infrarrojo cercano y en el infrarrojo de onda corta (^{Dong, 2018}).

El problema en la delimitación del área glaciar radica principalmente en la resolución espacial de las imágenes usadas, pues debido al tamaño del píxel, cuerpos glaciares de pequeña extensión no se toman en consideración. Aun así, cuerpos menores a esta área de cobertura (900 m2) localizados en los bordes del glaciar, son más susceptibles a los efectos causados por los aumentos de temperatura y por ende tienden a presentar fusión, por lo que posiblemente no sean constantes a lo largo del tiempo e imperceptibles para la escala de estudio.

De acuerdo con ^{Rekowsky et al. (2018)} la retracción de las áreas glaciares en la SRC se presentó hacia los bordes del cuerpo glaciar, lo cual es coherente con los resultados obtenidos en este estudio (Figura 11). Dicha variación está relacionada con que estos bordes se encuentran normalmente en alturas más bajas, donde los valores de temperatura son mayores y el espesor de hielo no es significativo. Así, varios de los cuerpos glaciares que conforman la SRC enfrentaron un proceso de fragmentación, como ya había sido documentado anteriormente (Figura 11) (^{Herrera y Ruiz, 2008}). En esta dinámica de retroceso es evidente que además de la disminución del área glaciar principal, se han extinto otros picos menores tales como el Alto Peña Colorada y la Sierra Nievesita, indicando que la extensión de los glaciares tropicales es un factor determinante en las tasas de retroceso (^{Florez, 1991}; ^{Herrera y Ruiz, 2008}). De acuerdo con estudios sedimentológicos de las morrenas asociadas a la SRC, se indica que la extensión glaciar para 1850 era de aproximadamente 148700 ha (Figura 11) con una cobertura dominante sobre la vertiente occidental.

Figura 11 Variación del área glaciar de la Sierra Nevada del Cocuy (SRC) para los años 1850, 1987, 1997, 2009 y 2019. 

Al occidente de la SRC se presentan pendientes más suaves, revelando que la topografía del terreno controla también el espesor de hielo y la estabilidad del glaciar (Figura 12), el cálculo de las pendientes permite definir que en áreas menos inclinadas de la SRC el hielo y la nieve se extienden favoreciendo su desplazamiento, mientras que en la vertiente oriental, se presentan pendientes más abruptas que acaban en depresiones, lo que permite la acumulación de la masa glaciar (^{IDEAM, 2001}; ^{Herrera y Ruiz, 2008}). Según la tendencia de retroceso glaciar calculada por ^{Herrera y Ruiz (2008)}, el glaciar de la SRC desaparecería aproximadamente en el año 2075, presentando mayores cambios en la vertiente occidental. Esta hipótesis concuerda con el trabajo de ^{Quintero (2018)} que define la existencia de corrientes húmedas sobre la vertiente oriental, la cual, a pesar de ser menor en extensión, parece ser más estable a lo largo del tiempo.

Figura 12 Mapa de pendientes de la Sierra Nevada del Cocuy (SRC). 

La SRC enfrenta un proceso de retroceso glaciar de acuerdo con estudios previos (Tabla 4), sin embargo, la tasa de retroceso es diferente para cada uno de estos, lo cual depende en gran medida de las imágenes seleccionadas, el tiempo de adquisición de las escenas, el método de interpretación, el intervalo de tiempo y las condiciones meteorológicas. Una ecuación lineal con pendiente negativa de -46,195, obtenida a través de la regresión generada con el conjunto de datos (Tabla 4, Figura 13) sugiere una tendencia de retroceso, presentando un coeficiente de determinación de ~0,85 lo cual indica que el modelo se ajusta muy bien a los datos analizados. Cabe resaltar que la mayor dispersión se presenta en las áreas calculadas antes de 1990 lo cual puede deberse a la disponibilidad de productos de sensores remotos (incluyendo la resolución espacial, espectral y temporal de estos) y la exactitud de los datos con que las áreas glaciares fueron calculadas cuando aún no se disponía de imágenes satelitales.

Tabla 4 Valores de retroceso glaciar de la Sierra Nevada del Cocuy (SRC). 

Figura 13 Retroceso glaciar de la Sierra Nevada del Cocuy (SRC) desde 1955-2021. 

De acuerdo con la tendencia de retroceso calculada a partir de todos los datos analizados en este estudio y los recopilados de otros autores (i.e., ^{Florez, 1991}; ^{Guillen et al., 2004}; ^{Morris et al., 2006}; ^{Ceballos et al., 2006}, ²⁰¹⁰; ^{Herrera y Ruiz, 2008}; ^{Peña, 2015}; ^{Bautista y Correa, 2018}; ^{Rekowsky et al., 2018}) se estima que la SRC perderá toda su cobertura glaciar aproximadamente en el año 2048. Esto se ajusta con los resultados obtenidos por ^{Guillen et al. (2004)} y ^{Morris et al. (2006)} en las tasas de retroceso glaciar calculadas ^{por Herrera y Ruiz (2008)}, en las que se prevé la desaparición del glaciar para los años 2053 y 2040 respectivamente. Si se tienen en cuenta únicamente los resultados obtenidos en este estudio, se prevé que la extinción del glaciar ocurrirá en el año 2058, lo cual concuerda perfectamente con los resultados obtenidos por ^{Ceballos et al. (2006)}. A pesar de la concordancia en los valores encontrados de retroceso, no puede descartarse la incidencia de otros factores en las tasas de pérdida glaciar para la SRC. Por ejemplo, los datos obtenidos por ^{Rodríguez (2015)} indican variaciones en los valores de precipitación, los cuales junto con la temperatura, topografía y extensión del glaciar podrían acelerar o retardar las tasas de retroceso. Esto sugiere que el retroceso de la SRC no depende únicamente del tiempo transcurrido entre dos escenas analizadas, sino que incluye ciertos parámetros climáticos.

Incidencia de los factores climáticos

Los datos de referencia utilizados en este estudio para el año de 1850 (Figura 11) podrían asociarse a los efectos de la Pequeña Edad de Hielo. La Pequeña Edad de Hielo fue un periodo de enfriamiento en la Tierra, ocurrido entre 1600-1850 (^{Matthews y Briffa, 2005}) caracterizado por cambios climáticos a escala regional, con efectos moderados en el hemisferio norte (^{Boretti, 2020}). Los glaciares ecuatoriales registran el máximo glaciar en el siglo XVII, y un retroceso constante del área glaciar durante los siglos XVIII a XIX (^{Rabatel et al., 2005}), esto es contrario a lo ocurrido en el hemisferio norte donde los glaciares continuaban expandiéndose en estos dos siglos (^{Rowan, 2017}). Por lo tanto, no se puede asegurar una relación directa entre la Pequeña Edad de Hielo y los datos registrados por el ^{IDEAM (2019)} para la SRC.

Los valores medios anuales de temperatura y precipitación mostrados por el ^{IDEAM (2010b)} para la estación Sierra Nevada de El Cocuy sugieren variaciones importantes en estos parámetros climáticos desde 1974 hasta 2007. Así, se evidencia una disminución progresiva en las tasas medias de precipitación (Figura 14A) y un aumento en las temperaturas medias anuales (Figura 14B). Estas variaciones se relacionan directamente con el cambio climático, que normalmente exhibe esta tendencia de aumento de las temperaturas debido al calentamiento global y el efecto invernadero asociado (^{Valdez- Cepeda et al., 2003}). Esta tendencia parece responder de manera más o menos lineal, aunque puede incluir a su vez ciertos parámetros de variabilidad climática que no se relacionan directamente con el cambio climático. Dentro de estos parámetros, la actividad de El Niño- Southern Oscillation (ENSO) controla en gran medida los valores de precipitación y temperatura en la SRC, como resultado del enfriamiento y calentamiento de las aguas superficiales del pacífico, el debilitamiento de los vientos alisios en el este y el desplazamiento del núcleo de convección profunda (^{IDEAM, 2007}).

Figura 14 Factores meteorológicos. A. Valores de precipitación media anual, estación Sierra Nevada El Cocuy (^{IDEAM, 2010b}). B. Valores de temperatura media anual, estación Sierra Nevada El Cocuy (^{IDEAM, 2010b}). C. Dinámica del ENSO. Los eventos del ENSO se clasifican en débiles (con una anomalía de Sea Surface Temperature (SST) de 0,5 a 0,9), moderados (1,0 a 1,4), fuertes (1,5 a 1,9) y muy fuertes (≥ 2,0), por encima de la anomalía de 0,5° para eventos cálidos (El Niño) y por debajo de -0,5° para eventos fríos (La Niña). Modificado de ^{Null (2021)}. 

El fenómeno de El Niño se considera como una época en la que predominan condiciones más secas y cálidas. Los años en los que se reportó una fuerte incidencia de El Niño, presentan grandes valores de retroceso glaciar en la SRC entre las escenas analizadas (e.g., 1987-1991, 1991-1997, 2000-2003, 2003-2009, 2014- 2019) (Figura 14C). Por el contrario, periodos donde el fenómeno de La Niña fue más fuerte indican una recuperación del área glaciar (e.g., 2009-2014, 2019- 2021) (Figura 14C), lo que se debe al predominio de condiciones más húmedas y mayores tasas de precipitación.

La variabilidad de la precipitación en los Andes tropicales está relacionada con la ocurrencia del ENSO, por lo cual el fenómeno de El Niño ocurrido en 1997/1998, indujo un balance de masa negativo debido a temperaturas más altas y a la reducción de las precipitaciones, seguido de un lento retroceso durante las fuertes condiciones de La Niña en 1999 (^{Ceballos et al., 2006}). En el caso de los glaciares colombianos, un aumento en la temperatura del aire y la radiación solar fueron los principales mecanismos detrás del rápido retroceso de los glaciares, más que los cambios en las precipitaciones (^{Euscategui, 2002}). Este supuesto debe ser validado considerando el hecho de que un aumento de temperatura puede desencadenar un cambio en la fase de precipitación de sólido (nevadas) a líquido (lluvia), lo que a su vez favorece el retroceso de los glaciares por reducción del albedo (^{Veettil, 2017}).