Introducción

El contacto interplaca en las zonas de subducción constituye una de las principales fuentes de peligro sísmico en el mundo debido que en estas regiones se generan los sismos de mayor magnitud. En este contacto se producen grandes esfuerzos que se acumulan a lo largo del tiempo y eventualmente se liberan durante sismos mayores. Se piensa que estos esfuerzos se acumulan en áreas de alta fricción, llamadas asperezas, en la zona de contacto y eventualmente se liberan (^{Kanamori, 1978}; ^{Lay y Kanamori, 1981}) irradiando la energía principal de las ondas sísmicas registradas en el campo lejano. En conjunto, las asperezas y las zonas débiles que las rodean representan el área de ruptura total (^{Kanamori, 1981}, ^{Lay et al., 1982}) que excita las ondas sísmicas de más largo período y define el momento sísmico del evento.

México se encuentra en una región de gran deformación tectónica debido a la interacción de las placas de Rivera, Cocos y Norteamérica. Esta área incluye la región de Michoacán-Colima-Jalisco donde varios terremotos destructivos han ocurrido en el último siglo. La Figura 1 muestra los sismos de magnitud Mw mayor a 7 que han ocurrido en la región desde el año 1900, basado en los momentos sísmicos calculados para los eventos (Tabla 1). Estos sismos generalmente resultan de fallas inversas en el contacto interplaca y han causado daños significativos tanto tierra adentro como cerca de la costa. Los eventos incluyen el terremoto de Mw 8.1 de Jalisco del 3 de junio de 1932, el sismo interplaca más grande registrado instrumentalmente hasta la fecha en la zona de subducción de México. ^{Singh et al. (1985)} estimaron una longitud de ruptura de 220 km para el evento basado en la distribución de réplicas. El sismo Mw 8.1 del 19 de septiembre de 1985 que resultó en daños y pérdidas sin precedentes en la Ciudad de México (^{Esteva, 1988}) también ocurrió en esta región. Recientemente, el evento del 9 de octubre de 1995 (Mw 8.0) generó daños importantes cerca al área epicentral en Manzanillo, Colima y en el noroeste del estado de Jalisco (^{Juárez-García et al., 1997}).

Figura 1 Sismos mayores a Mw 7 (estrellas) que han ocurrido en la zona de subducción de Michoacán-Colima- Jalisco desde 1900 entre las longitudes de -101.5º a -107.00º. Se muestran los mecanismos focales reportados por el Global Centroid Moment Tensor (www.globalcmt.org) para eventos ocurridos de 1973 a 2003. TMA es la trinchera mesoamericana. 

Se han publicado modelos de deslizamiento para varios de los sismos Mw > 7 ocurridos en la zona de subducción de Michoacán-Colima-Jalisco utilizando registros sísmicos (p. ej. ^{Mendoza y Hartzell, 1988}; ¹⁹⁸⁹; ^{Mendoza, 1993}; ^{Mendoza, 1995}; ^{Courboulex et al., 1997}; ^{Mendoza y Hartzell, 1999}; ^{Yagi et al., 2004}; ^{Santoyo et al. 2005}; ^{Mendoza et al., 2012}; ^{UNAM Seismology Group, 1986}; ²⁰¹⁵; ^{Mendoza y Martínez-López, 2017}). Estos modelos delinean la ruptura del terremoto e identifican áreas de concentración de deslizamiento que se han interpretado como asperezas en el contacto interplaca. Sin embargo, los estudios han utilizado metodologías con diferentes tipos de datos, incluyendo registros telesísmicos y en algunos casos, estaciones locales o regionales. Además, los parámetros utilizados para definir las fallas no han sido siempre consistentes entre sí. Los estudios más antiguos, por ejemplo, no han permitido suficiente flexibilidad en el tiempo que tarda el deslizamiento en llegar al valor final (duración de la dislocación), debido principalmente a la capacidad computacional disponible en el momento. Estos diferentes enfoques dificultan una comparación adecuada entre eventos y una evaluación apropiada del comportamiento de las asperezas.

Los sismos mayores a Mw 7 ocurridos desde el año 1980 en la zona de subducción de Michoacán Colima Jalisco han sido registrados a distancias telesísmicas por instrumentos digitales que proporcionan formas de onda de cuerpo de alta calidad que se pueden examinar para identificar las propiedades de la fuente sísmica. Estos registros digitales se han utilizado en algunos casos para la distribución del deslizamiento cosísmico. En este trabajo proponemos analizar estos datos globales para determinar de manera uniforme la historia de ruptura de cinco sismos recientes mayores a Mw 7 que han ocurrido en la zona de subducción de Michoacán-Colima-Jalisco desde el año 1980. Estos incluyen los sismos de Playa Azul del 25 de octubre de 1981 (Mw 7.2), de Michoacán del 19 de septiembre de 1985 (Mw 8.1), de Zihuatanejo del 21 de septiembre de 1985 (Mw 7.5), de Colima-Jalisco del 9 de octubre de 1995 (Mw 8.0), y de Tecomán del 22 de enero de 2003 (Mw 7.5). Se aplica el mismo método de inversión de la fuente utilizando parámetros similares de la falla con el fin de comparar distribuciones similares del deslizamiento cosísmico para los distintos eventos. Utilizamos el esquema de inversión de falla finita de ^{Hartzell y Heaton (1983}; ¹⁹⁸⁶), modificado por ^{Mendoza y Hartzell (2013)} para identificar directamente el nivel de estabilización y obtener la solución en un solo paso. Recientemente, ^{Mendoza y Martínez-López (2017)} aplicaron esta misma metodología a las ondas P y SH registradas a distancias telesísmicas para el sismo Mw 7.3 del 18 de abril 2014 en Papanoa, Guerrero.

La determinación consistente de los modelos telesísmicos nos permite además examinar las posibles incertidumbres debidas al conocimiento incompleto de los parámetros que se utilizan para definir la falla. Variaciones en estos parámetros afectan los resultados obtenidos (^{Hartzell, 1989}; ^{Hartzell y Langer, 1993}; ^{Beresnev, 2003}; ^{Lay et al., 2010}). Por ejemplo, ^{Hartzell y Langer (1993)} sugieren que se pueden generar resultados erróneos si se utiliza una parametrización de la falla que no permita suficiente variación espacial en la amplitud de deslizamiento, en la duración de la dislocación y en el tiempo de inicio de la ruptura. Errores en la geometría de la falla que se utiliza para modelar la ruptura también contribuyen a la incertidumbre del modelo de deslizamiento determinado. En un estudio telesísmico de los eventos del 9 de octubre 1995 (Mw 8.0) y del 22 de enero 2003 (Mw 7.5), ^{Mendoza et al. (2012)} observaron que los errores en el buzamiento y la profundidad de nucleación afectan más la precisión de la fuente inferida que los errores que se tienen en el rumbo y en el ángulo de deslizamiento. Por otra parte, Hartzell et al. (2013) invirtieron formas de ondas registradas a distancias telesísmicas y regionales para obtener distintos modelos de deslizamiento para el evento Mw 5.8 de 2011 occurrido en Mineral, Virginia, EUA. Ellos aplicaron a cada solución un procedimiento simple donde variaron el alineamiento entre datos observados y teóricos de manera sistematica asumiendo una distribución Gaussiana que les permitió calcular la variabilidad del deslizamiento. Mencionan que esta manera de abordar la variabilidad de los errores podría ser útil en el análisis de las incertidumbres variando otros parámetros del modelo (^{Hartzell et al., 2013}). En este trabajo adoptamos este procedimiento para examinar la incertidumbre en los modelos telesísmicos de deslizamiento de sismos recientes en Michocán-Colima-Jalisco debido a posibles errores en los parámetros de la falla. En nuestro caso, variamos el rumbo, el buzamiento, el ángulo de deslizamiento, la profundidad del hipocentro y la velocidad de ruptura utilizando una distribución uniforme de los parámetros. Los resultados de la redeterminación de patrones de ruptura y sus posibles incertidumbres permiten una evaluación de la ubicación e interacción entre asperezas en esta sección de la zona de subducción de México. Además aportan información consistente de las rupturas sísmicas recientes en esta zona que podrían ser de utilidad en la definición de leyes de escalamiento.

Datos

Para determinar los modelos de deslizamiento se invirtieron las ondas de cuerpo telesísmicas registradas en las estaciones digitales mundiales disponibles del Incorporated Research Institutions for Seismology (IRIS, http://www.iris.edu/). Utilizamos las ondas P registradas entre 25 y 95 grados y las ondas SH registradas entre 35 y 80 grados para minimizar los efectos de difracción del núcleo, la propagación en el manto superior y las triplicaciones del manto. Las ondas P se obtuvieron a partir de los registros verticales, y las ondas SH se recuperaron mediante la rotación de los componentes horizontales, corrigiendo los registros P y SH por la respuesta instrumental para obtener observaciones en desplazamiento. Utilizamos el ancho de banda más amplio posible para las observaciones. Para eventos ocurridos después del año 1990 estos corresponden a datos de banda ancha registrados por la Red Sísmica Global (GSN, por sus siglas en ingles). Antes de 1990, las estaciones digitales registraban señales en varias bandas de frecuencias para producir componentes de período largo, período intermedio y período corto. Para esos eventos se utilizaron los registros de periodo largo e intermedio. Sin embargo, había muy pocas estaciones de registro intermedio al inicio de la decada de los ochenta y se consideraron además registros de periodo corto para el sismo de Playa Azul de 1981. Para los registros de banda ancha y de periodo intermedio, aplicamos un filtro pasa banda en un intervalo de periodo entre 1- 60s y remuestreamos las formas de onda a 0.25s. Los registros de período largo se filtraron entre 10-80s y se remuestrearon a 1s. Los registros de periodo corto registrados para el sismo de Playa Azul de 1981 se filtraron entre 1-5s y se remuestrearon a 0.1s.

Metodología

Se aplicó la metodología de falla finita originalmente desarrollada por ^{Hartzell y Heaton (1983)} que se ha utilizado en el estudio de sismos mayores en la zona de subducción de la costa del Pacífico Mexicano (p. ej. ^{Mendoza y Hartzell, 1988}; ¹⁹⁸⁹; ^{Mendoza, 1993}; ^{Mendoza, 1995}; ^{Mendoza y Hartzell, 1999}; ^{Mendoza y Martinez-Lopez, 2017}). La metodología se basa en una parametrización cinemática de la falla para identificar la distribución del deslizamiento cosísmico que mejor reproduce las formas de onda registradas. En la aplicación del método se identifica un plano de falla con orientación y geometría basadas en el mecanismo focal del sismo. Se subdivide la falla en un número específico de celdas y se fija la ubicación del hipocentro. Posteriormente, se calculan los sismogramas sintéticos para cada celda asumiendo que cada una de estas está compuesta de fuentes puntuales distribuidas uniformemente a lo largo y ancho de sus dimensiones. Se asume que cada una de las fuentes puntuales se dispara cuando el frente de ruptura, que viaja a una velocidad constante a lo largo de la falla desde el hipocentro, llega a ese punto. Las respuestas para cada fuente puntual (funciones de Green) se calculan utilizando una función rectangular (boxcar) de duración fija usando un modelo de corteza basado en los resultados obtenidos por ^{Stolte et al. (1986)} para la región de Michoacán-Guerrero. Las funciones de Green se retardadan por el tiempo de ruptura y se suman para obtener el sismograma sintético de cada celda.

El problema numérico de inversión es completamente lineal y se construye poniendo los sismogramas generados para cada celda de todas las estaciónes una tras otra para formar las columnas de una matriz A de amplitudes sintéticas. El número de columnas entonces corresponde al número de celdas que se consideran en la inversión. Las formas de onda observadas en todas las estaciones similarmente se juntan una tras otra para formar un vector de datos b. Los detalles de la construcción de esta matriz y el vector de datos se puede consultar en el trabajo de ^{Hartzell y Heaton (1983)}. Los datos y los sintéticos forman un sistema sobredeterminado de ecuaciones lineales C_d ^-1 Ax=C_d ^-1 b donde C_d ^-1 es una matriz de covarianza de datos que normaliza cada registro de estación a su amplitud máxima. El vector solución x contiene los deslizamientos requeridos en cada celda para reproducir las observaciones. Para permitir flexibilidad en el tiempo de inicio de ruptura en cada celda, se añaden columnas adicionales a la matriz de coeficientes C_d ^-1 A donde los sintéticos normalizados de cada celda se retrasan por el ancho de la función rectangular que se utilizó para generar las funciones de Green. El número de veces que los sintéticos de cada celda se retrasan y se añaden a la matriz de coeficientes corresponde al número de ventanas de tiempo utilizadas para discretizar la duración de la dislocación sobre la falla. La inversión recupera el deslizamiento en cada celda para cada una de estas ventanas de tiempo. Típicamente el deslizamiento mayor se observa en las ventanas iniciales y a partir de cierto tiempo no hay más contribuciones al deslizamiento. Este tiempo define la duración de dislocación requerida por las observaciones.

El problema inverso se estabiliza mediante la adición de ecuaciones de restricción de la forma λFx= 0 al sistema lineal

donde F₁ representa la diferencia de deslizamiento entre celdas adyacentes, y F₂ es la matriz identidad. La primera restricción impone una transición suave de deslizamiento entre celdas, y la segunda restricción reduce la longitud del vector de solución x, minimizando el momento sísmico total. El parámetro λ controla el compromiso entre aplicar las restricciones y ajustar las observaciones y se estima usando la relación λ=90|a| _avg sugerida por ^{Mendoza y Hartzell (2013)}, donde |a| _avg es la media de los valores absolutos de los elementos de la matriz normalizada de coeficientes C_d ^-1 A. Esta relación se obtuvo a partir de una revisión de la variación de los valores de suavizamiento de varios sismos mayor a Mw 6 con modelos de deslizamiento determinados utilizando ondas P telesísmicas (^{Mendoza y Hartzell, 2013}).

Los modelos de deslizamiento obtenidos en este estudio se determinaron tomando el plano nodal con buzamiento noreste del mecanismo focal reportado por el gCMT Project (http://www.globalcmt.org). Se utilizó el epicentro mejor conocido para cada evento obtenido de estudios previos (Tabla 1). Para los eventos con magnitud Mw entre 7.2 y 7.5 se utilizaron 10 ventanas de tiempo con longitudes de 1 s cada una, permitiendo una dislocación máxima de 10s. Para los eventos con magnitudes Mw 8.0 se utilizaron 10 ventanas de tiempo con una longitud de 2s, permitiendo dislocaciones máximas de 20s. Las dimensiones de las celdas son de 5x5 km y de 10x10 km para sismos de Mw ~7 y de Mw 8, respectivamente. Se limito la profundidad máxima de la falla a 40 km para los eventos Mw 8.0 considerando el límite inferior de la zona de acoplamiento sismogénico determinado por ^{Martínez-López y Mendoza (2015)} para la región de Michoacán-Colima-Jalisco.

Para estimar las incertidumbres debido a posibles errores en los parámetros de entrada se aplicó un procedimiento similar al proceso de muestreo que utilizaron ^{Hartzell et al. (2013)} para analizar el sismo de 2011 en Virginia, EUA. En nuestro caso consideramos errores aleatorios en la geometría de la falla, la profundidad del hipocentro y la velocidad de ruptura. Los valores se varian en un rango de ±5° en el rumbo, ±5° en el buzamiento, ±5° en el ángulo de deslizamiento, ±4 km en la profundidad del hipocentro, y ±0.5 km/s en la velocidad de ruptura y distribuidos de manera uniforme con respecto a los valores de entrada. Para los eventos analizados se observó que tanto los residuales numéricos como los ajustes a los sismogramas son similares para modelos determinados con parámetros de entrada en estos rangos, por lo que consideramos que los rangos representan el mínimo error para cada uno de los parámetros. Se hicieron 300 inversiones independientes para cada evento utilizando las diferentes combinaciones de los parámetros de entrada. Los 300 modelos de deslizamiento obtenidos a partir de este proceso se utilizaron para calcular la desviación estándar del deslizamiento en cada celda.

Modelos de Deslizamiento

Sismo del 25 de octubre de 1981 (Mw 7.2)

Para el sismo de Playa Azul del 25 de octubre de 1981, ^{Mendoza (1993)} obtuvo un modelo de deslizamiento utilizando ondas P telesísmicas con una sola ventana de tiempo de 1s. En este estudio incrementamos el número de ventanas a 10 para permitir flexibilidad en la duración de la dislocación sobre el plano de falla. Utilizamos un conjunto de estaciones similar al que utilizó ^{Mendoza (1993)}, incluyendo además 4 registros de periodo corto y una forma de onda SH. Se invirtieron registros de 65s en longitud disminuyendo la amplitud exponencialmente a los 20s, 35s, y 40s después del inicio de los registros de período corto, período intermedio y período largo, respectivamente, para minimizar los efectos de propagación no relacionados a la fuente. Se utilizó una falla de 80 km por 80 km dividida en 256 celdas de 5 km x 5 km con un rumbo de 287°, un buzamiento de 20° y un ángulo de deslizamiento de 82°. El hipocentro se coloca en el centro de la falla y se utiliza una velocidad de ruptura de 2.6 km/s para calcular los sismogramas sintéticos de cada celda.

Se analizaron varias profundidades hipocentrales tomadas entre los rangos de profundidad reportados por diferentes autores. Estas pruebas indican que una profundidad de 17.5 km es la que mejor ajusta los datos observados. La Figura 2a muestra la distribución de deslizamiento cosísmico obtenida para el evento y la Figura 2b muestra los ajustes de registros para un momento sísmico total de 7.6x10¹⁹ Nm (Mw 7.2). En la zona de ruptura se puede observar una región de alto deslizamiento definida por valores mayores a 45% del pico de 132 cm. Esta zona de alto deslizamiento (A) tiene un deslizamiento promedio de 91 cm y un área de 1050 km². Si consideramos esta área equivalente a un área circular podemos calcular la caída de esfuerzo (∆σ) utilizando la relación ∆σ= (7/16)Mo/r³ desarrollada por ^{Eshelby (1957)}, donde Mo es el momento sísmico del área en dinas-cm y r es su radio en cm. Esto nos da un valor ∆σ de 0.3 MPa. Además, se puede identificar una duración de dislocación de 6s para esta misma zona tomando el promedio de los valores obtenidos para cada celda que se encuentra dentro del polígono de la Figura 2a. La Figura 2c muestra la desviación estándar del deslizamiento observado en cada celda a partir de las 300 inversiones. La mayor incertidumbre (30-40 cm) se observa cerca del hipocentro, donde los valores de deslizamiento en el modelo cosísmico son altos.

Figura 2 a) Modelo de deslizamiento obtenido para el evento del 25 de octubre de 1981 a partir de la inversión de datos telesísmicos de ondas P y SH. El rectángulo corresponde a la zona de ruptura definida con el procedimiento de ^{Somerville et al. (1999)}. El polígono incluye las celdas con valores mayores al 45% del deslizamiento máximo y la estrella representa el hipocentro. b) Ajustes entre datos observados (línea continua) y teóricos (línea discontinua) para un momento sísmico de 7.6 x 10¹⁹ Nm (Mw 7.2). Los numeros a la derecha de cada estación corresponden al pico de amplitud (en micrones) del registro observado. c) Desviación estandar del deslizamiento estimado para celda después de 300 inversiones independientes variando el rumbo, buzamiento, ángulo de deslizamiento, profundidad y velocidad de ruptura. 

Sismo del 19 de septiembre de 1985 (Mw 8.1)

^{Mendoza y Hartzell (1988};¹⁹⁸⁹⁾ obtuvieron un modelo de deslizamiento del evento de Michoacán del 19 de septiembre de 1985 a partir de una inversión conjunta de datos locales y formas de onda P telesísmicas donde permitieron una duración máxima de dislocación de 6s. En este trabajo se redeterminó el modelo de deslizamiento utilizando 10 ventanas de tiempo de 2s cada una, permitiendo así un intervalo máximo de dislocación de 20s. Se invirtieron 12 formas de onda P de período intermedio y período largo, y 3 formas de onda SH de período intermedio. Se utilizó una longitud de registro de 120s con una disminución exponencial de amplitud después de 60s y 70s para las ondas P y SH, respectivamente. Se utilizó un plano de falla con rumbo de 301°, buzamiento de 18°, y ángulo de deslizamiento de 105°. La falla tiene una longitud de 220 km y un ancho de 120 km, dividida en 264 celdas de 10 km x 10 km. La velocidad de ruptura es de 2.6 km/s.

La profundidad hipocentral que mejor ajusta los datos observados es de 18 km. En la Figura 3a se muestra el modelo de deslizamiento para esta profundidad. Se observan dos fuentes principales separadas por aproximadamente 30 km, una con un deslizamiento máximo de ~550 cm cerca al hipocentro, y otra al sureste con un deslizamiento máximo de ~300 cm. La zona A de alto deslizamiento que se localiza cerca del hipocentro tiene un deslizamiento promedio de 356 cm, un área de 1200 km² y una caída de esfuerzo de 0.95 MPa. En la zona B de alto deslizamiento al sureste del hipocentro el deslizamiento promedio es 277 cm sobre un área de 900 km² con una caída de esfuerzo de 0.85 MPa. La duración de dislocación en ambas zonas de máximo deslizamiento es 12s. En la Figura 3b se muestran los ajustes de los registros para un momento sísmico total de 9.8 x 10²⁰ Nm (Mw 8.0). La Figura 3c muestra la desviación estándar del deslizamiento observado en cada celda a partir de las 300 inversiones. Para la zona A los valores de desviación estándar varían entre 40 y 95 cm. Para la zona B los valores varían entre 40 y 75 cm.

Figura 3 a) Modelo de deslizamiento obtenido para el evento del 19 de septiembre de 1985 a partir de la inversión de ondas P y SH telesismicas. El rectángulo corresponde a las zonas de ruptura definida con el procedimiento de ^{Somerville et al. (1999)}. Los polígonos marcan las celdas con deslizamiento mayor al 45% del deslizamiento máximo de cada fuente principal. La estrella representa el hipocentro. b) Los ajustes entre registros observados (línea continua) y teóricos (línea discontinua) para un momento sísmico de 9.8 x 10²⁰ Nm (Mw 8.0). Los numeros a la derecha corresponden a la amplitud máxima (en micrones) del registro observado. c) Desviación estandar del deslizamiento estimado para celda después de 300 inversiones independientes variando el rumbo, buzamiento, ángulo de deslizamiento, profundidad y velocidad de ruptura. 

Sismo del 21 de septiembre de 1985 (Mw 7.5)

^{Mendoza (1993)} determinó un modelo de deslizamiento para el evento de Zihuatanejo del 21 de septiembre de 1985 usando formas de onda P y SH telesísmicas y una sola ventana de tiempo de 1s. En este estudio se redeterminó el modelo telesísmico de deslizamiento utilizando 10 ventanas de tiempo permitiendo una duración máxima de dislocación de 10s. Se invirtieron 13 registros de onda P de período largo e intermedio y 4 formas de onda SH de periodo intermedio. Se utilizó una longitud de registro de 75s con una disminución exponencial en amplitud a 40 y 50s para P y SH, respectivamente. Se uso una falla de 100 km x 100 km dividida en 400 celdas de 5 km x 5 km. Se utilizó un plano de falla con rumbo de 296°, buzamiento de 17°, y ángulo de deslizamiento de 85°. El hipocentro se colocó en el centro de la falla y se utilizó una velocidad de ruptura de 2.6 km/s.

Resultados de las pruebas utilizando diferentes profundidades indican que la profundidad de nucleación que mejor ajusta los datos observados es de 24 km. El modelo obtenido para esta profundidad se muestra en la Figura 4a. En esta se observa una zona de alto deslizamiento con un deslizamiento promedio de 245 cm, un área de 775 km² y una caída de esfuerzo de 0.8 MPa. La duración de dislocación en la zona de alto deslizamiento es de 6s. El momento sísmico total es de 1.7 x 10²⁰ Nm (Mw 7.4), correspondiente a los ajustes que se muestran en la Figura 4b. La Figura 4c muestra la desviación estándar observada en cada celda a partir de las 300 inversiones. En la zona de alto deslizamiento del modelo obtenido se observan desviaciónes estándar entre 40 y 120 cm.

Figura 4 a) Modelo de deslizamiento obtenido para el evento del 21 de septiembre de 1985 a partir de la inversión de ondas P y SH telesísmicas. El rectángulo corresponde a la zona de ruptura definida con el procedimiento de ^{Somerville et al. (1999)}. El polígono incluye las celdas con deslizamiento mayor al 45% del deslizamiento máximo y la estrella representa el hipocentro. b) Ajustes entre los datos observados (línea continua) y teóricos (linea discontinua) para un momento sismico de 1.7 x 10²⁰ Nm (Mw 7.4). Los numeros a la derecha corresponden al pico de amplitud (en micrones) del registro observado. c) Desviación estándar del deslizamiento estimado para celda después de 300 inversiones independientes variando el rumbo, buzamiento, ángulo de deslizamiento, profundidad y velocidad de ruptura. 

Sismo del 9 de octubre de 1995 (Mw 8.0)

^{Mendoza y Hartzell (1999)} determinaron un modelo de deslizamiento para el sismo de Colima-Jalisco del 9 de octubre de 1995 (Mw 8.0) usando formas de onda P telesísmicas de banda ancha y 10 ventanas de tiempo de 2s cada una. Aquí hacemos una inversión similar con el mismo número de ventanas, invirtiendo registros banda ancha incluyendo 31 formas de onda P y 9 formas de onda SH. Se utilizan registros de 120s aplicando una disminución exponencial en amplitud a partir de 65s y 70s para las ondas P y SH, respectivamente. Se tomó la geometría de falla del gCMT con rumbo de 302°, buzamiento de 9° y ángulo de deslizamiento de 92°. La longitud y ancho de la falla son 220 km y 100 km, respectivamente. La falla fue dividida en 220 celdas de 10 km x 10 km. La velocidad de ruptura utilizada es de 2.8 km/s.

Resultados de la inversión usando diferentes profundidades hipocentrales indican que la profundidad que mejor ajusta los datos observados es de 15 km. El modelo correspondiente se muestra en la Figura 5a. Se observan tres fuentes principales, una cerca del hipocentro (A) con un pico de deslizamiento de ~175 cm, una segunda fuente (B) 50 km al noroeste con un deslizamiento máximo de ~310 cm y una tercer fuente (C) más al noroeste con un deslizamiento máximo de ~210 cm. Las fuentes B y C se encuentran separadas por ~20 km. La zona A tiene un deslizamiento promedio de 160 cm, un área de 900 km² y una caída de esfuerzo de 0.5 MPa. La zona B tiene un deslizamiento promedio de 175 cm, una área de 1500 km² y una caída de esfuerzo de 0.4 MPa. La tercera fuente (C) tiene un deslizamiento promedio de 180 cm, una área de 600 km² y una caída de esfuerzo de 0.7 MPa. En la Figura 5b se muestran los ajustes a los registros observados y corresponden a un momento sísmico de 6.0 x 10²⁰ Nm (Mw 7.8). La Figura 5c muestra la desviación estándar del deslizamiento observado en cada celda a partir de las 300 inversiones. Para la fuente A se observa una desviación estándar de 30 a 60 cm. En la zona B los valores de desviación estándar varían entre 30 y 80 cm. El valor más alto de desviación estándar se observa en la celda de máximo deslizamiento de la Figura 5a. En la zona C los valores de desviación estándar se encuentran entre 55 y 75 cm.

Figura 5 a) Modelo de deslizamiento obtenido para el evento del 9 de octubre de 1995 a partir de la inversión de ondas P y SH telesismicas. El rectángulo corresponde a la zona de ruptura definida con el procedimiento de ^{Somerville et al. (1999)}. Los polígonos marcan las celdas con deslizamiento mayor al 45% del deslizamiento máximo de cada fuente principal. La estrella representa el hipocentro. b) Ajustes entre las formas de onda observadas (línea continua) y las teóricas (línea discontinua) para un momento sísmico de 6.0 x 10²⁰ Nm (Mw 7.8). Los numeros a la derecha corresponden al pico de amplitud (en micrones) del registro observado. c) Desviación estandar del deslizamiento estimado para cada celda despues de 300 inversiones independientes variando el rumbo, buzamiento, angulo de deslizamiento, profundidad y velocidad de ruptura. 

Sismo del 22 de enero de 2003 (Mw 7.5)

^{Yagi et al. (2004)} obtuvo un modelo de ruptura para el evento de Tecomán del 22 de enero de 2003 a partir de una inversión de datos locales y ondas P y SH telesísmicas. Ellos utilizaron 13 ventanas de tiempo de 1s cada una aplicando una metodología similar a la que se utiliza en este estudio. Aquí utilizamos 10 ventanas de tiempo de 1s cada una para invertir 31 registros banda ancha de onda P y 11 de onda SH. Se invirtió una longitud de registro de 80s con una disminución exponencial de amplitud después de 40s y 48s para ondas P y SH, respectivamente. Se utilizaron dimensiones de la falla de 100 km de largo y 80 km de ancho, dividida en 320 celdas de 5 km x 5 km. El hipocentro se ubica a 30 km del borde sureste y 50 km de la parte superior de la falla. Se tomó un plano de falla con rumbo de 308°, buzamiento de 12° y ángulo de deslizamiento de 110°. Se utilizó una velocidad de ruptura de 3.0 km/s.

Se realizaron pruebas para determinar la profundidad que mejor ajusta los datos observados, identificando una profundidad hipocentral de 22 km. La Figura 6a muestra el modelo de deslizamiento para esta profundidad. Se observan dos fuentes principales, una en la parte superior de la falla (A) con un deslizamiento máximo de 263 cm y otra (B) en la parte inferior con un deslizamiento máximo de 270 cm. En la zona A el deslizamiento promedio es de 205 cm dentro de un área de 250 km² con una caída de esfuerzo de 1.2 MPa. En la zona B el deslizamiento promedio es de 190 cm en un área de 500 km² con una caída de esfuerzo de 0.8 MPa. La duración de dislocación en las zonas A y B es de 6s y 4s, respectivamente. En la Figura 6b se muestran los ajustes que corresponden al momento sísmico estimado de 1.2 x 10²⁰ Nm (Mw 7.3). La Figura 6c muestra la desviación estándar del deslizamiento observado en cada celda a partir de las 300 inversiones. Se observan valores altos de desviación estándar (~90 cm) en la parte superior de la falla. En la zona A se observan desviaciones estándar que varían de 35 a 90 cm. En la zona B se observan desviaciones entre 35 y 100 cm.

Figura 6 a) Modelo de deslizamiento obtenido para el evento del 22 de enero de 2003 a partir de la inversión de ondas P y SH telesismicas. El rectángulo corresponde a la zona de ruptura definida con el procedimiento de ^{Somerville et al. (1999)}. Los polígonos incluyen las celdas con deslizamiento mayor al 45% del deslizamiento máximo observado para cada fuente principal. La estrella representa el hipocentro. b) Ajustes a los registros observados (línea continua) para un momento sismico de 1.2 x 10²⁰ Nm (Mw 7.3). Los numeros a la derecha corresponden a la amplitud máxima (en micrones) del registro observado. c) Desviación estándar del deslizamiento estimado para celda después de las 300 inversiones independientes variando el rumbo, buzamiento, ángulo de deslizamiento, profundidad y velocidad de ruptura. 

Distribución e incertidumbres de las asperezas

Los modelos telesísmicos de deslizamiento obtenidos en este estudio identifican zonas de alto deslizamiento que se interpretan como asperezas en el contacto interplaca. Estas asperezas se muestran en la Figura 7, representadas por los polígonos marcados en las Figuras 2 a 6 para cada evento. Se muestran además zonas de ruptura obtenidas utilizando el procedimiento de ^{Somerville et al. (1999)}. Este procedimiento descarta filas y/o columnas de los modelos que no contribuyen de manera significativa al momento sísmico total y sirve para delinear de manera general la ubicación de la ruptura en la zona de subducción. En general se observa que las rupturas de los sismos Mw ≥ 7 recientes (desde 1980) no se empalman, con excepción del evento de Playa Azul 1981 que se encuentra dentro de la zona de ruptura del sismo de Michoacán 1985. Sin embargo, las dos fuentes principales del evento de Michoacán 1985 excluyen la ruptura del evento de 1981 (p. ej., ^{Astiz et al., 1987}; ^{Mendoza y Hartzell, 1989}; ^{Mendoza, 1993}), consistente con la distribución de asperezas que se observa para los dos eventos en la Figura 7. En esta región también ocurrió el sismo histórico de Mw 7.6 del 7 de junio de 1911 pero la localización de este evento no está bien definida (^{Singh et al., 1981}, ^{Astiz et al., 1987}) y es difícil examinar su relación con los eventos de 1981 y 1985.

Figura 7 Asperezas en la zona de subducción Michoacán-Colima-Jalisco definidas a partir de los polígonos de alto deslizamiento identificados en los modelos de deslizamiento determinados en este estudio. Los rectángulos corresponden a las zonas de ruptura marcadas en las figuras 2 a 6. Se muestran además las áreas de réplicas que se conocen para los sismos del 3 de junio de 1932 (1932-1), del 18 de junio de 1932 (1932-2), del 15 de abril de 1941 y del 30 de enero de 1973 [tomadas de ^{Ramírez-Herrera et al. (2010)} y de ^{Kelleher et al. (1973)} ]. 

La Figura 7 también muestra las áreas de réplicas que se conocen para sismos históricos ocurridos desde el año 1900. Estas áreas parecen coincidir con algunas de las zonas de ruptura de sismos recientes. Por ejemplo, la ruptura del evento del 9 de octubre de 1995 cubre la parte sur del área de réplicas del evento del 3 de junio de 1932 y las réplicas del sismo del 18 de junio de 1932 en la costa de Jalisco. Sin embargo, las réplicas de los eventos de 1932 se encuentran fuera de las asperezas del sismo de 1995. Esto es consistente con la conclusión de ^{Pacheco et al. (1997)} quienes sugieren que el evento de 1995 no corresponde a una repetición de los sismos de 1932 basado en una comparación de los sismogramas registrados para los eventos. En cambio, las réplicas del evento Mw 7.6 de 1973 si se encuentran dentro del área de réplicas determinada por ^{Kelleher et al. (1973)} para el sismo Mw 7.6 de 1941. La distribución de réplicas del sismo de 1941, sin embargo, esta definida con muy pocos epicentros (ver ^{Kelleher et al., 1973}) y es difícil evaluar la relación entre los dos eventos a pesar de que el evento de 1973 cuenta con un modelo preliminar de deslizamiento determinado por ^{Santoyo et al. (2006)}.

En la Figura 7 se puede observar que las asperezas de sismos recientes de Mw > 7 forman un patrón entrelazado de zonas de alto deslizamiento. Por ejemplo, las asperezas de los eventos de 1995 y 2003 en la costa de Colima y también las asperezas de los sismos de 1981 y 1985 en Michoacán se ubican en lugares adyacentes en la zona de contacto interplaca. Esta observación es consistente con los resultados de ^{Mendoza y Martínez López (2017)} quienes observaron que las zonas de máximo deslizamiento del sismo de Papanoa del 2014 (Mw 7.3) en la costa de Guerrero no coinciden con las del evento de Petatlán de 1979 (Mw 7.4). Esta distribución de asperezas adyacentes podría implicar que los espacios entre asperezas de sismos anteriores identifican zonas de mayor deslizamiento en sismos futuros a lo largo del contacto interplaca (^{Mendoza y Martínez-López, 2017}). Otra alternativa es que las asperezas se mantienen fijas en el tiempo y se deslizan en sismos recurrentes, aunque los sismos recientes en la zona de Michoacán-Colima-Jalisco no muestran este comportamiento.

La Tabla 2 muestra las propiedades de las asperezas, representadas por valores promedio calculados de las 300 soluciones para los polígonos definidos anteriormente en las Figuras 2 a 6 para cada evento. Se incluyen valores promedio de 1) D_max, el deslizamiento máximo observado para cada aspereza, 2) D_p, el deslizamiento promedio dentro del polígono, y 3) T_d, la duración promedio de dislocación dentro de este mismo polígono. La Tabla 2 incluye la desviación estándar calculada para cada uno de estos tres parámetros. Las 300 soluciones independientes también se examinaron de manera individual para identificar variaciones en el área de las asperezas. Para cada evento, se identificaron las dimensiones mínimas y máximas de las asperezas y estas se utilizaron para identificar el rango de variación en el área. Basado en estos rangos, se calculo además la variación en la caída de esfuerzo de cada aspereza. Estos rangos en área y caída de esfuerzo se incluyen en la Tabla 2.

Los parámetros presentados en la Tabla 2 sirven para definir la variabilidad que puede haber en las propiedades de las asperezas tomando en cuenta los posibles errores en los parámetros de entrada que se utilizan en el proceso de inversión. Se puede observar que el deslizamiento máximo y promedio de las asperezas generalmente se mantienen en un rango de ~1-4 m con excepción del sismo de Playa Azul 1981 que tiene valores menores y la aspereza A del sismo de Michoacán 1985 que tiene valores mayores. Las duraciones de dislocación son parecidas para eventos de similar magnitud. Para los sismos de Mw ~7.5 los valores son de ~5s y para los sismos Mw 8 el valor es ~12s con una incertidumbre del orden del ancho de la ventana de tiempo, indicando que la duración de la dislocación en las asperezas se escala con el tamaño del sismo.

La variabilidad en el área de las asperezas es similar en general con excepción del sismo de Tecomán 2003 que muestra áreas más pequeñas con menor variabilidad. Este evento del 2003 tiene la mejor cobertura azimutal de estaciones, y la mejor definición en el área de las asperezas podría estar ligada a la distribución de estaciones. Se examinaron las ubicaciónes de las asperezas tomando en cuenta las incertidumbres en el área y se encontro que las asperezas siguen estando en lugares adyacentes a lo largo del contacto interplaca sin un empalme significativo, incluso para eventos que ocurrieron cerca de cada uno (p.ej. los sismos de 1995 y 2003 en Jalisco y los de 1981 y 1985 en Michoacán). Las incertidumbres en las dimensiones de las asperezas sugieren caídas de esfuerzo menores a 2.0 MPa para las asperezas, similares a los valores de 3.0 MPa que se observan en general para sismos de subducción (p.ej. ^{Allmann y Shearer, 2009}).

En algunos casos las soluciones individuales obtenidas de las 300 inversiones muestran más de una aspereza en el mismo lugar. Por ejemplo, en algunas de las soluciones del evento del 9 de octubre de 1995 las asperezas B y C se juntan para formar una sola fuente del tamaño aproximado de las dos asperezas. Similarmente, en algunas soluciones del evento del 22 de enero de 2003 las dos asperezas se juntan para formar una sola que cubre aproximadamente la misma área. Se observa que el momento sísmico total de estas áreas separadas corresponde a la suma del momento sísmico de las asperezas individuales. Estas observaciones reflejan el nivel de incertidumbre que existe en la aplicación de la metodología de falla finita utilizando los datos telesísmicos para estos eventos. Sin embargo, las ubicaciones de las zonas de alto deslizamiento en la zona de ruptura no varían significativamente, indicando que las asperezas de los sismos se pueden recuperar de manera adecuada. Estos resultados además indican que el procedimiento utilizado para muestrear los errores de los parámetros de entrada ofrece una alternativa para determinar las incertidumbres de las asperezas utilizando los datos telesísmicos.

Conclusiones y discusión

Se determinaron de manera sistemática los modelos de deslizamiento de los sismos de magnitud Mw mayor a 7 ocurridos desde 1980 en la zona de subducción de Michoacán-Colima-Jalisco utilizando ondas P y SH registradas a distancias telesísmicas. Se aplicó una metodología de inversión de falla finita utilizando una parametrización consistente que permite una amplia flexibilidad en la duración de dislocación. Los modelos de deslizamiento obtenidos de la inversión se utilizaron para definir las asperezas de los eventos en la zona de subducción. Se observa que las asperezas de los sismos recientes se entrelazan entre si a lo largo del contacto interplaca, consistente con los resultados de ^{Mendoza y Martínez-López (2017)}, quienes encontraron que las zonas de alto deslizamiento del sismo Mw 7.3 de Papanoa en 2014 no coinciden con las áreas de máximo deslizamiento del sismo Mw 7.4 de Petatlán de 1979. ^{Park y Mori (2007)} observaron un patrón similar de empalme limitado entre asperezas para los zona de subducción de Nueva Bretaña en la región del sur-Pacífico.

Nuestros resultados indican además que el contacto interplaca en la región de Michoacán-Colima-Jalisco se ha roto en secciones individuales en los últimos 40 años y quizás los espacios entre asperezas de sismos recientes de Mw ~7 podrían identificar lugares de deslizamiento principal en eventos futuros como lo menciona ^{Mendoza y Martínez López (2017)}. También es posible que las asperezas se mantienen fijas en el tiempo. Sin embargo, los eventos analizados en este estudio no muestran ese comportamiento, quizás debido al corto tiempo de evaluación. Sería importante realizar estudios más completos que incluyan una mayor cantidad de eventos para ver si el comportamiento observado en este trabajo se mantiene para distintos ciclos sísmicos.

En la zona de estudio no se puede descartar la posibilidad de que se rompa todo el segmento en un solo sismo. Esto se ha observado recientemente en otras zonas de subducción como Chile y Japón (^{Delouis et al., 2010}; ^{Simons et al., 2011}) donde la interface de la placa tectónica se había estado rompiendo en eventos de magnitud Mw ~7-8 y no se esperaba un sismo de magnitud 9.0. En México se tienen reportes de un sismo histórico de magnitud Mw 8.6 en Oaxaca el 28 de marzo de 1787, en una región donde solo se han generado sismos de magnitud entre 7.3 a 8.2 en los últimos 100 años (^{Suárez y Albini, 2009}). Es importante entonces seguir investigando el comportamiento de asperezas en esta y otras zonas de subducción para examinar el potencial sísmico real que pueda existir a nivel regional.

Se aplicó además un procedimiento simple donde se corrieron 300 inversiones independientes para identificar las incertidumbres en los modelos telesísmicos de cada evento tomando en cuenta los posibles errores en la geometría de la falla, la profundidad del hipocentro y la velocidad de ruptura. Este ejercicio permitio identificar la variabilidad en las propiedades específicas de las asperezas incluyendo el deslizamiento máximo, el deslizamiento promedio, la duración de dislocación, el área y la caída de esfuerzo. Las incertidumbres que se observan en las dimensiones de las asperezas no son suficientemente amplias para generar un empalme significativo entre asperezas, indicando que los posibles errores en la parametrización de la falla no afecta nuestra conclusión de asperezas entrelazadas para sismos recientes en la zona de subducción de Michoacán-Colima-Jalisco.

Nuestros resultados indican que a grandes rasgos el deslizamiento máximo, el deslizamiento promedio, y las dimensiones de las asperezas de los diferentes eventos son generalmente similares. Esto tendría implicaciones importantes para comprender el proceso de generación de sismos interplaca en zonas de subducción y también para la simulación realista del movimiento del suelo esperado. Sería importante entonces estudiar otros eventos para ver si estas similitudes generales se mantienen. Las duraciones de dislocación de las asperezas están relativamente bien restringidas, con incertidumbres limitadas al ancho de la ventana de tiempo utilizada en el proceso de inversión. Se observan valores de duración para las asperezas que se escalan con el tamaño del evento (~5s y ~12s para sismos de Mw ~7.5 y de Mw 8.0, respectivamente). Las caídas de esfuerzo calculadas para todas las asperezas se mantienen menores a 2.0 MPa, consistente con valores estimados previamente para sismos interplaca en zonas de subducción. Los resultados de la investigación indican que un simple muestreo de los parámetros de entrada ofrece una alternativa para examinar las propiedades e incertidumbres de las asperezas.

Los resultados de este estudio también tienen implicaciones importantes para la estimación de leyes de escalamiento que relacionan los diferentes parámetros de la fuente sísmica. Se han desarrollado leyes de escalamiento para zonas de subducción a partir de modelos de deslizamiento recuperados para eventos de fallamiento inverso (p. ej. ^{Somerville et al., 2002}; ^{Murotani et al., 2008}; ²⁰¹³; ^{Ramirez-Gaytan et al., 2014}). Sin embargo, los modelos de deslizamiento se han determinado con diferentes metodologías, datos y parametrizaciones de la falla lo cual dificulta una definición confiable y consistente de las zonas de alto deslizamiento. Los modelos de deslizamiento que se utilizan además tienen incertidumbres que no se han tomado en cuenta y sería de gran utilidad poder incorporar estas incertidumbres para considerar la posible variabilidad en los parámetros de la fuente. Para la zona de subducción de México, por ejemplo, sería útil examinar las asperezas de todos los eventos de magnitud Mw mayor a 7 que han ocurrido desde 1980 siguiendo el procedimiento desarrollado en este trabajo. Esto resultaría en una definición mas realista de las leyes de escalamiento, lo cual mejoraría el conocimiento actual del peligro sísmico.