Clase 6. EL PROGRAMA SACUDIDA PARA LA ESTIMACION DE LA PELIGROSIDAD SISMICA (primera parte)

Consideraciones generales

La versión presente del programa SACUDIDA fue confeccionada atendiendo a todos los aspectos descritos en la clase anterior y permite realizar los cálculos, tanto con el modelo de Riznichenko, modificado para evitar el uso de mapas de actividad sísmica, como con el modelo de tratamiento probabilístico de la atenuación de McGuire-Bender. Otros modelos estadísticos no han sido incluidos.

La evaluación de las integrales se realiza mediante su transformación en sumatorias:

(1)

(2)

donde los índices i,j corresponden a las coordenadas horizontales y el k a la coordenada vertical, mientras que el término f_ijk corresponde al factor de normalización espacial de los gráficos magnitud-frecuencia. Esto implica una subdivisión del espacio V de integración de la forma siguiente:

• en las coordenadas x,y se crea un reticulado con espaciamiento Dx,Dy

• en la vertical se toman las profundidades de cada zona fuente de terremotos como parámetros fijos y se preparan mapas con la representación en superficie, sobre este reticulado, de las zonas fuente mencionadas, con la condición de que sean no interceptantes. Cuando dos o más zonas se interceptan o superponen deben representarse en mapas diferentes. Si todas las zonas fuente se pueden representar en un mapa, k=1, si no, se preparan tantos mapas como sean necesarios para que se cumpla esa condición, con el correspondiente incremento del valor de k.

El factor f_ijk representa la contribución de la celda ijk de dimensiones Dx,Dy, a una profundidad h_k, en la sacudibilidad de la región de interés (en superficie), la cual también esta dividida en cuadrículas de dimensiones Dx,Dy, a cuyos centros se le asignan los valores calculados. De hecho esto no es mas que la normalización espacial de los gráficos magnitud-frecuencia que mencionamos en la segunda clase. El paso de integral a sumatoria implica una aproximación, ya que tanto los datos iniciales como los resultados finales están referidos a cuadrículas de dimensiones Dx,Dy. Con respecto a los datos iniciales, las zonas fuente, al ser superpuestas al reticulado pueden ser aproximadas correctamente solo en el caso de Dx,Dy muy pequeñas (lo que conspira contra el tiempo de cómputo). Con respecto a los resultados, por limitaciones técnicas, el proceso de interpolación para pasar de los valores en el reticulado a mapas de isolíneas debe ser realizado por el usuario. Para ello el programa da 3 posibilidades:

• uso directo de los valores de los períodos de recurrencia/retorno para intensidades dadas, los cuales pueden ser impresos en forma de tablas

• uso de una impresión simbólica (un símbolo por cuadrícula) de estos valores agrupados en 10 intervalos, lo que permite que la interpolación se realice manualmente; donde el error en el trazado de las isolíneas debe ser del orden de (Dx+Dy)/8

• procesamiento ulterior de los valores de los períodos de recurrencia (o retorno) con programas especializados de interpolación y trazado de isolíneas; para ello existe la opción de grabarlos en un fichero.

Además de esta opción de obtención de resultados en forma reticular, el programa permite la realización de los cálculos para puntos aislados a partir de sus coordenadas geográficas, las cuales en general son no coincidentes con el centro de las cuadrículas de la red. En este caso el error de aproximación se limita sólo al mapa de zonas fuente, lo que garantiza resultados más precisos.

Los factores f_ijk se calculan, para cada zona fuente, de la siguiente forma. Sobre un mapa de zonas fuente superpuesto al reticulado en (x,y) se determina que cuadrículas son cubiertas por cada zona fuente y en que grado, se calcula igualmente el área equivalente (en número de cuadrículas) cubierta por la zona fuente, y los factores de normalización se obtienen como el cociente del grado de cubrimiento de una cuadrícula entre el número equivalente de cuadrículas cubierto por la zona. Se debe señalar que los factores de normalización calculados de esta forma consideran equiprobable la generación de terremotos en toda la extensión de una zona fuente. En caso de existir lugares preferentes para la generación de terremotos, ésta se debe subdividir y asociar a cada subzona un nuevo valor del parámetro a de la forma:

(3)

donde p_i es el peso relativo de la ocurrencia de terremotos en la subzona i. Esta fórmula se emplea igualmente para subdividir una zona fuente de gran espesor (en profundidad) en varias subzonas con profundidades diferentes.

Debe señalarse que la forma implícita de división en ejes coordenados corresponde a grados de latitud y longitud en una proyección Mercator, con valores positivos para la longitud oeste y la latitud norte. No obstante, existe la opción de fijar el sistema de coordenadas que se desee, pero debe tenerse en cuenta que el eje x debe ser orientado en sentido contrario a lo habitual en coordenadas cartesianas.

Estructura y principios de trabajo del programa

Está compuesto del programa principal y numerosas subrutinas. Un grupo grande de subrutinas se encarga de la determinación de los valores iniciales del parámetro que se van a utilizar y de la distancia máxima de búsqueda en función de las magnitudes máximas posibles en la región de estudio. Este parámetro puede ser intensidad, aceleración, velocidad, desplazamiento, etc., cada uno de los cuales posee diferentes fórmulas de atenuación. Otro grupo de subrutinas se encarga del modelo de isosistas, del cálculo de las frecuencias acumulativas y de los estimados probabilísticos. También existe un conjunto de subrutinas para impresión de resultados, incluyendo algunas que grafican simbólicamente datos y resultados. La versión que ustedes reciben esta confeccionada en Fortran 77, y compilada con el MS Fortran 5.2

El principio de trabajo del programa es el siguiente. Se parte de la región dividida en cuadrículas elementales y de un conjunto de mapas de zonas fuente no interceptantes. Después de leer los datos para el primer mapa de zonas fuente se fija el valor de x (intensidad, aceleración, velocidad, etc.). Para cada cuadrícula de la zona de interés se calcula la sacudibilidad. El proceso consiste en buscar todas las zonas fuente que pueden influenciar la cuadrícula, tomando en consideración el modelo de isosistas y determinar M_x . Si M_x<M_max+DM_max (ó M_x<M_c+DM_c/2) se calcula la frecuencia acumulativa N (M_x ) o la esperanza matemática E[N (x)] (en función del modelo seleccionado). El proceso se repite para todos los valores de x. En caso de más de un mapa de zonas fuente, se leen los datos de los siguientes y se repite todo el proceso. Al final se calculan los períodos de recurrencia resultantes y los estimados probabilísticos para todos los valores de x.

Existe la opción de obtener impresiones parciales de los resultados para cada mapa de zonas fuente. El número de estos últimos no está limitado; por ello, si se desea analizar la influencia particular de cada zona fuente se pueden preparar mapas independientes para cada una de ellas, o en mapas independientes aquellas cuyo efecto se desea analizar por separado, colocando las demás en el número mínimo de mapas que cumpla la condición de no intersección de zonas.

Los datos iniciales son, de forma general:

• parámetros de control

• parámetros de cada una de las zonas fuente.

• mapas de zonas fuente.

Los resultados son:

• impresión opcional de datos iniciales.

• períodos de recurrencia/retorno de las sacudidas para diferentes valores del parámetro x. Impresión opcional tabulada o simbólica.

• probabilidad de que no sean igualados o excedidos diferentes valores del parámetro x para tiempos de espera de 1, 10, 20, 50, 70, 100, 200, 500 y 1000 años. Impresión opcional tabulada o simbólica.

Fórmulas de atenuación implícitas y modelo de isosistas

El programa tiene incluidas diferentes fórmulas de atenuación. Cada una de ellas es caracterizada, además de por los coeficientes correspondientes, por 3 arreglos unidimensionales de 8 elementos, que son:

• valores del parámetro a utilizar en los cálculos (ordenados de menor a mayor)

• símbolos (4 caracteres) que corresponden a como se desee imprimir ese valor en los distintos títulos que escribe el programa

• distancias a que un terremoto con M = 8 provoca esos valores

Los valores implícitos considerados son:

a_v - 2%, 4%, 10%, 20%, 30%, 40%, 60% y 80% de g

a_h - 4%, 10%, 20%, 30%, 40%, 60%, 80% y 100% de g

v_v - 2.5, 5. 10., 15., 20. 30., 40. y 60. cm/seg

v_h - 2.5, 5. 10., 15., 25., 40., 60. y 90. cm/seg

d_v - 1.5, 3., 5., 7., 10., 13., 16. y 20 cm

d_h - 3., 5., 7., 10., 13., 16., 20. y 25 cm

I - III, IV, V, VI, VII, VIII, IX y X grados MSK

- III-IV, IV-V, V-VI, VI-VII, VII-VIII, VIII-IX, IX-X y X-XI ídem

En caso de no desearse utilizar estos valores implícitos deben suministrarse con los datos los referidos arreglos. Igualmente ocurre con las fórmulas de atenuación.

Para el caso de la intensidad sísmica sólo está incluida la fórmula (23) con los coeficientes determinados por Fedotov y Shumilina (1971) para Kamchatka:

(4)

Como parámetros ingenieriles se utilizaron aceleración vertical (a_v) y horizontal (a_h), velocidad vertical (v_v) y horizontal (v_h), así como desplazamiento vertical (d_v) y horizontal (d_h). Fueron incluidos los diversos casos de fórmulas de atenuación que se presentaron en la clase 3. Aparte de estas fórmulas implícitas pueden ser utilizadas otras, siempre y cuando sean del tipo

(5)

(6)

y que representen tanto uno de los parámetros ya considerados como otro diferente.

El programa tiene implícito el uso del modelo de isosistas elípticas presentado en la clase 3. En el caso de isosistas circulares se usarán los valores A/B=1.0 y h =0.

Preparación de los datos iniciales.

El paso más difícil en la utilización de un programa de estimación de la peligrosidad sísmica es la preparación de los datos iniciales. En ella hay dos aspectos claves: los mapas de zonas fuente y los parámetros de cada una de ellas. Cada zona fuente está caracterizada por un conjunto amplio de parámetros que corresponden al modelo de isosistas, al régimen sísmico y al factor de normalización espacial. El programa requiere que cada zona fuente sea homogénea en todo ese espectro de variables, por lo que si se presentan inhomogeneidades en algunas de ellas, las zonas deben subdividirse hasta que la condición se cumpla. Así, un mapa con n_i zonas, en el proceso de preparación de datos se transformará en otro con n_f zonas, donde n_f³n_i. Los parámetros de las zonas fuente deben incluirse en los datos de forma independiente para cada una de las zubzonas en que quede dividido el mapa correspondiente .

En la versión presente del programa SACUDIDA estos mapas deben ser cifrados manualmente, con valor 0 para la cuadrícula donde no exista zona fuente, y cuando esta exista, con el número de orden para dicha zona. Esto obliga a numerar las zonas fuente en orden consecutivo a partir de 1, lo cual es un aspecto muy importante, ya que el programa siempre se referirá a cualquier zona de acuerdo a ese orden. La numeración se debe realizar no sobre la base de las n_i originales, sino a partir de las n_f obtenidas como resultado de la subdivisión anteriormente mencionada.

Un aspecto que casi siempre conduce a la subdivisión de las zonas es la determinación de los factores de normalización. Dicho proceso no está automatizado por lo que es conveniente explicarlo por partes. Veamos un ejemplo sencillo, una zona fuente circular sobre un reticulado bastante espaciado, de forma tal que la misma quede circunscrita en cuatro cuadrículas (Fig. 1). Como en ellas está exactamente la misma área cubierta, el factor de normalización es de 0.25 para todas y no es necesario subdividir la zona (n_i=n_f=1), por lo que todas las cuadrículas que ella cubre se numerarán con el número 1. Debemos señalar que esta es una aproximación grosera, ya que como resultado del proceso de transformación de integral a sumatoria una zona circular se ha convertido en cuadrada.

Si esa misma zona se superpone a un reticulado el doble más denso (Fig. 2), tendremos tres casos de cubrimiento de cuadrículas: cuatro centrales cubiertas totalmente (f=1/N₀), ocho a ambos lados de los ejes principales con un cubrimiento x (f=x/N₀), y cuatro en los vértices con un cubrimiento menor y (f=y/N₀). N₀ es el número equivalente de cuadrículas cubiertas por la zona fuente, y su valor viene dado por N₀=4x1+8.x+4.y. En este caso la zona fuente se ha desdoblado en 3 subzonas (n_i=1, n_f=3), las cuales se deben simbolizar con los números de orden de las 3 subzonas resultantes (1-3). Con respecto al caso anterior la aproximación es mucho mejor, ya que aunque se mantiene la conversión de zona circular en cuadrada, el peso de cada cuadrícula disminuye en dirección a los vértices. Una densificación adicional del reticulado mejoraría aún más la aproximación.

En la fig. 3 se presenta un caso de 2 zonas fuente. Para ambos casos el cubrimiento de las cuadrículas es de dos tipos, por lo que el problema se convierte en uno de 4 zonas fuente; para la primera tenemos dos casos de cubrimiento: dos cuadriculas con nivel x (simbolizadas con 1) y 2 con nivel y (simbolizadas con 2); para la segunda zona se tienen 2 cuadrículas con cubrimiento v (simbolizadas con 3) y dos con w (simbolizadas con 4). Nótese que para cumplir el requisito de numeración consecutiva, tan pronto terminan de numerarse las subzonas con diferente nivel de cubrimiento de la zona de la izquierda se continúa la numeración con las subzonas correspondientes de la zona de la derecha. Los factores de normalización serán, [f=(x/N₀,y/N₀), N₀=2.x+4.y] para la primera y [f=(v/N₀,w/N₀), N₀=2.v+2.w] para la segunda.

Mapas más complejos se preparan con el mismo principio, pero surge una dificultad relacionada con el hecho de que raras veces las zonas fuente son tan regulares y ocupan una disposición espacial respecto al reticulado tan simétrica como las de los ejemplos anteriores, por lo que es necesario realizar aproximaciones. Resulta útil, cuando el reticulado es bastante fino, fijar arbitrariamente varios niveles de cubrimiento de cuadrículas, y determinar para todas las zonas fuente del mapa los factores de normalización correspondientes. De tal forma, si se fijan los niveles de cubrimiento {v, w, x, y, z, 1}, los factores de normalización serán:

(7)

Bajo este principio, un cubrimiento real u se aproxima a uno de los 5 prefijados y se simplifica mucho el trabajo. Por ejemplo, se obtienen buenos resultados considerando sólo los 4 niveles de cubrimiento (0.25, 0.50, 0.75, 1.0) a los que corresponden los factores de normalización (0.25/N₀, 0.5/N₀, 0.75/N₀, 1.0/N₀).

Ya se vio que por los factores de normalización las zonas fuente se desdoblan, igualmente ocurre con respecto al régimen sísmico y al modelo de isosistas. Veamos el caso del régimen sísmico. Sea una zona fuente lineal (fig 4), en la cual los terremotos ocurren preferentemente en su región central(BC).

Es obvio, que si se construye el gráfico magnitud-frecuencia para la zona como un todo, se distribuirá artificialmente la actividad sísmica y se perderá la especificidad señalada. La variante óptima sería la construcción de estos gráficos para cada uno de los sectores (AB, BC, CD), en cuyo caso esta zona se desdoblaría en 3 subzonas. No obstante, no siempre es posible realizar esto, ya que no se dispone de una estadística suficiente, por lo que deben tomarse soluciones de compromiso. Es muy común construir un gráfico para varias subzonas, considerar el parámetro b como constante para todas y determinar el parámetro a para cada una de ellas de la forma:

(8)

donde "total" corresponde a todo el conjunto de datos, "subzona" a la subzona específica, y M₁ es el centro del primer intervalo representativo de magnitud en el gráfico magnitud-frecuencia acumulativo. En el ejemplo anterior tendríamos:

(9)

donde

La solución del problema sera casuística, en algunos casos se construirán gráficos para zonas fuente individuales y en otros se construirán para varias zonas agrupadas.

Con respecto a la M_max, por lo general se estima de forma independiente para cada zona o zubzona; ya sea por métodos geológicos, sismoestadísticos, de reconocimiento de patrones, etc. La profundidad es igualmente otro de los parámetros del régimen sísmico que se determina mediante el análisis de la sismicidad. Se debe tener presente que en un mapa de zonas fuente, cada una de ellas debe tener una profundidad fija. Si ocurre el caso de que una zona fuente de gran espesor deba ser subdividida, esto implica la preparación de varios mapas para cumplir el requisito de que las zonas sean no intersectantes; en este caso la intersección en el plano corresponde a una superposición en profundidad. Esta situación se presenta también cuando bajo un mismo punto se encuentran 2 o más zonas fuente diferentes, como el caso de las regiones de subducción, donde aparecen, tanto la sismicidad superficial asociada a la cadena volcánica, como la profunda asociada a la capa subducente.

El modelo de isosistas es otra de las fuentes de variación de los parámetros de las zonas fuente. En principio puede variar cualquiera de los parámetros que lo definen, aunque en el programa no está prevista la variación de los coeficientes de la fórmula (5) de una zona a otra. Por ejemplo en la zona representada en la fig 5 tenemos que el ángulo con respecto a la dirección EW es diferente en los sectores AB, BC y CD, lo que necesariamente (si A/B desigual de 0) implica un desdoblamineto de la zona original en 3 subzonas.

Además, pudiera presentarse el caso de que en el sector AA' la excentricidad de las elipses (razón A/B de los semiejes) fuera diferente que en el A'B; en tal caso debe subdividirse adicionalmente la zona. También pudiera ocurrir que en el sector CC' la dirección de medición del radio efectivo fuera el semieje mayor A, mientras que para el sector C'D fuera el radio medio , en cuyo caso también debe subdividirse dicha zona.

En resumen se tendrán tantas subdivisiones de las zonas fuente como casos en que se diferencie al menos uno de los parámetros

• régimen sísmico (a, b, h, M_max)

• modelo de isosistas (D_e, A/B, h, n)

• factor de normalización

• modelo del terremoto característico [M_c, DM_c, n(M_c)]

Todas estas subdivisiones deben numerarse en orden consecutivo, y dichos números son los que se deben asignar a las cuadrículas del mapa (o mapas) de zonas fuente que se prepare.

Finalmente, la determinación de los parámetros dM y s_I que reflejan el error estimado de las formulas de atenuación (aspecto no tratado en el curso) no siempre resulta fácil. En el caso más simple de campo de isolíneas circulares, el parámetro s_I se calcula cuando por análisis de regresión se determinan los coeficientes de la fórmula de atenuación, siendo posible estimar el dM a partir de ese mismo material, lo cual le da un sentido físico mayor que el correspondiente al simple semiintervalo de agrupamiento de los datos para la confección de los gráficos magnitud- frecuencia. Sin embargo, cuando se usa el modelo de isosistas elípticas el proceso es mucho más complejo dadas las características del mismo. El procedimiento para la solución de este problema debe partir del análisis del ajuste del modelo con los parámetros medios para una zona fuente con respecto a cada uno de los terremotos de dicha zona, independientemente de los valores específicos con que fueran ajustados originalmente, proceso que debe realizarse para cada zona fuente. En este análisis deben considerarse varios factores, como son: puntos incluidos en isosistas de grado diferente al reportado, rango de variación de la magnitud que admite cada caso para el resto de los parámetros fijos, así como grado de confiabilidad de los valores de intensidad reportados. No obstante, como este proceso es muy complejo puede resultar conveniente tomar soluciones de compromiso; en el epígrafe siguiente se presenta un análisis práctico de como influyen los valores seleccionados de estos parámetros sobre los estimados de peligrosidad sísmica.

Notas:

1. Los pies de figura se refieren a “zonas OT” en lugar de “zonas fuente”. Esta sigla corresponde a zonas “de origen de terremotos”, nomenclatura en uso en la antigua escuela soviética, equivalente a la más extendida de zonas fuente, que es se ha empleado a lo largo de este curso.

2. Los alumnos deben descargar también el programa de cálculo. Viene en el fichero sacudida.zip que contiene el programa (.exe), unas instrucciones de como preparar los datos (.doc) y un juego de datos de prueba (.dat). Se recomienda que lo prueben con ese juego de datos a ver si trabaja bien. A partir de ahí, pueden modificar algunos parámetros (siguiendo las instrucciones mencionadas) y familiarizarse con su uso. En la segunda parte de esta clase se realizarán algunos ejercicios con el programa.