1. Introducción

El 16 de abril del año 2016 a las 23:58 GMT, 18:59 hora local, se produjo un terremoto de 17.4 km de profundidad en la provincia de Manabí, Ecuador. El epicentro del sismo ocurrió aproximadamente a 27 km de las ciudades de Pedernales y Muisne, con una magnitud registrada de Mw = 7.8. Posterior al terremoto se contabilizaron 661 personas fallecidas, 6274 heridos y 28,678 desplazados sin contar con el aumento de estas cifras debido a las réplicas (^{Goretti et al., 2017}).

El catálogo sísmico del Instituto Geofisico de la Escuela Politecnica Nacional (IG) está constituido por los sismos con registro instrumental en la zona de Ecuador desde el inicio del siglo XX. Se suman además los datos de los sismos históricos desde el año 1587. Para los sismos de magnitud menor a M = 4.5 y mayor a M = 3.5 el periodo de completitud comienza en 1997 (^{Beauval et al., 2013}). La completitud para las magnitudes hasta M = 3.0 se alcanzó en el 2008 (^{Alvarado et al., 2018}).

Desde la instalación del sistema SeisComP3 el IG reporta magnitudes de momento sísmico (^{Alvarado et al., 2018}). Los eventos previos a la instalación de este sistema han sido homogenizados para tener todo el catálogo en magnitud de momento (^{Beauval et al., 2013}).

El catálogo del IG contiene la totalidad de los sismos localizados incluyendo precursores y réplicas. Se registran todos los eventos de las diversas fuentes sísmicas que existen en Ecuador: corticales, interfase o intraslab. En el Ecuador las fuentes más profundas de sismos en el slab están a 250 km. El terremoto de Manabí forma parte del catálogo de sismos homogeneizado. Los 14,474 datos del catálogo tomados desde julio de 1901 hasta noviembre de 2021 (Figura 1) muestran la función de potencia entre la frecuencia y la magnitud según la ley Gutenberg-Richter (GR), Figura (3 a y 3 b), que permite calcular la probabilidad de amenaza sísmica (^{Lee, 2002}).

Donde N es el número total de terremotos que exceden la magnitud m normalizado por el lapso de tiempo total cubierto por el catálogo de terremotos. Los valores a y b son constantes empíricas positivas específicas, que varían dependiendo de cada región (da Silva, 2021). En la Figura 1 se muestra la distribución de eventos sísmicos registrados en la zona alrededor de Manabí, influenciada por el terremoto de Pedernales. Se destacan los sismos de magnitud mayor a M = 5, así como la ubicación de la provincia de Manabí y los cantones Pedernales, Jama y Sucre. Además, se indica la posición del terremoto ocurrido el 16 de abril de 2016.

Figura 1 Eventos sísmicos registrados en la zona alrededor de Manabí cerca de los cantones Pedernales, Jama y Sucre, influenciada por el terremoto de Pedernales de 2016. 

Figura 2 Bosquejo de la placa de Nazca que entra en subducción debajo de la placa Sudamericana. El movimiento de la placa está bloqueado por asperezas de diversos tamaños que interaccionan todas entre sí siguiendo la dinámica de un sistema complejo. La ruptura de las asperezas origina los sismos. 

Figura 3 a. Se muestra el ajuste modificado GR utilizando los datos del catálogo sísmico del IG para la zona alrededor de Manabí. Figura b. Histograma de frecuencias que muestra la magnitud del momento sísmico vs el número de eventos para la zona alrededor de Manabí. 

La ley GR es la base para el cálculo de la amenaza sísmica en una región. El aporte de la complejidad mediante a esta ley se observa en las magnitudes menores a M = 4 ya que para magnitudes mayores el ajuste se aproxima bastante a la clásica recta de la ley GR. Para el cálculo del peligro sísmico en Ecuador se toman los eventos con magnitud mayor a M = 4 (^{Beauval, 2013}). Pero estructuras civiles más sensibles requieren incluir magnitudes menores en el análisis y esto es posible mediante el ajuste basado en los sistemas complejos.

En nuestro artículo no llegamos a un análisis probabilista de peligro sísmico. Sin embargo, existen evidencias que la aplicación de la estadística no extensiva al análisis de peligrosidad en la región de Teherán en Irán, incrementan los valores de la peligrosidad sísmica respecto a los cálculos clásicos (^{Motaghed et al., 2023}).

La ley GR únicamente toma en cuenta la frecuencia de los sismos y la magnitud de momento sísmico obtenida de los equipos de medición, más no el mecanismo físico que genera los terremotos.

El modelo de interacción fragmento aspereza planteado por ^{Sotolongo-Costa y Posadas (2004)} modifica esta ley GR al considerar el movimiento de los bloques que conforman la falla y las asperezas presentes entre los bloques como un sistema complejo (^{Araujo, 2017}).

En la Figura 2 hemos representado la dinámica compleja de la placa de Nazca. La placa oceánica entra en subducción debajo de la placa Sudamericana y su movimiento se bloquea debido a numerosas asperezas de varios tamaños. Estas asperezas forman un sistema complejo debido a que todas interaccionan entre sí ya que la ruptura de alguna de ellas afecta en mayor o menor grado a los esfuerzos tectónicos sobre las demás. Cuando una aspereza se rompe se genera un sismo y la magnitud del terremoto está condicionada por las dimensiones lineales de la aspereza. Hay muchas más asperezas pequeñas que grandes por los que la distribución de las magnitudes de los sismos va a seguir una ley de potencia. La hipótesis fundamental de la aplicación de la estadística no extensiva de Tsallis a los terremotos es que las asperezas forman un sistema complejo (^{Sotolongo-Costa & Posadas, 2004}).

Hay que advertir que, si bien hemos representado en la Figura 2 las asperezas de la zona de interfase entre las placas, el sistema complejo que estudiamos incluye también todas las fallas en la corteza que también tienen su propia estructura compleja. El análisis estadístico no extensivo en zonas de subducción permite por lo tanto tener en cuenta toda la variedad de mecanismos mediante los cuales se producen los movimientos sísmicos.

Este tipo de sistemas se caracterizan por presentar correlaciones de largo rango tanto en el espacio como en el tiempo, procesos no Markovianos, es decir que tienen memoria a largo plazo, ruidos aditivos y multiplicativos en sus ecuaciones mesoscópicas tipo Langevin, caos débil con exponentes de Lyapunov que tienden a cero, geometría multifractal, largo rango de interacción en sistemas de muchos cuerpos, ecuaciones mesoscópicas de Fokker-Planck no lineales y/o no homogéneas (^{Soria, 2018}). Por esta razón la introducción de la mecánica estadística no extensiva se muestra como una solución ad hoc para brindar un sustento teórico a la ley GR.

Otra de las diferencias entre el ajuste GR y el de sistemas complejos es que para efectuar el primero se requiere eliminar las réplicas ya que en principio se toman todos los sismos como independientes. Desde los sistemas complejos en cambio la hipótesis de partida es precisamente la contraria: todos los sismos forman un único sistema y están todos interrelacionados entre sí. No tendría por lo tanto sentido hacer un análisis mediante la entropía de Tsallis de los sismos sin sus réplicas.

En el caso de existir una brecha sísmica en la zona de estudio se manifestaría por una falta de actividad sísmica magnitud mayor a 5 o 6. Los sismos de magnitud menor a 5 o 4 sin embargo, si registran en el catálogo y por lo tanto incorporan a la zona de brecha sísmica al sistema complejo analizado.

Hay que recordar que para escoger nuestra área de estudio hemos partido del hecho que toda la región forma un único sistema tectónico (^{Nocquet et al., 2017}). Un análisis a una menor escala solo considerando los sismos de interfase entre la placa de Nazca y la placa Sudamericana, no es posible por el momento ya que si bien hay avances en el conocimiento de la interfase entre las placas usando tomografía sísmica (Lynner et al., 2019), no se tienen todavía estudios de este tipo al norte de la posición del terremoto de Pedernales.

De acuerdo con el formalismo no extensivo, se sustituye la entropía de Boltzmann por la entropía no aditiva S _q propuesta por ^{Tsallis en (Tsallis, 2009)}, dada por

Donde W es el número total de configuraciones, p en este caso es la probabilidad de encontrar un fragmento de superficie relativa σ, y k es la constante de Boltzmann. La función logaritmo q introducida en la entropía no aditiva es una modificación de la función logaritmo en la que el valor de q modifica la curva de la función (^{Tsallis, 2009}).

De acuerdo con (^{Vallianatos et al., 2016}), la entropía de Tsallis es una generalización de la entropía de Boltzmann para casos en los que q = 1.

A partir de la entropía de Tsallis, en su forma continua representada por

Donde k es una constante positiva que en el caso de la termodinámica estadística es igual a la constante de Boltzman y para la teoría de la información es igual a 1. La variable continua σ representa a todos los microestados que recorre el sistema cada uno de los cuales tiene probabilidad p(σ). El parámetro q representa lo complejidad del sistema y es igual a uno si la complejidad es nula.

A partir de la ecuación (5) se obtiene la función de distribución de probabilidad utilizando el principio de máxima entropía sujeta a 2 restricciones (^{Sotolongo-Costa & Posadas, 2004}). La primera restricción corresponde a la normalización de la función de probabilidad expresada mediante:

En la que p(σ) representa la función de probabilidad y la suma continua de las probabilidades de todos los micro estados es igual a 1. La segunda restricción se relaciona con el valor medio-q de la función de probabilidad, representado de forma continua mediante:

En este caso se observa una modificación al valor medio para casos en los que q es diferente de 1. Tomando en cuenta la proporcionalidad entre la energía sísmica liberada y la magnitud del momento sísmico se obtiene una versión modificada para la ley GR en la que también se presenta el parámetro q (^{Sotolongo-Costa & Posadas, 2004}):

El mecanismo físico de la generación de terremotos aún no ha sido desarrollado por completo por lo que aún no es posible hacer predicciones certeras (^{Vallianatos et al., 2016}). El parámetro de complejidad q obtenido a partir del modelo fragmento aspereza aplicado al estudio de terremotos representa una medida de complejidad y estabilidad sísmica (^{Posadas & Sotolongo-Costa, 2023}), y ha sido calculado para algunas regiones sísmicas. El catálogo de Ecuador presenta un valor de q =1.74 para sismos de profundidad intermedia, entre 75 km y 300 km, q =1.72 para los sismos superficiales menor a 75 km y q = 1.71 para la totalidad del catálogo de sismos del Ecuador (^{Araujo, 2017}). Un resultado análogo se obtiene para el catálogo de sismos de Panamá, con un valor de q = 1.69 (^{Samudio et al., 2015}). Para el catálogo sísmico de Japón se realizó una evaluación de la evolución temporal de este parámetro en fechas cercanas al terremoto de Kobe de 1995 y se interpretó el valor de q como un indicador de estabilidad de esta zona (^{Papadakis et al., 2015}).

El modelo fragmento aspereza se aplica con éxito al catálogo sísmico del Ecuador donde la sismicidad se produce mediante la deformación cortical ocasionada por la subducción oblicua. Esta subducción se produce con una velocidad aproximada de 58 mm/yr (^{Gailler et al., 2007}). La dinámica entre las 2 placas produce sismos de subducción y terremotos corticales de magnitud mayores a 5 y hasta M = 8.5 (^{Gailler et al., 2007}).

En este estudio se examina el posible cambio del parámetro q en relación al tiempo para los precursores y réplicas del terremoto de Pedernales, y evaluar si existen cambios abruptos al presentarse sismos de gran magnitud M~7. El análisis se lleva a cabo para los periodos anterior y posterior al terremoto del 16 de abril de 2016, sin embargo, para el cálculo se utiliza el catálogo sísmico total delimitado para la zona de influencia cercana a la provincia de Manabí.

2. Datos y cálculo del conjunto de parámetros

Los datos empleados para este estudio corresponden al subconjunto del catálogo de sismos correspondiente a la zona de mayor influencia para el terremoto de Manabí. El análisis se realiza para 3609 sismos corticales superficiales de profundidad máxima de 75 km. La ubicación de este enjambre de sismos se encuentra entre latitud máxima de 1.25º, latitud mínima -1.5º, longitud máxima -79.25º y longitud mínima de -82.0º, de acuerdo a la ruptura de un segmento de 100 km y deslizamiento sísmico expuesto en (^{Nocquet et al., 2017}) . La Figura 1, la cual muestra los sismos con magnitud mayor o igual a M = 5 en la zona de estudio. No podemos analizar el cambio del parámetro q únicamente para estos sismos debido a que el parámetro solo tiene sentido desde la física estadística y necesitaríamos un gran número de estos eventos para alcanzar un análisis significativo.

El parámetro q se obtiene a partir de la ecuación (8) mediante la resolución del problema inverso por mínimos cuadrados no lineales desde una perspectiva Bayesiana. Esto se aplica para casos en los que se obtiene información sobre los valores de algunas cantidades medidas, y se busca emplear una relación teórica para obtener información sobre los valores del conjunto de parámetros (^{Tarantola & Valette, 1982}).

El problema directo d=g(p), donde d son los datos medidos y p es un vector que contiene los parámetros se resuelve a través de la convergencia del siguiente algoritmo (^{Tarantola & Valette, 1982}).

En donde Cp0p0 y Cd0d0 son matrices de covarianza, G es la matriz de derivadas parciales y G ^T es su matriz transpuesta.

La separación geográfica del catálogo sísmico se realizó en el lenguaje Fortran 90 y la resolución del problema inverso se modeló a través del programa Matlab. Para obtener la evolución temporal del coeficiente de complejidad, se calculó q para cada mes utilizando todas las magnitudes previas en una secuencia que se extiende hasta completar el catálogo. Por último, también se obtuvo el ajuste con la relación GR modificada y el histograma de frecuencia-magnitud para el conjunto de datos correspondientes a la zona delimitada de Manabí.

3. Resultados y Discusión

Se obtuvo el ajuste de los valores de magnitud vs el logaritmo del número de eventos y el histograma de frecuencia-magnitud utilizando la ecuación (8) para todo el catálogo de la provincia de Manabí desde el año 1900. Analizando el catálogo total se observa un rango amplio de los eventos sísmicos, con una mayor presencia de los sismos de menor magnitud M<4.

El ancho de las clases en los histogramas se calcula como la diferencia entre las magnitudes máxima y mínima divididas para un valor constante. Nuestro conjunto de datos tiene un número de eventos del orden de los miles por lo que el valor de la constante es 15.

El valor de q obtenido para la sección del catálogo total correspondiente a Manabí es q = 1.71. El valor de q calculado utilizando los precursores, que constan de 2674 datos de magnitud de momento sísmico, desde el primer dato registrado el 31 de enero de 1906 hasta el terremoto de Manabí, es de q = 1.7584. Por otro lado, en el caso de las réplicas, que comprenden 935 datos posteriores al terremoto del 16 de abril, se observa un valor de q = 1.687. También se obtuvo el ajuste GR modificado para los precursores y las réplicas y las gráficas frecuencia vs magnitud.

Tomando un intervalo de día a día de forma acumulativa para el cálculo de cada valor de q (Figura 4) se observan 3 cambios de pendiente importantes, uno en el 2016 correspondiente a la fecha del terremoto, otro cambio menor cercano al año 2010, y otro en el año 2005. Los cambios de pendiente se toman desde el año 1995 debido a que la insuficiencia de datos previos a esta fecha no permite la realización del estudio estadístico. La información digital de la RENSIG está disponible a partir de 1988, pero la cantidad de datos comienza a ser significativa desde el año 1995. Aun así, para el cálculo del parámetro q se toman todos los datos previos desde el año 1900.

Figura 4 Evolución temporal de q desde 01-Enero-1995 hasta 01-Enero-2018 y los importantes cabios de pendiente. 

Figura 5 Relación q vs tiempo al rededor del 16 de abril de 2016, fecha del terremoto de M = 7.8. 

En la gráfica (Figura 4) se observa un cambio de pendiente el 16 de diciembre de 1996 asociado a un sismo de magnitud M = 5.5 mar adentro, cerca de la orilla entre los cantones Jama y Sucre, provincia de Manabí longitud: -80.80º, latitud: -0.12º, otro cambio se produce el 4 de agosto de 1998 debido a un sismo de magnitud de M = 7.1 en las coordenadas longitud: -80.66º, latitud: -0.49º, se muestra también un cambio de pendiente el 28 de septiembre del 2000 por un sismo de magnitud M = 6.4. En enero del 2005 se presenta un cambio de pendiente aún más pronunciado que coincide con un enjambre de sismos de magnitudes entre 5 y 6 cerca de las costas de los cantones Montecristi, Manta, Jipijapa y Puerto López en la provincia de Manabí. El 10 de diciembre del 2007 hay un cambio de pendiente correspondiente a un sismo de magnitud M = 5.1, y el 28 de febrero de 2010 se muestra un ligero cambio de pendiente en relación a un sismo de magnitud M = 5.2. La ubicación de los sismos registrados, así como la fecha de ocurrencia se obtienen a partir del catálogo de sismos. El momento más significativo corresponde al terremoto del 16 de abril de 2016 que muestra un incremento del valor de q asociado al mismo terremoto, tal como se observa en la Figura 4. Este incremento representa una recarga del sistema y un aumento instantáneo de complejidad, sin embargo, luego de este sismo el valor de q continúa descendiendo.

Esta observación de incremento de q coincide con estudios realizados para el terremoto de Kobe de 1995 analizados en (^{Papadakis et al., 2015}) y los datos sísmicos del terremoto de L´Aquila estudiados en (^{Telesca, 2010}) donde se observa un incremento de la complejidad antes de originarse estos sismos.

Las Figuras 6 a, 6 b, 6 c y 6 d muestran la relación frecuencia-magnitud y el ajuste utilizando la ecuación (8) para los sismos previos y posteriores al terremoto de Pedernales. Se observa una reducción en el rango de eventos desde antes hasta después del terremoto del 16 de abril, lo cual concuerda con la disminución de q, indicando una simplificación del sistema. El valor de q en la Figura 4 decrece debido al aumento de estabilidad sísmica que gana el sistema, sin embargo el valor de q aún sigue siendo relativamente elevado en comparación con los estudios de (^{Papadakis et al., 2015}) donde q inicia con un valor de 1.46 aproximadamente previo y posterior al terremoto y 1.55 durante el evento sísmico, y (^{Telesca, 2010}) el cual es 1.48 previo al terremoto de L´Aquila, alcanza un pico de 1.74 durante el sismo y luego desciende al valor de 1.69. El incremento anómalo del parámetro q también coincide con un estudio realizado en la región de Iquique, al norte de Chile, durante 2014 y asociado al sismo de magnitud M = 8.1, donde se introduce el concepto de "mutabilidad" para entender la información en la secuencia sísmica. La mutabilidad mostró un aumento antes del evento principal, señalando cambios significativos en las configuraciones accesibles en el espacio de estados. La correlación entre la disminución de la entropía de Tsallis y el aumento de la mutabilidad sugiere una relación compleja entre la acumulación de energía y su liberación durante los eventos sísmicos. Este hallazgo respalda la idea de que la entropía de Tsallis responde a condiciones en el espacio real, mientras que la mutabilidad está influenciada por los estados accesibles en el espacio de configuración. (^{Pasten et al., 2023}).

Figura 6 a. Ajuste mediante estadística no extensiva para las réplicas del terremoto de Pedernales del 16-abril-2016. Figura 6 b. Histograma de frecuencia-magnitud para las réplicas del terremoto de Pedernales del 16-abril-2016. 

Figura 6 c. Ajuste mediante estadística no extensiva para los precursores del terremoto de Pedernales del 16-abril-2016. Figura 6 d. Histograma de frecuencia-magnitud para los precursores del terremoto de Pedernales del 16-abril-2016. 

4. Conclusiones

El valor de q para el catálogo de Manabí entre 1900 y 2021 fue de q = 1.71; q = 1.75 y q = 1.68 para los sismos obtenidos como precursores y réplicas respectivamente. Los cambios de los valores de q a lo largo de los años presentan una relación con los eventos sísmicos y las magnitudes.

Es importante destacar que un cambio importante de pendiente se presenta asociado a la ocurrencia de un sismo de gran magnitud M>5, en el caso de Manabí este cambio ocurrió durante el terremoto de Pedernales del 16 de abril de 2016. El cambio de q a través de los años entre 1995 y 2021 muestra un descenso en el valor de q lo que implica una reducción de la complejidad del sistema. La evolución temporal de q presenta una similitud con la evolución estudiada para otros catálogos sísmicos en Japón, Europa y Chile.

La recomendación para futuros análisis es la recolección de una mayor cantidad de mediciones a escala de horas con el fin mejorar el pronóstico sísmico en función del cambio de q.