Modelamiento Geoestadístico de datos de temperatura y resistividad eléctrica para la caracterización del sistema geotérmico del Volcán Azufral

Área de estudio: Área geotérmica del Volcán Azufral
Variables analizadas: Sondeos superficiales de temperatura (SST) y Sondeos Eléctricos Verticales (SEV)
Autor: Gilbert Fabian Rodriguez Rodriguez
Fecha: 21/11/2025

---

0. Planteamiento del problema

Colombia, por su ubicación en el cinturón de fuego del Pacífico, concentra una importante anomalía térmica que favorece la presencia de recursos geotérmicos. Reconocimientos tempranos realizados por OLADE (1982) y CHEC (1983) identifican al menos diez regiones de interés, entre las cuales el volcán Azufral se destaca como una de las áreas prioritarias. Este edificio volcánico activo, cuya historia eruptiva reciente está respaldada por dataciones de depósitos con edades entre 17.500 y 280 años, alberga un sistema hidrotermal de alta entalpía cuya dinámica está controlada por la interacción entre el vulcanismo, las estructuras tectónicas y la circulación de fluidos, lo que lo convierte en un escenario estratégico para la exploración geotérmica.

Desde finales de la década de 1990, el Servicio Geológico Colombiano (SGC) ha adelantado estudios geológicos, geoquímicos y geofísicos en el Azufral, a partir de los cuales se ha planteado modelos conceptuales preliminares del sistema. Este modelo propone la existencia de un flujo ascendente (upflow) de fluidos hidrotermales en dirección noreste, que se dispersa lateralmente (outflow) favorecido por estructuras con la misma orientación y se acumula en un reservorio profundo a 6–8 km del cráter.

El desarrollo geotérmico en Colombia se enmarca en políticas energéticas que promueven el uso de fuentes no convencionales, respaldadas por leyes, decretos y planes de desarrollo. En este contexto, sistemas como el Nevado del Ruiz, Paipa y Tufiño–Chiles–Cerro Negro han registrado avances significativos, mientras que en el Azufral persiste la necesidad de integrar nuevas técnicas de exploración que fortalezcan la confiabilidad del modelo conceptual existente. La ausencia de un modelo geoestadístico geofisico robusto integrado con datos de temperatura genera incertidumbre en la identificación de alguna de las componentes del sistema geotérmico como la capa sello, las zonas de acumulación de fluidos para la evaluación de su potencial, por lo que resulta indispensable aplicar métodos y análisis de la estadística espacial para optimizar la caracterización del sistema y favorecer un aprovechamiento eficiente de este recurso energético.

---

1. Descripción del conjunto de datos

El conjunto de datos espaciales utilizados corresponde a dos variables observadas en el área geotérmica del volcán Azufral, ubicado en la zona sur-occidental de Colombia, en cercanías al municipio de Túquerres, departamento de Nariño. Estas observaciones se obtienen mediante dos métodos de exploración del subsuelo.

La primera variable corresponde a los sondeos superficiales de temperatura , que consisten en mediciones realizadas a dos profundidades de investigación: 20 cm y 150 cm bajo la superficie terrestre. Estas mediciones se llevan a cabo mediante sondas equipadas con sensores en un extremo, diseñados para registrar las variaciones de temperatura en el subsuelo.

El segundo método corresponde a los sondeos eléctricos verticales (SEV), empleados para determinar la resistividad eléctrica de las rocas. Este procedimiento se realiza mediante un arreglo de electrodos instalados en la superficie, a través de los cuales se hace circular una corriente eléctrica con el fin de medir la diferencia de potencial generada entre ellos. A partir de estas mediciones, es posible estimar la resistividad eléctrica del subsuelo en cada punto de observación. Para cada SEV se consideran profundidades de investigación diferenciadas: una fija a 20 cm y otra que llega a los 150 cm de profundidad.

La finalidad de estas observaciones es caracterizar algunos de los componentes del sistema geotérmico del volcán Azufral, tales como la presencia de una capa sello, la acumulación de fluidos o la identificación de un posible reservorio geotérmico.

---

2. Variables a utilizar

Para la selección de variables se utilizó inicialmente el software ArcGIS Pro, donde se integraron dos conjuntos de datos principales: los sondeos superficiales de temperatura y los sondeos eléctricos verticales (SEV). Al realizar la visualización en el mapa, se evidenció que los SEV presentaban una menor densidad espacial en comparación con los registros de temperatura superficial. Con el fin de establecer una relación entre ambos, se aplicaron herramientas de análisis espacial que permitieron desarrollar dos procedimientos complementarios.

En primer lugar, se generaron áreas de influencia mediante buffers de aproximadamente 800 metros alrededor de cada estación de SEV, representadas en el mapa por circunferencias de color azul. Posteriormente, se implementó un análisis de proximidad para seleccionar y asociar, dentro de cada buffer, los valores de temperatura superficial más cercanos a las estaciones de resistividad. Este procedimiento permitió consolidar la información en una única tabla de atributos que integra las variables de ambos métodos.

Variable	Unidades	Nº de muestras	Tipo de variable
Temperatura a 20 cm	°C	100	Dependiente (numérica)
Temperatura a 150 cm	°C	100	Independiente (numérica)
Resistividad a 20 cm	Ω·m	100	Independiente (numérica)
Resistividad a 150 cm	Ω·m	100	Independiente (numérica)
Altura de cada punto	m	100	Independiente (numérica)

Como resultado, se obtuvo un nuevo conjunto de datos en el que cada estación de SEV combina la información de resistividad eléctrica y temperatura, ambas medidas a dos profundidades (20 cm y 150 cm). Adicionalmente, se incorporó la variable de altitud, registrada en campo mediante GPS, con lo cual se conforma una base de datos integrada que facilita la correlación espacial entre las propiedades eléctricas del subsuelo y el régimen térmico superficial.

Tabla de sondeos eléctricos verticales y temperatura superficial

ID_SEV	ID_SST	Fecha_SEV	Fecha_SST	X_Magna	Y_Magna	Res_20cm	T_20cm	Prof_1	Res_150cm	T_150cm	Prof_2	Altura
1	120	25/03/2009	20/08/2015	935860.000	612035.000	148.278	12.100	0.200	511.883	12.500	1.500	3212.276
2	19	25/03/2009	05/02/2015	936818.000	611469.000	16.283	15.400	0.200	432.481	15.700	1.500	3179.040
3	128	25/03/2009	07/02/2015	938660.000	610491.000	85.014	15.700	0.200	128.919	16.400	1.500	3006.000
4	121	24/03/2009	20/08/2015	934697.000	611855.000	251.914	12.400	0.200	468.876	12.100	1.500	3360.951
5	120	24/03/2009	20/08/2015	935566.000	611458.000	440.922	12.100	0.200	213.525	12.500	1.500	3212.276
6	46	24/03/2009	06/08/2015	936073.000	610990.000	153.882	14.100	0.200	111.777	14.800	1.500	3086.023
7	39	25/03/2009	05/08/2015	938036.000	610092.000	184.797	15.300	0.200	254.437	15.500	1.500	3017.716
8	136	26/03/2009	11/09/2015	934944.000	611190.000	548.893	13.000	0.200	1781.703	12.900	1.500	3230.964
9	41	26/03/2009	05/08/2015	937281.000	609447.000	383.526	14.400	0.200	162.909	14.100	1.500	2993.027
10	40	26/03/2009	05/08/2015	937806.000	608822.000	31.593	18.000	0.200	51.249	17.000	1.500	2939.580
11	130	27/03/2009	11/09/2015	937902.000	607514.000	299.654	19.800	0.200	180.444	18.900	1.500	2908.624
12	41	28/03/2009	05/08/2015	936380.000	609492.000	104.885	14.400	0.200	138.956	14.100	1.500	2993.027
13	131	27/03/2009	11/09/2015	936755.000	607916.000	55.849	17.200	0.200	164.279	16.700	1.500	2976.957
14	130	02/04/2009	11/09/2015	937735.000	606715.000	82.029	19.800	0.200	314.152	18.900	1.500	2908.624
15	74	27/03/2009	11/08/2015	935787.000	608946.000	134.885	14.000	0.200	65.164	14.200	1.500	3036.717
16	132	27/03/2009	11/09/2015	935844.000	608139.000	461.219	16.600	0.200	144.041	15.800	1.500	2977.930
17	139	28/03/2009	12/09/2015	934798.000	607561.000	128.295	13.800	0.200	57.279	13.800	1.500	3047.113
18	213	02/04/2009	23/09/2015	934855.000	602548.000	50.833	15.400	0.200	72.367	14.300	1.500	3132.181
19	299	31/03/2009	22/09/2016	929936.000	604575.000	932.287	13.600	0.200	511.035	14.100	1.500	3074.168
20	197	31/03/2009	21/09/2015	929060.000	605323.000	208.216	13.600	0.200	9

3. Mapa de la zona de estudio

En la imagen se observó la distribución espacial de las 100 estaciones correspondientes a los sondeos eléctricos verticales y a las mediciones superficiales de temperatura realizadas en el área geotérmica del volcán Azufral. Este conjunto de datos constituyó la base para el desarrollo del análisis geoestadístico orientado a la interpolación espacial de las variables estudiadas. Las estaciones se encontraron distribuidas tanto en la franja oriental como en la franja occidental del edificio volcánico, cuyo cráter se localizó aproximadamente en las coordenadas X: 930.000 y Y: 608.000. Esta distribución espacial permitió evaluar de manera adecuada la variabilidad lateral de los parámetros geofísicos asociados al sistema geotérmico del Azufral.

4. Análisis geoestadístico univariado y bivariado de al menos dos variables espaciales.

Con el propósito de identificar la presencia de tendencias globales en las variables aleatorias —resistividad a 20 cm, temperatura a 20 cm, resistividad a 150 cm, temperatura a 150 cm y altura— se construyeron diagramas de dispersión en función de las coordenadas espaciales X y Y. Este análisis exploratorio permitió evidenciar visualmente que la altura, así como la temperatura a 20 cm y a 150 cm de profundidad, mostraron un comportamiento tendencial dentro del dominio de estudio. En particular, las temperaturas superficiales y subsuperficiales tendieron a incrementarse hacia el sector oriental, mientras que la altura presentó una tendencia decreciente respecto a la coordenada X. De manera similar, al evaluar cada variable en función de la coordenada Y, se observó una tendencia positiva tanto para la altura como para las temperaturas a 20 cm y 150 cm de profundidad, lo que sugirió un patrón espacial consistente que debía ser considerado en las siguientes etapas del análisis geoestadístico.

5. Análisis y modelamiento de cada una de las medias

Con el fin de confirmar o refutar las tendencias globales identificadas preliminarmente mediante los diagramas de dispersión, se plantearon cinco modelos de regresión lineal multivariada, cada uno asociado a una variable aleatoria distinta como variable dependiente. En particular, el modelo 1 tomó como respuesta la resistividad a 20 cm de profundidad; el modelo 2 empleó la temperatura a 20 cm; el modelo 3 utilizó la resistividad a 150 cm; el modelo 4 incorporó la temperatura a 150 cm; y el modelo 5 consideró la altura. Tras ajustar y validar cada modelo, se evaluó el aporte individual de las variables espaciales X magna y Y magna mediante el análisis de la significancia estadística de sus coeficientes.

Los resultados evidenciaron que, para el modelo 1, únicamente el coeficiente asociado a la coordenada X magna resultó estadísticamente significativo. En los modelos 2 y 3, ambas coordenadas mostraron coeficientes significativos; no obstante, en el modelo 3 el intercepto no alcanzó significancia estadística. En los modelos 4 y 5, tanto X magna como Y magna presentaron altos niveles de significancia. Este análisis se fundamentó en la evaluación de los p-valores, los cuales, para las variables significativas, se mantuvieron por debajo del umbral del 5%, permitiendo rechazar la hipótesis nula (β = 0) y confirmar la hipótesis alternativa de que los parámetros β fueron distintos de cero y, por tanto, relevantes para la explicación de la variabilidad espacial de cada respuesta.

Adicionalmente, se definió desde el inicio que los modelos 1 y 3 requerían la aplicación de una transformación logarítmica a sus variables dependientes. Esta decisión se sustentó en la observación de que, tras generar los mapas de interpolación mediante Kriging Simple, las predicciones originales en dichos modelos arrojaban valores negativos, lo cual carece de sentido físico, dado que la resistividad del subsuelo es estrictamente positiva. La transformación logarítmica garantizó que las estimaciones obtenidas en los mapas interpolados mantuvieran coherencia geofísica, preservando valores positivos y permitiendo una interpretación adecuada dentro del contexto del sistema geotérmico del volcán Azufral.

6. Selección de los semivariogramas empíricos individuales, más adecuados para la estimación del modelo teórico

Una vez validados y aceptados los modelos de regresión lineal multivariada definidos en el paso anterior, se procedió a trabajar exclusivamente con los residuales provenientes de dichos modelos. Esto se debió a que, mediante la regresión, ya se había removido la tendencia global presente en cada una de las variables aleatorias, garantizando que los residuales representaran únicamente la estructura local de variabilidad espacial. Por lo tanto, estos residuales se utilizaron como insumo base para la construcción de los semivariogramas empíricos individuales.

Con el fin de caracterizar de manera integral la dependencia y autocorrelación espacial de cada variable, se generó un conjunto amplio de semivariogramas exploratorios. En primera instancia, se construyeron semivariogramas omnidireccionales para los residuales de las cinco variables analizadas. Posteriormente, se generaron semivariogramas direccionales en las orientaciones típicamente empleadas en geoestadística —45°, 90° y 135°— con el propósito de identificar posibles anisotropías geométricas o zonales en la estructura espacial.

Esta primera etapa se desarrolló empleando el estimador clásico de momentos (Matheron), el cual es sensible a la presencia de valores atípicos. Con el fin de evaluar la robustez de la estructura espacial identificada, se replicó el procedimiento utilizando dos estimadores resistentes: el estimador de Cressie-Hawkins, basado en un promedio robusto, y posteriormente un estimador adicional basado en la mediana. En ambos casos se reconstruyeron los semivariogramas omnidireccionales y direccionales para cada conjunto de residuales. Esta estrategia permitió comparar cómo distintos niveles de resistencia a valores extremos afectaban la forma, el alcance y el comportamiento del semivariograma, fortaleciendo la interpretación de la estructura espacial antes de proceder al ajuste de modelos teóricos.

La selección final de los semivariogramas empíricos se fundamentó en un criterio técnico claramente definido, basado en el comportamiento esperado de la estructura espacial. Para cada una de las combinaciones generadas, se evaluó si el semivariograma presentaba un incremento coherente de la semivarianza con la distancia, evitando tendencias erráticas o decrecientes que contradijeran el principio de continuidad espacial. Adicionalmente, se valoró que la semivarianza inicial permaneciera en niveles bajos y que, en lo posible, se evidenciara de manera clara la estabilización de la curva hacia un nivel de varianza estructural o sill, lo cual es indicativo de un proceso espacial intrínsecamente estacionario. Estos criterios permitieron aceptar únicamente los semivariogramas cuya morfología reflejara un comportamiento físico y estadístico plausible, y descartar aquellos que presentaban distorsiones producto del ruido, de la dirección del muestreo o de la influencia de observaciones atípicas. Gracias a esta metodología, se garantizó que los modelos seleccionados representaran adecuadamente la dependencia espacial inherente a cada variable aleatoria.

Después de analizar de manera detallada la totalidad de los semivariogramas empíricos generados para cada conjunto de residuales, se seleccionaron los modelos que mostraron el comportamiento más estable, interpretable y coherente con la estructura espacial observada en los datos. Esta selección permitió definir, para cada variable aleatoria, el semivariograma empírico más adecuado para proceder posteriormente al ajuste del modelo teórico.

En el caso de la resistividad medida a 20 cm de profundidad (residual 1), se determinó que el semivariograma empírico más apropiado fue el omnidireccional estimado mediante el método clásico de momentos, dado que presentó una estructura espacial limpia, sin fluctuaciones abruptas y con un incremento progresivo acorde con un comportamiento geoestadísticamente consistente. Para la temperatura a 20 cm de profundidad (residual 2), también se seleccionó el semivariograma omnidireccional obtenido mediante el estimador clásico, ya que mostró una relación distancia–semivarianza congruente con la tendencia espacial esperada y sin una influencia considerable de valores extremos.

En contraste, para la resistividad a 150 cm de profundidad (residual 3) se encontró que el estimador clásico resultaba demasiado sensible a valores atípicos, lo que distorsionaba la estructura espacial. Por esta razón, se seleccionó como modelo óptimo el semivariograma omnidireccional construido mediante el estimador robusto de Cressie-Hawkins, el cual logró estabilizar la nube de semivarianzas y permitió identificar un patrón espacial más claro y confiable. Una situación similar se evidenció en la temperatura a 150 cm de profundidad (residual 4), para la cual también se seleccionó el semivariograma omnidireccional estimado mediante el método de Cressie, dada su mayor robustez frente a valores extremos y su mejor definición del comportamiento espacial.

Finalmente, para la variable altura (residual 5) se observó nuevamente que el estimador de Cressie ofrecía una caracterización más estable de la estructura espacial, por lo que se seleccionó el semivariograma omnidireccional resistente a datos atípicos como el más adecuado para el ajuste del modelo teórico. En conjunto, estas decisiones permitieron garantizar que la etapa posterior de modelación teórica se realizara sobre estructuras espaciales adecuadamente representadas, maximizando la confiabilidad de las interpolaciones mediante kriging.

7. Estimación de los semivariogramas individuales teóricos a sentimiento

Con el fin de estimar los semivariogramas teóricos individuales, se inició definiendo valores preliminares para los parámetros del modelo: efecto pepita (nugget), meseta (sill) y rango (range). Esta selección inicial es necesaria para alimentar adecuadamente los modelos teóricos considerados. En este caso, el valor preliminar del efecto pepita se tomó como el mínimo valor reportado en la función de semivarianza del variograma empírico. Para la meseta, se seleccionó el valor máximo observado en dicha función, mientras que el rango preliminar se definió como el valor máximo de la distancia registrada en el semivariograma empírico, dividido entre tres. Este mismo esquema de selección de valores iniciales para nugget, sill y range se aplicó a cada uno de los demás modelos analizados.Una vez establecidos los valores preliminares, se procedió a probar cinco modelos teóricos distintos con el fin de evaluar cuál ofrecía el mejor ajuste a cada semivariograma empírico. Los modelos considerados fueron: exponencial, esférico, gaussiano, lineal y Matérn.

Para la estimación de los parámetros se emplearon tres métodos: Mínimos Cuadrados Ordinarios (OLS), utilizando los valores iniciales de nugget, sill y range para obtener una primera aproximación del ajuste entre el modelo teórico y el semivariograma experimental.Mínimos Cuadrados Ponderados (WLS), con el propósito de mejorar la precisión del ajuste al considerar la variabilidad asociada a cada punto del semivariograma empírico.Máxima Verosimilitud (ML), implementado mediante funciones disponibles en paquetes de R, con el objetivo de obtener estimaciones óptimas de los parámetros del efecto pepita, meseta y rango para los modelos teóricos que mostraron mejor comportamiento frente al variograma experimental.Este procedimiento completo se aplicó de manera sistemática para cada una de las variables aleatorias analizadas en esta etapa del estudio.

8. Estimación de los semivariogramas individuales y del modelo teórico usando MCO, MCP, MV, CL.

Una vez obtenidos los resultados derivados de las distintas combinaciones entre modelos teóricos y métodos de estimación, se procedió a seleccionar aquellos modelos que presentaron el mejor ajuste frente a sus respectivos semivariogramas empíricos. Como resultado, se identificaron los siguientes modelos teóricos óptimos:Resistividad a 20 cm de profundidad: modelo exponencial.Temperatura a 20 cm de profundidad: modelo gaussiano.Resistividad a 150 cm de profundidad: modelo exponencial.Temperatura a 150 cm de profundidad: modelo exponencial. Altura: modelo exponencial.Es importante resaltar que todos estos modelos fueron estimados utilizando el método de Mínimos Cuadrados Ordinarios (OLS), el cual proporcionó el ajuste más adecuado para las variables analizadas.

9. Mapas de predicción de kriging acompañados de cada uno de los mapas respectivos de varianza del error de predicción

Una vez seleccionados los modelos teóricos de semivariograma que mejor se ajustaron a cada uno de los semivariogramas experimentales, se extrajeron los valores estimados de sus parámetros: la silla, el rango y el efecto pepita. Con esta información, se procedió a aplicar el método de Kriging simple para generar los mapas finales de interpolación y los correspondientes mapas de varianza de los errores de predicción.

Para cada una de las cinco variables analizadas se construyó un mapa de predicción mediante Kriging simple, acompañado de su respectivo mapa de varianza, garantizando así una caracterización completa tanto de la estimación espacial como de su incertidumbre asociada.

## [using simple kriging]

## [using simple kriging]

## [using simple kriging]

## [using simple kriging]

## [using simple kriging]

Para la construcción de los mapas de interpolación se definió una malla regular con una resolución espacial de 5 metros. Con esta resolución, la grilla generada estuvo compuesta por 2.937 filas y 2.510 columnas, lo que produjo un total de 7.371.870 puntos de predicción. Esta malla sirvió como soporte para aplicar el método de Kriging simple y obtener las superficies interpoladas y los mapas de varianza de los errores de predicción para cada una de las cinco variables analizadas.

Los resultados obtenidos mediante la aplicación del Kriging simple permiten identificar patrones espaciales relevantes en las variables geofísicas medidas en el área geotérmica del volcán Azufral.El mapa de resistividad evidencia la presencia de una capa superficial con fluidos o minerales que favorecen la reducción de resistividad a nivel somero. Esta respuesta es coherente con la hipótesis de que aguas subterráneas circulantes o niveles freáticos poco profundos pueden afectar la conductividad del suelo en los primeros 20 cm. Respecto al error de predicción, los valores más bajos se concentran en la zona cercana a las estaciones de medición, mientras que los errores incrementan progresivamente hacia áreas sin observaciones, lo que refleja la densidad y distribución espacial de los datos utilizados en la interpolación.

En la parte noroccidental del área geotérmica, los valores de temperatura son inferiores a 14 °C, fenómeno que se asocia con la presencia de fluidos superficiales y con la respuesta de resistividad observada. La correlación espacial entre ambas variables sugiere que las zonas de mayor conductividad coinciden con temperaturas más bajas, lo que refleja la influencia de la autocorrelación espacial en los niveles someros del terreno.El mapa de error muestra mínimos de predicción en las proximidades de las estaciones y un aumento de error en regiones sin información directa.

A esta profundidad, la resistividad promedio en la mayor parte del área geotérmica, especialmente en la zona baja del volcán, es inferior a 200 Ω·m. Este patrón es consistente con la conductividad observada a 20 cm y sugiere que la infiltración de aguas superficiales continúa influyendo en capas ligeramente más profundas, reafirmando la relación entre la litología local y la respuesta eléctrica. Los errores de predicción siguen la misma tendencia que a 20 cm: mínimos alrededor de las estaciones y mayores en zonas sin datos directos.

El comportamiento térmico a 150 cm mantiene las tendencias observadas a 20 cm, destacando temperaturas más bajas en el mismo sector noroccidental. Se registran valores próximos a 12 °C, particularmente sobre el cráter del volcán, zona caracterizada por páramo y altitudes superiores a 4.000 m. El mapa de error indica que los menores errores se concentran alrededor de las estaciones de medición, evidenciando la robustez del modelo de Kriging simple en las zonas con información directa.

El mapa interpolado de la variable altura refleja con precisión la morfología del volcán: Las áreas de mayor elevación, representadas con tonalidades verdes, coinciden con el cráter y las cumbres del volcán. Las zonas más bajas, representadas con colores claros, corresponden a valles y planicies, como el sector del municipio de Túquerres, con altitudes inferiores a 3.000 m s. n. m. El mapa de error asociado confirma la confiabilidad de la interpolación, con predicciones más precisas en áreas cercanas a las estaciones y resultados razonables en los sectores intermedios.

10. Mostrar ejemplos de estimación de modelos de semivariograma distintos a los que aparecen en los paquetes geoR y gstat, usando paquetes de regresión lineal y no lineal, usando matrices de ponderación fijas y también matrices de ponderación que varían en cada iteración.

Para resolver la solicitud de ajustar modelos de semivariograma distintos a los disponibles en los paquetes geoR y gstat, se implementaron explícitamente dos estructuras teóricas no estándar: el semivariograma power acotado y el rational quadratic. Estos modelos no existen como opciones predefinidas dentro de las funciones vgm() o variofit(), por lo que fue necesario definir sus ecuaciones manualmente y construir funciones que devolvieran su valor teórico para cualquier distancia ℎ.

A partir de estas expresiones analíticas, se utilizó el paquete minpack.lm, mediante la función nlsLM(), para realizar un ajuste inicial por mínimos cuadrados no lineales, empleando como insumos los lags del variograma empírico, los valores de semivarianza observados y el número de parejas por bin.

Para mejorar la estabilidad y otorgar mayor influencia a los bins con mayor número de pares, se construyó una matriz de ponderación fija, definida como es el número de parejas por bin , y en donde 𝜀 es un valor pequeño que evita inestabilidades numéricas. Esta matriz de pesos se incorporó tanto en nlsLM() como en la función objetivo empleada en optim(). Con ello se implementó explícitamente un procedimiento de Mínimos Cuadrados Ponderados (WLS) con pesos fijos. En esta etapa también se utilizó el método de optimización L-BFGS-B, permitiendo imponer restricciones como la no negatividad del nugget y de la sill.

\[ w_i = \frac{N_i}{\gamma_i^2 + \varepsilon} \]

Posteriormente, se desarrolló un esquema iterativo de reponderación, donde los pesos ya no permanecen fijos, sino que se recalculan en cada iteración utilizando el semivariograma teórico evaluado con los parámetros actualizados. Este enfoque constituye un ejemplo de Iteratively Reweighted Least Squares (IRWLS) aplicado al ajuste de variogramas no estándar.

En cada ciclo se optimiza nuevamente la función objetivo, se actualizan los parámetros (nugget,sill,𝑎),se calcula el error ponderado, y se recomputa la matriz de pesos basándose en la semivarianza teórica.

Este procedimiento genera un proceso completamente dinámico hasta alcanzar convergencia. Finalmente, se almacenan y reportan las trayectorias de los parámetros y del error ponderado, evidenciando cómo la ponderación evoluciona en cada iteración.

##          nug         sill            a 
## 5.429715e-01 5.567032e+01 1.199450e+09

## [1] TRUE

## [1] 1

##     it       nug     sill          a     wRSS
## nug  1 0.5429715 55.67032 1199450038 28.31577

11. Determinar de manera simbólica general las expresiones a optimizar en MCO, MCP, MV, CL y programarlas usando la función optim. Comparar los resultados con los obtenidos en R, geoR o gstat y geoModels usando los datos con los que cada grupo está desarrollando su informe.

En este capítulo se desarrollaron tres procedimientos complementarios orientados a formalizar el proceso de estimación de los parámetros estructurales de los modelos teóricos de semivariograma. El primero de ellos consistió en la formulación simbólica de las expresiones matemáticas asociadas a los métodos clásicos y modernos de inferencia geoestadística, específicamente los criterios de mínimos cuadrados ordinarios (OLS), mínimos cuadrados ponderados (WLS), máxima verosimilitud (ML) y máxima verosimilitud compuesta (CML). La construcción explícita de estas expresiones permite fundamentar, desde un punto de vista analítico, los procedimientos utilizados posteriormente para el ajuste de las estructuras de dependencia espacial, dado que cada uno de estos enfoques incorpora supuestos distintos sobre la distribución de los residuales, la estructura de covarianza y el comportamiento de la semivarianza empírica. Con base en lo anterior, a continuación se presentan las ecuaciones teóricas que sustentan cada uno de estos criterios de estimación.

11.1. Mínimos Cuadrados Ordinarios (MCO)

La función objetivo para MCO es:

\[ Q_{\text{MCO}}(\theta) = \sum_{k=1}^{m} \left[ \hat{\gamma}(h_k) - \gamma(h_k; \theta) \right]^2 \] ### 11.2. Mínimos Cuadrados Ponderados (MCP / WLS)

La función objetivo general para MCP es:

\[ Q_{\text{MCP}}(\theta) = \sum_{k=1}^{m} w_k \left[ \hat{\gamma}(h_k) - \gamma(h_k; \theta) \right]^2 \]

Pesos típicos:

\[ w_k = N(h_k) \]

o

\[ w_k = \frac{N(h_k)}{\hat{\gamma}(h_k)^2 + \varepsilon} \]

11.3. Máxima Verosimilitud (ML)

Suponiendo que los residuos siguen una distribución normal multivariada:

\[ z \sim N(0, \Sigma(\theta)), \]

la función de log-verosimilitud negativa (a minimizar) es:

\[ \mathcal{L}_{\text{ML}}(\theta) = \frac{1}{2} \left[ \ln|\Sigma(\theta)| + z^\top \Sigma(\theta)^{-1} z \right] \]

(La constante \(n \ln(2\pi)\) se omite porque no depende de \(\theta\)).

11.4. Máxima Verosimilitud Compuesta (CL)

Para pares \((i,j) \in P\), la log-verosimilitud compuesta es:

\[ CL(\theta) = \sum_{(i,j)\in P} \log f_{ij}(z_i, z_j; \theta) \]

donde cada término bivariado es:

\[ \log f_{ij} = -\frac{1}{2} \left[ \ln|\Sigma_{ij}(\theta)| + (z_{ij} - \mu_{ij})^\top \Sigma_{ij}(\theta)^{-1} (z_{ij} - \mu_{ij}) \right]. \]

El segundo procedimiento se centró en la implementación computacional de los distintos métodos de estimación teórica, incluyendo los esquemas de mínimos cuadrados ordinarios, mínimos cuadrados ponderados y máxima verosimilitud, con el propósito de contrastar sus resultados frente a los obtenidos mediante las funciones disponibles en los paquetes especializados de R, tales como geoR, gstat y GeoModels. La programación explícita de estos métodos permitió evaluar de manera controlada el comportamiento de los algoritmos de optimización, verificar la consistencia de los parámetros estimados y analizar las diferencias que pueden surgir entre una implementación manual y las rutinas preestablecidas de estas librerías, especialmente en la estimación del efecto pepita, el sill y el rango para diversas estructuras teóricas del semivariograma.

Estas rutinas programadas se aplicaron posteriormente a todas las combinaciones posibles de modelos teóricos previamente definidos. En particular, para el método de mínimos cuadrados ordinarios se evaluaron las estructuras de semivariograma exponencial, gaussiana, lineal, esférica y Matérn, generando un conjunto completo de ajustes para cada una de las variables analizadas. La misma batería de modelos teóricos se empleó en la implementación del esquema de mínimos cuadrados ponderados, con el fin de mantener la comparabilidad entre métodos y evaluar el efecto de la ponderación basada en el número de parejas o en la variabilidad dentro de cada bin. De manera análoga, para el procedimiento de máxima verosimilitud se utilizaron las correspondientes funciones de covarianza asociadas a los modelos exponencial, Matérn, lineal, gaussiano y esférico, permitiendo contrastar su desempeño frente a los métodos basados en mínimos cuadrados y analizar la estabilidad y eficiencia de cada estructura teórica bajo criterios de optimización estocástica.

Tablas de comparacion de los modelos 1 al 5 (Resistividad a 20cm) entre el metodo OLS usando la paqueteria de R VS OLS programado usando la funcion Optim

##        Modelo Nugget_OLS_Prog Sill_OLS_Prog Range_OLS_Prog Nugget_OLS_Opt
## 1 Exponencial       0.0000000     0.6652130       355.7894      0.4825230
## 2   Gaussiano       0.0000000     0.6627693       413.5594      0.5420099
## 3      Lineal       0.1929791     0.4627210       637.2907      0.4709126
## 4    Esférico       0.0000000     0.6628024       917.1256      0.3889861
## 5      Matérn       0.0000000     0.6652130       355.7894      0.0000000
##   Sill_OLS_Opt Range_OLS_Opt
## 1 0.2219466717      2077.427
## 2 0.1329346962      2077.427
## 3 0.0000922197      2077.426
## 4 0.2797083617      2077.427
## 5 0.6934277447      2077.426

##        Modelo Nugget_WLS_Prog Sill_WLS_Prog Range_WLS_Prog Nugget_WLS_Opt
## 1 Exponencial       0.5724466   1.505948325     54909.5333       1.061305
## 2   Gaussiano       0.6588004   0.000000000      2077.4266       1.321599
## 3      Lineal       0.6530795   0.001483825      1963.9760       1.368640
## 4    Esférico       0.0000000   0.657535581       838.4591       1.071230
## 5      Matérn       0.5724467   1.505951350     54909.6448       1.320940
##   Sill_WLS_Opt Range_WLS_Opt
## 1 0.9816822847      2077.432
## 2 0.7678457753      3279.497
## 3 0.0002141063      2077.407
## 4 0.7598654999      2077.431
## 5 0.9753778397      2077.426

##        Modelo Nugget_OLS_Prog Sill_OLS_Prog Range_OLS_Prog Nugget_OLS_Opt
## 1 Exponencial       0.0000000     0.9681727       494.8499      0.0000000
## 2   Gaussiano       0.2328200     0.7207784       574.1737      0.3068998
## 3      Lineal       0.3683180     0.5955331      1345.6380      0.5654861
## 4    Esférico       0.2685342     0.6933206      1560.1305      0.2687858
## 5      Matérn       0.0000000     0.9681727       494.8499      0.0000000
##   Sill_OLS_Opt Range_OLS_Opt
## 1 0.9681778320      494.8868
## 2 0.6502390114      651.4655
## 3 0.0001946664     2077.4248
## 4 0.6930836704     1560.9080
## 5 1.0235883150     2077.4254

##        Modelo Nugget_OLS_Prog Sill_OLS_Prog Range_OLS_Prog Nugget_OLS_Opt
## 1 Exponencial               0  1.228920e+05       5077.347      0.0000000
## 2   Gaussiano               0  1.265290e+00        915.680      0.0000000
## 3      Lineal               0  1.267984e+00       1645.370      0.1019251
## 4    Esférico               0  1.267145e+00       2122.739      0.0000000
## 5      Matérn               0  1.385372e+00       1308.235      0.0000000
##   Sill_OLS_Opt Range_OLS_Opt
## 1 1.3853234564     1308.0145
## 2 1.2610662303      884.8056
## 3 0.0005733594     2077.4242
## 4 1.2671541106     2122.8728
## 5 1.4673233367     2077.4260

##        Modelo Nugget_OLS_Prog Sill_OLS_Prog Range_OLS_Prog Nugget_OLS_Opt
## 1 Exponencial         0.10000      3000.000       8000.000         0.0000
## 2   Gaussiano       290.70451      4409.203       4516.776       350.0410
## 3      Lineal         0.00000      4762.291       7467.484         0.0000
## 4    Esférico         0.00000     47856.340     112848.764         0.0000
## 5      Matérn        47.28353    435513.472     651209.394       296.0338
##   Sill_OLS_Opt Range_OLS_Opt
## 1 1.269475e+07  19345417.579
## 2 4.795645e+03      4982.638
## 3 6.570098e-01      6256.905
## 4 1.938853e+06   4425496.206
## 5 5.443947e+03      1711.254

Tablas de comparacion de los modelos 1 al 5 (Resistividad a 20cm) entre el metodo WLS usando la paqueteria de R VS WLS programado usando la funcion Optim

##        Modelo Nugget_WLS_Prog Sill_WLS_Prog Range_WLS_Prog Nugget_WLS_Opt
## 1 Exponencial       0.5724466   1.505948325     54909.5333      0.5663954
## 2   Gaussiano       0.6588004   0.000000000      2077.4266      0.6144480
## 3      Lineal       0.6530795   0.001483825      1963.9760      0.6364683
## 4    Esférico       0.0000000   0.657535581       838.4591      0.5591221
## 5      Matérn       0.5724467   1.505951350     54909.6448      0.5891228
##   Sill_WLS_Opt Range_WLS_Opt
## 1   0.09109027      2077.598
## 2   0.02665240      2077.628
## 3   0.00000000      2076.456
## 4   0.08032067      2077.475
## 5   0.09734505      2077.520

##        Modelo Nugget_WLS_Prog Sill_WLS_Prog Range_WLS_Prog Nugget_WLS_Opt
## 1 Exponencial       0.9118615     1.3080239       2788.470      0.6837183
## 2   Gaussiano       1.2452099     0.8203835       2425.285      1.1437745
## 3      Lineal       1.1573680     0.8797922       4608.953      0.8233237
## 4    Esférico       0.9947979     1.0270214       5213.782      0.9437814
## 5      Matérn       0.9118616     1.3080239       2788.471      1.1643764
##   Sill_WLS_Opt Range_WLS_Opt
## 1 1.3665379427      2077.679
## 2 0.8917011241      3075.866
## 3 0.0004547074      2076.439
## 4 1.0475269380      5641.378
## 5 1.1660094769      2077.584

##        Modelo Nugget_WLS_Prog Sill_WLS_Prog Range_WLS_Prog Nugget_WLS_Opt
## 1 Exponencial       0.6825854     0.3161665       1549.602      0.6954182
## 2   Gaussiano       0.7737058     0.2044435       1802.381      0.7648460
## 3      Lineal       0.6936068     0.2723778       1976.944      0.6660582
## 4    Esférico       0.4875145     0.4759414       1818.016      0.4153641
## 5      Matérn       0.6825853     0.3161666       1549.601      0.8167379
##   Sill_WLS_Opt Range_WLS_Opt
## 1 0.3027496746      2077.548
## 2 0.1988654865      2077.561
## 3 0.0001350076      2077.103
## 4 0.5267611571      1713.942
## 5 0.2259675909      2077.478

##        Modelo Nugget_WLS_Prog Sill_WLS_Prog Range_WLS_Prog Nugget_WLS_Opt
## 1 Exponencial      0.80445767   134.0667887    1181803.558      0.0000000
## 2   Gaussiano      0.88996483     0.6713424       4403.972      0.0000000
## 3      Lineal      0.39075128     0.8726137       1954.551      0.1212866
## 4    Esférico      0.00385043     1.2553026       2034.701      0.0000000
## 5      Matérn      0.80445771   134.0679090    1181809.976      0.0000000
##   Sill_WLS_Opt Range_WLS_Opt
## 1 1.3875620563     1658.6450
## 2 1.1951184195     1008.9971
## 3 0.0005128826     2077.0765
## 4 1.2122581873     2649.6691
## 5 1.2439050146      554.8771

##        Modelo Nugget_WLS_Prog Sill_WLS_Prog Range_WLS_Prog Nugget_WLS_Opt
## 1 Exponencial          0.0000    502039.526     733465.285         0.0000
## 2   Gaussiano        407.2337      4371.359       4741.435       233.1037
## 3      Lineal          0.0000      3778.985       5757.327         0.0000
## 4    Esférico          0.0000     56465.080     133644.557         0.0000
## 5      Matérn          0.0000    502042.725     733468.145       179.4841
##   Sill_WLS_Opt Range_WLS_Opt
## 1 3.751140e+06   6240048.667
## 2 3.459875e+03      3683.643
## 3 6.005828e-01      6581.329
## 4 1.571849e+07  39259459.173
## 5 4.817134e+03      2416.353

Tablas de comparacion de los modelos 1 al 5 (Resistividad a 20cm) entre el metodo ML usando la paqueteria de R VS ML programado usando la funcion Optim

##        Modelo Nugget_ML_Prog Sill_ML_Prog Range_ML_Prog Nugget_ML_Opt
## 1 Exponencial      0.6494326  0.022443785      2077.427 -1.421085e-14
## 2   Gaussiano      0.6625620  0.006819926      2077.427  4.312915e+02
## 3      Lineal      0.6782158 -0.009816632      2077.427  7.503059e+01
## 4    Esférico      0.6781736 -0.009615656      2077.427  6.607950e-01
## 5      Matérn      0.6494326  0.022443785      2077.427  1.233842e+02
##    Sill_ML_Opt Range_ML_Opt
## 1 11014.129891    11400.222
## 2  8923.066601     3455.590
## 3     1.352149     2097.489
## 4     0.000000     2077.425
## 5  4762.290785     2077.427

##        Modelo Nugget_ML_Prog Sill_ML_Prog Range_ML_Prog Nugget_ML_Opt
## 1 Exponencial      0.6046454    1.4507207      2077.427 -1.421085e-14
## 2   Gaussiano      1.0807126    0.7195626      2077.427  4.312915e+02
## 3      Lineal      0.9031862    0.8377420      2077.427  7.503059e+01
## 4    Esférico      0.2581466    1.5494396      2077.427  6.607950e-01
## 5      Matérn      0.6046454    1.4507207      2077.427  1.233842e+02
##    Sill_ML_Opt Range_ML_Opt
## 1 11014.129891    11400.222
## 2  8923.066601     3455.590
## 3     1.352149     2097.489
## 4     0.000000     2077.425
## 5  4762.290785     2077.427

##        Modelo Nugget_ML_Prog Sill_ML_Prog Range_ML_Prog Nugget_ML_Opt
## 1 Exponencial      0.5945484    0.3563385      2077.427 -1.421085e-14
## 2   Gaussiano      0.6987660    0.2160619      2077.427  4.312915e+02
## 3      Lineal      0.5871735    0.3249316      2077.427  7.503059e+01
## 4    Esférico      0.4064470    0.5106158      2077.427  6.607950e-01
## 5      Matérn      0.5945484    0.3563385      2077.427  1.233842e+02
##    Sill_ML_Opt Range_ML_Opt
## 1 11014.129891    11400.222
## 2  8923.066601     3455.590
## 3     1.352149     2097.489
## 4     0.000000     2077.425
## 5  4762.290785     2077.427

##        Modelo Nugget_ML_Prog Sill_ML_Prog Range_ML_Prog Nugget_ML_Opt
## 1 Exponencial      0.3710857    1.1935714      2077.427 -1.421085e-14
## 2   Gaussiano      0.7287749    0.6605486      2077.427  4.312915e+02
## 3      Lineal      0.5317620    0.7576317      2077.427  7.503059e+01
## 4    Esférico      0.1701055    1.1190475      2077.427  6.607950e-01
## 5      Matérn      0.3710857    1.1935714      2077.427  1.233842e+02
##    Sill_ML_Opt Range_ML_Opt
## 1 11014.129891    11400.222
## 2  8923.066601     3455.590
## 3     1.352149     2097.489
## 4     0.000000     2077.425
## 5  4762.290785     2077.427

##        Modelo Nugget_ML_Prog Sill_ML_Prog Range_ML_Prog Nugget_ML_Opt
## 1 Exponencial     -604.55061     4464.318      2077.427 -1.421085e-14
## 2   Gaussiano      354.51855     4360.418      2077.427  4.312915e+02
## 3      Lineal       79.73659     2890.905      2077.427  7.503059e+01
## 4    Esférico     -890.49995     4285.028      2077.427  6.607950e-01
## 5      Matérn     -604.55061     4464.318      2077.427  1.233842e+02
##    Sill_ML_Opt Range_ML_Opt
## 1 11014.129891    11400.222
## 2  8923.066601     3455.590
## 3     1.352149     2097.489
## 4     0.000000     2077.425
## 5  4762.290785     2077.427

A partir de los resultados obtenidos en cada una de las tablas, tanto para los modelos ajustados mediante las funciones de los paquetes especializados como para aquellos estimados con la función optim(), fue posible establecer criterios comparativos que permiten discernir cuál método de ajuste es más adecuado en cada situación. Aunque ambos procedimientos implementan el mismo enfoque de mínimos cuadrados ordinarios, las diferencias en los algoritmos de optimización generan variaciones apreciables en la estimación de los parámetros, lo que obliga a realizar una revisión detallada de su comportamiento.

El primer criterio de selección se basó en identificar casos en los que, independientemente del modelo teórico de semivariograma utilizado, la función optim() produce exactamente el mismo valor para el rango. Este patrón sugiere que el proceso de optimización no logró converger hacia una solución adecuada y quedó atrapado en un punto no óptimo del espacio de parámetros. La evidencia se hace explícita cuando dicho valor se repite sistemáticamente en todos los modelos teóricos considerados, lo que indica una falla en la búsqueda del mínimo de la función objetivo y, por tanto, una estimación poco confiable del parámetro.

El segundo criterio, más flexible, consiste en evaluar la proximidad entre las estimaciones generadas por ambos métodos. Cuando los valores de rango, sill y efecto pepita obtenidos mediante la paquetería y mediante optim() son prácticamente equivalentes —como ocurre en el modelo número 3 y, de forma similar, en el modelo número 4—, cualquiera de las dos alternativas proporciona resultados comparables sin afectar la interpretación geoestadística ni la estructura de dependencia espacial. En estas situaciones, la elección entre uno u otro método carece de impacto práctico en el ajuste final del semivariograma.

Considerando estos dos criterios, se optó en varios casos por privilegiar las estimaciones provenientes de las funciones incorporadas en los paquetes geoestadísticos, dado que la función optim() mostró dificultades recurrentes para identificar el valor óptimo de ciertos parámetros, especialmente cuando repetía soluciones idénticas entre modelos teóricos distintos. Sin embargo, en situaciones donde las estimaciones fueron estables y prácticamente coincidentes, ambos métodos resultan equivalentes y pueden utilizarse indistintamente sin que se esperen diferencias sustantivas en los resultados.

10. Kriging espacio-tiempo para una de las variables. Repetir para este contexto todos los pasos en 1, análisis de la media, estimación empírica y teórica de los modelos de semivarianza utilizando todos los métodos y comparando todos los modelos y generar varios mapas de predicción horizontes cercanos.

En este nuevo capítulo, dedicado a la construcción del kriging espacio–temporal mediante análisis geoestadístico, se inició el proceso cargando la base de datos correspondiente a 350 estaciones de medición de temperatura superficial, registradas específicamente a 150 cm de profundidad. Para desarrollar el análisis, se consideró como variable aleatoria únicamente la temperatura observada; adicionalmente, para cada medición se incorporaron las coordenadas espaciales X y Y en columnas independientes, así como la fecha de adquisición de cada registro.

Dado que la fecha estaba originalmente en formato calendario (año–mes–día), se requirió transformarla a una escala temporal continua. Para ello, se definió como fecha de referencia el 31 de enero de 2015, correspondiente al primer día en que se inició la adquisición de este conjunto de datos. A partir de esa fecha base, se calculó el número de horas transcurridas hasta el momento de cada una de las 350 observaciones. Este valor en horas permitió representar el tiempo como una variable cuantitativa continua, requisito fundamental para la aplicación de técnicas geoestadísticas espacio–temporales. Con esta conversión se completó la primera etapa del análisis, dejando la base de datos estructurada adecuadamente para los pasos posteriores del modelamiento espacio–temporal.

Posteriormente, se construyó el diagrama de dispersión entre la temperatura medida a 150 cm de profundidad, ubicada en el eje Y, y el tiempo expresado en horas desde la fecha de referencia, ubicado en el eje X. Este gráfico permitió visualizar el comportamiento temporal de la variable y evidenció la presencia de cuatro agrupamientos o acumulaciones diferenciadas de puntos. Dichos grupos surgen debido a que las mediciones no se realizaron de manera continua a lo largo del tiempo —aunque sí fueron continuas en el espacio—, lo cual generó intervalos sin registro entre los distintos periodos de muestreo.

Posteriormente, se procedió a construir un modelo de regresión lineal multivariado, en el cual la variable dependiente correspondió al valor de la temperatura medida a 150 cm de profundidad. Como variables independientes se emplearon las coordenadas espaciales x y y, así como el tiempo expresado en horas desde la fecha de referencia, tal como se describió previamente en el proceso de conversión temporal.

A partir del ajuste del modelo y del análisis del resumen estadístico de los coeficientes, se observó que los tres parámetros estimados (los coeficientes asociados a x,y y al tiempo) resultaron estadísticamente significativos. Esto indica que las variables independientes explican de manera adecuada la variabilidad observada en la variable dependiente.

Con base en estos resultados, se decidió utilizar los residuales del modelo de regresión como insumo para el análisis geoestadístico posterior. El propósito de este procedimiento fue eliminar la tendencia global presente en los datos y conservar únicamente la estructura de dependencia espacial y temporal de carácter local, necesaria para la construcción del modelo de kriging espacio–temporal.

Después de esta etapa, se llevó a cabo un análisis exploratorio independiente para determinar, a partir de los rezagos temporales, cuáles podrían considerarse adecuados para el conjunto de datos. De forma análoga, se realizó el análisis univariado de los rezagos espaciales con el fin de establecer su comportamiento y definir posteriormente un esquema de discretización apropiado. A partir de estos resultados preliminares se concluyó que era necesario fijar un intervalo de búsqueda para cada rezago, tanto espacial como temporal; en el caso espacial se estableció un ancho de 100 metros alrededor del valor central del rezago, mientras que para el componente temporal se definió un intervalo de una hora alrededor del rezago temporal correspondiente.

Una vez definidos estos intervalos, se procedió con el conteo del número de parejas de observaciones que caían simultáneamente dentro de cada ventana espacial y temporal. Para ello se adoptó un criterio unívoco: una pareja se contabilizaba únicamente si cumplía simultáneamente con el intervalo espacial asignado y con el intervalo temporal correspondiente. Este procedimiento se aplicó de manera exhaustiva para todas las combinaciones posibles de rezagos espacio–temporales.

Posteriormente, se determinó el valor representativo tanto del rezago espacial como del rezago temporal, dado que estos valores eran necesarios para la construcción del diagrama tridimensional del semivariograma espacio–temporal y para la correcta segmentación de los ejes. En el caso espacial, dicho valor se obtuvo mediante el cálculo de la media aritmética del número de parejas que satisfacían simultáneamente las condiciones del intervalo espacial y temporal. De forma equivalente, para el rezago temporal se identificaron las parejas que cumplían con ambos criterios y se calcularon los tiempos de adquisición (en horas) asociados a esas observaciones, a partir de los cuales se obtuvo también el promedio aritmético como valor representativo.

Con estos valores definidos, se procedió al cálculo de la semivarianza para cada conjunto de parejas en cada rezago espacio–temporal. Para ello se empleó el estimador clásico por momentos del semivariograma, el cual consiste en obtener la diferencia entre los valores observados de la variable de interés, elevarla al cuadrado y dividirla entre dos veces el número total de parejas válidas para ese rezago. De esta manera se obtuvo la estimación empírica de la semivarianza correspondiente a cada combinación de rezago espacial y temporal.

Finalmente, con el propósito de representar simultáneamente los tres componentes del semivariograma —los rezagos espaciales, los rezagos temporales y la semivarianza— se construyó una visualización tridimensional utilizando la biblioteca plotly, lo que permitió obtener el gráfico mostrado previamente. En esta representación se aprecia que la mayor parte de los puntos presentan valores de semivarianza reducidos, lo cual se evidencia en la concentración de tonalidades oscuras, asociadas a semivarianzas inferiores a cinco unidades. Este comportamiento indica una elevada similitud entre las observaciones a distancias cortas y medias, reflejando una marcada autocorrelación espacial y temporal en la variable analizada.

En términos geoestadísticos, esta estructura sugiere que la variable aleatoria —la temperatura medida a 150 centímetros de profundidad— exhibe un cambio gradual tanto en el dominio espacial como en el temporal, lo que implica la presencia de un rango relativamente amplio. Dicho de otra manera, la dependencia espacio–temporal persiste a lo largo de distancias y lapsos considerables, lo que coincide con el patrón observado en la gráfica tridimensional.

## convergence: 0

## objective: 72456.9

Ahora bien, para ajustar un modelo teórico de semivariograma espaciotemporal al semivariograma empírico obtenido previamente, fue necesario realizar un análisis preliminar orientado a determinar cuál familia de modelos resultaba más adecuada. En geoestadística espaciotemporal existen dos enfoques generales para caracterizar la dependencia conjunta: los modelos separables y los modelos no separables. Por lo tanto, el primer paso consistió en establecer cuál de estas dos estructuras describía de manera más fiel el comportamiento observado en los datos.

Con este propósito se generaron distintos cortes del semivariograma empírico. En una primera etapa se fijaron varios rezagos espaciales y, para cada uno de ellos, se examinó la evolución del semivariograma temporal. Posteriormente, se aplicó el procedimiento inverso: se seleccionaron diferentes rezagos temporales y se evaluó cómo variaba el semivariograma espacial bajo cada uno de esos cortes. Este análisis permitió comparar la morfología de las curvas en ambos dominios y detectar patrones consistentes o divergentes entre sí.

Los resultados mostraron que las curvas de semivarianza obtenidas mediante cortes espaciales y las derivadas de cortes temporales no presentaban un comportamiento homogéneo. En particular, se observó que la forma del semivariograma temporal cambiaba a medida que aumentaba la distancia espacial, y de manera análoga la estructura del semivariograma espacial se modificaba cuando los rezagos temporales eran mayores. Esta interacción entre dominios indica que la dependencia no puede descomponerse como un producto o suma de componentes puramente espaciales y puramente temporales, lo cual descarta la hipótesis de separabilidad.

En consecuencia, se concluyó que la estructura de dependencia espaciotemporal del proceso solo puede ser descrita adecuadamente mediante un modelo no separable, ya que las componentes espacial y temporal se encuentran acopladas. Por esta razón se optó por emplear el modelo no separable propuesto por Cressie y Huang, cuya formulación permite capturar explícitamente esta interacción entre escalas espaciales y temporales.

Finalmente, se construyó la superficie ajustada al semivariograma espaciotemporal empírico empleando el modelo no separable seleccionado. La representación tridimensional resultante permitió visualizar la estructura conjunta de la dependencia espacio–tiempo. En dicha superficie se aprecia una inclinación dirigida hacia la región central del dominio, mientras que los valores de semivarianza aumentan progresivamente en dirección a rezagos espaciales y temporales más grandes. Este comportamiento es consistente con un proceso que exhibe mayor disimilitud entre observaciones conforme crece la separación tanto en el espacio como en el tiempo, lo cual respalda la validez del modelo teórico adoptado para describir la estructura de correlación espaciotemporal.

## tibble [350 × 11] (S3: tbl_df/tbl/data.frame)
##  $ ID                : num [1:350] 1 2 3 4 5 6 7 8 9 10 ...
##  $ X_Magna           : num [1:350] 927039 927786 928562 929306 928059 ...
##  $ Y_Magna           : num [1:350] 594302 595082 595686 596361 598005 ...
##  $ T_150cm           : num [1:350] 16.9 16 14.9 13.8 15 15 14.7 14.7 14.1 15.1 ...
##  $ Fecha             : Date[1:350], format: "2015-01-31" "2015-01-31" ...
##  $ Hora              : POSIXct[1:350], format: "1899-12-31 14:03:00" "1899-12-31 14:42:00" ...
##  $ Horas_desde_ref   : num [1:350] 0 0 0 0 48 48 48 48 48 72 ...
##  $ Hora_sola         : 'hms' num [1:350] 14:03:00 14:42:00 16:10:00 17:23:00 ...
##   ..- attr(*, "units")= chr "secs"
##  $ Hora_frac         : num [1:350] 14.1 14.7 16.2 17.4 10.3 ...
##  $ Tiempo_total_horas: num [1:350] 14.1 14.7 16.2 17.4 58.3 ...
##  $ residual          : Named num [1:350] 2.288 1.381 0.247 -0.872 0.702 ...
##   ..- attr(*, "names")= chr [1:350] "1" "2" "3" "4" ...

Finalmente, se procedió a realizar la etapa de predicción mediante Kriging Simple Espaciotemporal, con el objetivo de estimar los valores de la temperatura en un dominio simultáneamente continuo en el espacio y en el tiempo. Para ello se definió una malla de predicción estructurada, conformada por 701 nodos espaciales, y para cada uno de estos puntos se generaron 301 predicciones temporales, lo que permitió obtener una superficie espaciotemporal densa y adecuadamente discretizada. Esta estructura de grilla facilitó la evaluación completa del modelo no separable ajustado, al proyectar el comportamiento de la variable en un dominio espacio–tiempo continuo y garantizar una caracterización detallada de la dinámica térmica a 150 cm de profundidad.

Posteriormente, con base en el cubo de predicciones obtenido mediante el Kriging Simple Espaciotemporal —el cual representa la evolución conjunta de la temperatura en el dominio espacio–tiempo— se decidió generar cortes temporales para visualizar la distribución espacial de la variable en fechas específicas. Para ello se extrajeron mapas correspondientes al 4 de febrero de 2015, 4 de junio de 2015 y 4 de octubre de 2016, con el propósito de analizar la dinámica térmica superficial en distintos momentos del periodo de estudio.

Los mapas obtenidos muestran un comportamiento espacial consistente: las temperaturas más bajas se concentran sistemáticamente en el sector noroccidental del área geotérmica del volcán Azufral. Además, en las predicciones correspondientes a junio de 2015 y octubre de 2016 se observa una disminución térmica adicional en la zona central del dominio, coincidiendo espacialmente con la ubicación del edificio volcánico. Esta distribución es coherente desde el punto de vista geofísico, ya que las mediciones se realizaron a una profundidad somera de aproximadamente 1.5 metros, en una región de páramo localizada sobre un edificio volcánico cuya altitud alcanza los 4.000 metros sobre el nivel del mar. En consecuencia, la presencia de valores térmicos más bajos en dichas zonas concuerda con la topografía y con el régimen climático característico de áreas de alta montaña, validando así la consistencia física de los resultados generados por el modelo espaciotemporal.

11. Repetir todos los pasos 1.1-1.6 para datos una variable aleatoria espacial funcional con argumento en una dimensión, esto es, curvas. Este ejercicio hacerlo preferiblemente con la misma variable usadas para el kriging, pero con series largas de datos. (de tiempo, frecuencia, profundidad, etc, según sea el caso).

Para iniciar la construcción de los datos funcionales asociados a la temperatura medida a 150 centímetros de profundidad, se partió de un principio fundamental del análisis de datos funcionales: las funciones deben derivarse de observaciones discretas que presenten continuidad en el dominio, ya sea espacial, temporal o ambos. Sin embargo, las mediciones originales obtenidas en el volcán Azufral consistían únicamente en una observación de temperatura por cada punto espacial, lo cual impedía, en su forma original, la aplicación directa de técnicas funcionales.

Esta limitación se resolvió utilizando los resultados generados por el Kriging Simple Espaciotemporal, ya que a partir de este modelo se obtuvo un cubo de predicciones con continuidad explícita en el espacio y en el tiempo. En particular, se dispuso de 701 localizaciones espaciales en la malla de predicción, y para cada una de estas posiciones se generaron 301 estimaciones temporales de la temperatura. De esta manera, cada punto espacial quedó asociado a una trayectoria temporal continua y suficientemente densa, cumpliendo así con las condiciones necesarias para la representación y posterior análisis funcional de la variable aleatoria temperatura.

Es así como se decidió, después de extraer esta información correspondiente con la interpolación optimidad del Krig en espacio temporal, se decidió utilizar como función base el método de BSP-Line y con este se construyeron 700 curvas continuas que representan la variabilidad de la temperatura en cada uno de los puntos espaciales observados.

Luego de definir las funciones representativas de cada serie temporal en las estaciones, se implementó el método de λ–Kriging para realizar la interpolación funcional mediante kriging ordinario. En este caso, las predicciones corresponden directamente a funciones completas, ya que el método trabaja sobre el espacio funcional sin necesidad de reducir la dimensión.

Posteriormente, se aplicó el segundo enfoque de interpolación funcional: el Score–Kriging. Este método requiere un paso adicional, dado que consiste en proyectar las funciones originales sobre una base reducida obtenida mediante Análisis de Componentes Principales Funcionales (FPCA). Tras estimar los scores principales, se aplicó nuevamente kriging ordinario sobre estos scores y, finalmente, se reconstruyeron las funciones predichas mediante la combinación lineal de las componentes principales retenidas.

Finalmente, con las funciones predichas por ambos métodos (λ–Kriging y Score–Kriging), se realizaron cortes temporales específicos para analizar la evolución espacio–temporal. Se generaron mapas interpolados para las siguientes fechas seleccionadas: 1 de febrero de 2015, 1 de agosto de 2015, 1 de septiembre de 2015, 1 de octubre de 2015, 1 de noviembre de 2015 y 1 de septiembre de 2016.

## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]

## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]
## [using ordinary kriging]

Los resultados obtenidos a partir de los cortes temporales realizados sobre las superficies interpoladas —tanto para el método λ-Kriging como para el método Score-Kriging— permiten observar patrones espacio-temporales consistentes en la evolución de la temperatura superficial.

Para las fechas cercanas al inicio del periodo de muestreo, es decir, alrededor del 31 de enero de 2015, se evidencia una tendencia leve hacia temperaturas más bajas en la zona noroccidental del área de estudio. En contraste, en la zona central aparecen núcleos bien definidos de temperaturas elevadas, con valores cercanos a 18 °C, lo cual sugiere condiciones térmicas localmente más cálidas en este sector durante el periodo inicial.

Sin embargo, en los cortes correspondientes a agosto, septiembre, octubre y noviembre de 2015, así como en el corte de septiembre de 2016, emerge un patrón diferente. En el intervalo comprendido aproximadamente entre las coordenadas 610 000 y 928 000 (eje X), se observa de manera consistente una franja de temperaturas bajas, con valores por debajo de los 10 °C. Este comportamiento coincide espacialmente con la ubicación del Volcán Azufral, situado a una altitud cercana a los 4000 m s. n. m., lo cual respalda la interpretación de que la topografía y la elevación influyen directamente en la distribución espacial de la temperatura reconstruida mediante kriging funcional.

Adicionalmente, en las fechas como octubre de 2015, se identifican zonas más cálidas hacia el oriente del área, alrededor de la coordenada 940 000 en el eje X. Estas áreas corresponden a las partes más bajas del sistema geotérmico del Azufral, abarcando poblaciones como Túquerres, Ospina y Sapuyes. En este sector, las superficies interpoladas muestran temperaturas por encima de los 14 °C, lo que contrasta con los valores significativamente menores observados en la parte alta del volcán.

12. Conclusiones

Una conclusión fundamental de este trabajo es que los distintos enfoques de interpolación empleados —Kriging simple espacial, Kriging espacio–temporal y Kriging funcional— convergen en la identificación de un patrón térmico espacialmente coherente alrededor del volcán Azufral. En el caso del Kriging simple espacial aplicado a la temperatura medida a 150 cm de profundidad, se observa de manera consistente que el intervalo comprendido aproximadamente entre las coordenadas 608 000 y 930 000 m presenta temperaturas predichas más bajas, incluso en sectores sin observaciones directas. Este comportamiento coincide con la ubicación y elevación del volcán Azufral, cuya altitud cercana a los 4000 m s. n. m. explica físicamente la reducción térmica en esta zona.De forma concordante, los resultados derivados del Kriging espacio–temporal muestran que, para fechas específicas como el 4 de agosto de 2016, las temperaturas a 150 cm de profundidad no solo son más bajas en la zona del volcán, sino que además cubren un área más extensa y presentan un desplazamiento hacia el noroccidente, reforzando la estabilidad temporal del patrón identificado.Finalmente, las superficies obtenidas mediante el Kriging funcional —tanto con el método λ como con el método Score— confirman la presencia de temperaturas significativamente más bajas en el sector del volcán durante fechas como agosto y septiembre de 2015. Simultáneamente, en los tres enfoques de interpolación se observa que las temperaturas más altas tienden a ubicarse de manera sistemática hacia el sureste del área de estudio, donde las altitudes son menores y se localizan poblaciones como Túquerres, Ospina y Sapuyes.En conjunto, la coherencia entre los tres métodos respalda sólidamente la validez de las interpolaciones realizadas y evidencia la influencia conjunta de la topografía, la elevación y la dinámica espacio–temporal en la distribución térmica del sistema geotérmico del volcán Azufral.

Otra conclusión relevante es que el Kriging espacio–temporal no solo permite predecir la evolución de la temperatura en cada punto del espacio a lo largo del tiempo, sino que también constituye un insumo fundamental para la construcción de datos funcionales. Al generar series temporales completas en cada ubicación espacial —coherentes con la dinámica geotérmica y la geología del área de estudio— este método proporciona la continuidad temporal necesaria para modelar la variable como una función. De este modo, el Kriging espacio–temporal se convierte en la base que permite aplicar posteriormente técnicas de análisis funcional, asegurando que las funciones obtenidas tengan un comportamiento físico plausible y consistente con el contexto geológico del volcán Azufral.

En el proceso de construcción de un modelo de Kriging espacial, resulta fundamental realizar un análisis previo detallado que permita comprender adecuadamente la estructura de autocorrelación espacial de la variable que se desea interpolar. Para ello, no es suficiente elaborar un único semivariograma empírico omnidireccional; es indispensable complementar este análisis con semivariogramas direccionales y con la evaluación explícita de la isotropía o anisotropía presente en los datos. La comparación entre semivariogramas calculados en distintas direcciones facilita identificar patrones espaciales preferenciales, cambios en el rango o variaciones en la estructura de dependencia según el ángulo de análisis. De esta manera, es posible determinar si el fenómeno se comporta de manera homogénea en todas las direcciones (isotropía) o si, por el contrario, presenta direcciones de continuidad espacial más marcadas (anisotropía). Incorporar esta información en la selección y ajuste del modelo teórico asegura una representación más precisa y robusta de la dependencia espacial de la variable aleatoria estudiada, lo cual se traduce en predicciones de Kriging más confiables.

En la construcción de modelos de regresión lineal simple o múltiple, es fundamental realizar un análisis exhaustivo de los supuestos asociados al modelo, especialmente cuando los residuales serán utilizados posteriormente en un análisis geoestadístico. Este proceso de validación incluye verificar que: La esperanza de los residuales sea cero, es decir, que no exista un sesgo sistemático. Los residuales sean independientes e idénticamente distribuidos (i.i.d.), lo cual garantiza que no persistan estructuras no explicadas por el modelo.Exista homocedasticidad, es decir, que la varianza de los residuales permanezca constante en todo el dominio.Los residuales sigan una distribución normal, requisito importante para diversos procedimientos inferenciales.La evaluación de estos supuestos permite determinar si el modelo de regresión utilizado representa adecuadamente la variabilidad de la variable dependiente. Solo cuando los supuestos se cumplen es posible afirmar que la tendencia global ha sido correctamente removida y que los residuales resultantes constituyen un insumo adecuado para modelar la estructura de dependencia espacial de la variable aleatoria mediante técnicas geoestadísticas.

13. Referencias

• Olade, Icel y Geotérmica Italiana S.R.L. (1982). Estudio de Reconocimiento de los Recursos Geotérmicos de la República de Colombia. Bogotá: CONTECOL. https://catalogo.sg c.gov.co/cgi-bin/koha/opac-detail.pl?biblionumber=12940
• CHEC. (1983). Estudios de prefactibilidad de desarrollo geotérmico en las áreas del Complejo Volcánico Nevado del Ruiz