1Sismología e imágenes de la Tierra

El 22 de mayo de 1960, a las 19:11 hora universal, el mayor terremoto registrado desde 1900 rompió la zona de subducción frente al sur de Chile, con una magnitud de 9.5.

Sus ondas P llegaron a estaciones situadas al otro lado del planeta unos veinte minutos después; sus ondas superficiales dieron varias vueltas al globo; y durante días los instrumentos de periodo más largo registraron a toda la Tierra vibrando con periodos cercanos a una hora, como una campana golpeada. Cada una de estas señales transporta información sobre las rocas que atravesó. A partir de los tiempos de llegada de millones de estas ondas, registrados en miles de estaciones, los sismólogos han construido no solo el modelo estratificado de la Tierra que se recuerda más abajo, sino también imágenes tridimensionales del manto, en las que las losas que se hunden y las plumas que ascienden aparecen como regiones donde las ondas viajan algo más deprisa o más despacio que la media. Este capítulo muestra cómo: cómo viajan los rayos por una Tierra esférica, cómo las ondas superficiales sondean la corteza y el manto a través de su dispersión, cómo se describen las oscilaciones libres de todo el planeta, y cómo unos pocos miles de tiempos de recorrido se convierten en una imagen mediante la resolución de un problema inverso lineal.

Lo que ya sabes

En un sólido elástico isótropo se propagan dos tipos de ondas de cuerpo: las ondas P, a la velocidad $\alpha = \sqrt{(K + \tfrac43\mu)/\rho}$, y las ondas S, a la velocidad $\beta = \sqrt{\mu/\rho}$, donde $K$ y $\mu$ son los módulos de compresión y de cizalla y $\rho$ la densidad; las ondas S no atraviesan un líquido (el volumen de tercer año de los libros de física, en el capítulo sobre elasticidad de medios continuos). Una onda cuya velocidad de fase depende de su frecuencia es dispersiva; su energía se propaga a la velocidad de grupo $U = \mathrm d\omega/\mathrm dk$ (el volumen de segundo año de los libros de física, sobre ondas). Un cuerpo elástico acotado vibra libremente según un conjunto discreto de modos normales (mismo volumen). La distancia epicentral $\Delta$ de una estación es el ángulo en el centro de la Tierra entre el epicentro y la estación; la curva de tiempos de recorrido da el tiempo de primera llegada en función de $\Delta$, y el núcleo proyecta una zona de sombra para las ondas P directas (libro 2, año 11). El modelo de Tierra de referencia, las discontinuidades sísmicas y la onda de cabeza se estudiaron en el volumen de primer año, en su capítulo sobre la Tierra estratificada. Los mínimos cuadrados y la ley gaussiana de errores son materia de los libros de matemáticas.

1.1 Rayos en una Tierra esférica

Un frente de onda sísmico, lejos de su fuente y en un medio cuyas propiedades varían lentamente a lo largo de una longitud de onda, puede seguirse mediante sus rayos: las curvas perpendiculares en todo punto al frente de onda, a lo largo de las cuales se propaga la energía. En un medio de velocidad uniforme los rayos son rectos; allí donde la velocidad cambia, se curvan según la ley de Snell de la refracción. En la Tierra las velocidades dependen sobre todo de la profundidad, de modo que un primer modelo es esféricamente simétrico: la velocidad $v(r)$ de la onda considerada (P o S) es función únicamente del radio $r$.

Teorema 1.1 (Ley de Snell en una Tierra esférica)

En una Tierra con simetría esférica de velocidad de onda $v(r)$, sea $i$ el ángulo entre un rayo y el radio en el punto de radio $r$. A lo largo de todo el rayo,

$$\frac{r\sin i}{v(r)} = p$$

es constante.

Definición 1.2 (Parámetro del rayo, punto de retorno)

La constante $p = r\sin i/v(r)$ de un rayo en una Tierra con simetría esférica es su parámetro del rayo, en segundos por radián. En la superficie vale $p = R\sin i_0/v_0$, donde $i_0$ es el ángulo de salida y $v_0$ la velocidad en la superficie. El punto donde el rayo es horizontal ($i = 90{}^{\circ}$) es su punto de retorno (de un rayo sísmico), situado en el radio $r_t$ donde

$$\frac{r_t}{v(r_t)} = p .$$

Por debajo de él el rayo vuelve a ascender, simétrico respecto al radio que pasa por el punto de retorno.

Se escribe $\eta(r) = r/v(r)$. Un rayo de parámetro $p$ solo puede existir donde $\eta \ge p$, ya que $\sin i = p/\eta \le 1$. Si $\eta$ disminuye con la profundidad, como ocurre en gran parte del manto, todo rayo retorna a la profundidad en que $\eta$ ha descendido hasta $p$, y los rayos que salen más verticales (menor $p$) retornan a mayor profundidad. Donde $\eta$ aumenta con la profundidad, porque la velocidad disminuye bruscamente con la profundidad (así ocurre en la frontera núcleo–manto, donde la velocidad P desciende de 13.7 a $8.0\,\mathrm{km}\,\mathrm{s}^{-1}$), ningún rayo puede retornar: los rayos se curvan hacia abajo, atraviesan la capa y emergen muy lejos. Un salto de este tipo proyecta una sombra en la superficie.

El parámetro del rayo no es solo una constante del rayo: puede leerse en la curva de tiempos de recorrido.

Teorema 1.3 (Pendiente de la curva de tiempos de recorrido)

Sea $T(\Delta)$ el tiempo de recorrido de un rayo desde una fuente en superficie hasta una estación en superficie a la distancia epicentral $\Delta$ (en radianes). Entonces

$$\frac{\mathrm dT}{\mathrm d\Delta} = p ,$$

el parámetro del rayo que llega a $\Delta$.

El teorema es la forma esférica de la relación $\mathrm dT/\mathrm dx = 1/c$ entre la pendiente de una curva de tiempos de recorrido y la velocidad aparente $c$ a lo largo de la superficie. Las dos funciones $T(p)$ y $\Delta(p)$ pueden escribirse como integrales sobre el perfil de velocidades.

Proposición 1.4 (Tiempo de recorrido y distancia como integrales)

Con $\eta = r/v(r)$ decreciente con la profundidad, el rayo de parámetro $p$ que retorna en $r_t$ alcanza

$$\Delta(p) = 2\int_{r_t}^{R}\frac{p\,\mathrm dr}{r\sqrt{\eta^2 - p^2}}, \qquad T(p) = 2\int_{r_t}^{R}\frac{\eta^2\,\mathrm dr}{r\sqrt{\eta^2 - p^2}},$$

y el tiempo de intercepción $\tau(p) = T - p\Delta$ es $\tau(p) = 2\int_{r_t}^{R}\sqrt{\eta^2 - p^2}\,\mathrm dr/r$.

Estas integrales son la forma en que se calcula una tabla de tiempos de recorrido a partir de un modelo. El script de trazado de rayos de este capítulo realiza la misma suma sobre capas de 2 km del modelo de referencia PREM, y lo comprueba de dos maneras: en una esfera de velocidad uniforme, donde los rayos son cuerdas y $T = 2R\sin(\Delta/2)/v$ exactamente (ejercicio 5), y lo contrasta con el Teorema 1.3, que no utiliza.

La figura anterior muestra lo que predicen los teoremas. El rayo del manto más inclinado que apenas roza el núcleo llega a $98.5{}^{\circ}$ en PREM; el rayo siguiente, más inclinado en una fracción de grado, entra en el núcleo, se refracta hacia el radio y sale por el lado opuesto, más allá de $140{}^{\circ}$. Al representar el tiempo frente a la distancia de cada rayo se obtienen las curvas de tiempos de recorrido que aparecen a continuación.

Los tiempos P calculados a través de PREM son $273.5\,\mathrm{s}$ a $20{}^{\circ}$, $369.6\,\mathrm{s}$ a $30{}^{\circ}$, $607.2\,\mathrm{s}$ a $60{}^{\circ}$ y $779.7\,\mathrm{s}$ a $90{}^{\circ}$: trece minutos para atravesar un cuarto del globo.

Los arribos reales se apartan de esa curva en unos pocos segundos como mucho, y esos pocos segundos son los datos de la tomografía (Sección 1.5).

Definición 1.5 (Tabla de tiempos de recorrido)

Una tabla de tiempos de recorrido proporciona, para cada fase sísmica (P, S, PKP, PcP, …) y cada profundidad del foco, el tiempo de recorrido en función de la distancia epicentral, calculado a partir de un modelo de referencia con simetría esférica cuyas velocidades se ajustaron a su vez para reproducir los tiempos de un gran conjunto de terremotos bien localizados. Sus diferencias respecto de un arribo observado son el punto de partida de la localización de terremotos y de la tomografía.

Las tablas en uso hoy se calculan a partir de modelos de referencia construidos con ese fin. En el modelo ak135, por ejemplo, la corteza tiene dos capas, con velocidades de las ondas P y S de 5.80 y $3.46\,\mathrm{km}\,\mathrm{s}^{-1}$ por encima de $20\,\mathrm{km}$ y 6.50 y $3.85\,\mathrm{km}\,\mathrm{s}^{-1}$ hasta la discontinuidad de Moho, a $35\,\mathrm{km}$, donde el manto empieza con 8.04 y $4.48\,\mathrm{km}\,\mathrm{s}^{-1}$; la frontera núcleo–manto se sitúa a $2891.5\,\mathrm{km}$ y la frontera del núcleo interno a $5153.5\,\mathrm{km}$.

Sus propios autores describen el manto superior como una construcción artificial diseñada para ajustar los tiempos hasta $30{}^{\circ}$: un promedio, no la estructura de ningún lugar concreto.

Proposición 1.6 (Inversión de Herglotz–Wiechert)

Supóngase que $\eta = r/v(r)$ decrece estrictamente con la profundidad en los radios muestreados y que la curva de tiempos de recorrido $T(\Delta)$ es conocida, de modo que $p(\Delta) = \mathrm dT/\mathrm d\Delta$ es conocida. Entonces el rayo de parámetro $p_1$, que llega a la superficie en $\Delta_1$, retorna en el radio

$$r_1 = R\,\exp\!\left(-\frac1\pi\int_0^{\Delta_1} \operatorname{arccosh}\frac{p(\Delta)}{p_1}\,\mathrm d\Delta\right),$$

donde la velocidad es $v(r_1) = r_1/p_1$. La inversión de Herglotz–Wiechert inversión recupera $v(r)$ de este modo, punto a punto, a partir de la curva de tiempos de recorrido.

La fórmula procede de la integral de $\Delta(p)$ en la Proposición 1.4, que es una ecuación integral de Abel para la incógnita $r(\eta)$; su inversión se admite aquí. El ejercicio 12 la comprueba en la esfera uniforme. Dos condiciones la limitan. Se necesita que $\eta$ decrezca con la profundidad: una capa de baja velocidad produce un hueco en la curva de tiempos de recorrido, y la inversión no puede ver el interior de la capa. Y se necesita la curva completa, incluidas las ramas posteriores de cada triplicación, que son difíciles de identificar en un sismograma. Estas limitaciones explican por qué los modelos de referencia se construyen hoy mediante modelización directa y mínimos cuadrados (Sección 1.4).

1.2 Ondas superficiales

Una superficie libre atrapa ondas. Cerca de la superficie de un semiespacio, los movimientos P y S en el plano vertical de propagación pueden combinarse en una onda que viaja a lo largo de la superficie y decae con la profundidad; y donde una capa lenta yace sobre roca más rápida, las ondas de cizalla horizontales quedan atrapadas en la capa por reflexión total.

Definición 1.7 (Onda de Rayleigh, onda de Love)

Una onda de Rayleigh es una onda superficial cuyo movimiento de partícula se encuentra en el plano vertical que contiene la dirección de propagación: movimientos P y S polarizados verticalmente acoplados por la superficie libre, con un movimiento de partícula elíptico, retrógrado en la superficie de un semiespacio uniforme, cuya amplitud decae con la profundidad a lo largo de aproximadamente una longitud de onda. Una onda de Love es una onda superficial cuyo movimiento de partícula es horizontal y perpendicular a la dirección de propagación: ondas S polarizadas horizontalmente atrapadas en capas más lentas que la roca situada debajo.

En un semiespacio uniforme con $\lambda = \mu$ (un sólido de Poisson) la onda de Rayleigh viaja a $0.9194\,\beta$ y no se dispersa (ejercicio 6). Una onda de Love no puede existir en un semiespacio uniforme en absoluto: necesita la estratificación, y su velocidad está comprendida entre las velocidades S de la capa y de la roca subyacente.

Teorema 1.8 (Dispersión de las ondas de Love en una capa sobre un semiespacio)

Una capa de espesor $H$, módulo de cizalla $\mu_1$ y velocidad S $\beta_1$ descansa sobre un semiespacio de módulo de cizalla $\mu_2$ y velocidad S $\beta_2 > \beta_1$. Una onda de Love de frecuencia angular $\omega$ y número de onda $k$ se propaga con la velocidad de fase $c = \omega/k$, $\beta_1 < c < \beta_2$, si y solo si

$$\tan\!\left(kH q_1\right) = \frac{\mu_2 q_2}{\mu_1 q_1}, \qquad q_1 = \sqrt{\frac{c^2}{\beta_1^2} - 1},\quad q_2 = \sqrt{1 - \frac{c^2}{\beta_2^2}} .$$

El modo fundamental es la raíz con $0 < kHq_1 < \pi/2$; los sobretonos tienen $kHq_1$ en $(n\pi, n\pi + \pi/2)$.

Para un periodo dado la ecuación se resuelve para $c$ por bisección; la velocidad de grupo se obtiene de $U = \mathrm d\omega/\mathrm dk$. A periodos cortos ($kH$ grande) la raíz está cerca de $c = \beta_1$: la onda vive en la capa. A periodos largos ($kH \to 0$) debe cumplirse $kHq_1 \to 0$ con el lado derecho también pequeño, lo que obliga a $c \to \beta_2$: la onda llega hasta el interior del semiespacio y viaja a su velocidad. Entre ambos extremos, la velocidad de fase aumenta de forma constante con el periodo.

Definición 1.9 (Dispersión de las ondas superficiales)

La dispersión de las ondas superficiales es la dependencia de las velocidades de fase y de grupo de una onda de Rayleigh o de Love respecto de su periodo. En una Tierra estratificada surge porque una onda de periodo más largo alcanza mayor profundidad (aproximadamente un tercio de una longitud de onda para el modo fundamental) y promedia roca más rápida; las curvas de dispersión $c(T)$ y $U(T)$ medidas entre un terremoto y una estación, o entre dos estaciones, se invierten para obtener la velocidad S en función de la profundidad.

En este modelo la velocidad de fase a $20\,\mathrm{s}$ es $3.79\,\mathrm{km}\,\mathrm{s}^{-1}$ y a $100\,\mathrm{s}$ $4.444\,\mathrm{km}\,\mathrm{s}^{-1}$; las velocidades de grupo son 3.384 y $4.333\,\mathrm{km}\,\mathrm{s}^{-1}$. En un sismograma registrado a unos pocos miles de kilómetros de distancia, los periodos largos llegan primero y las ondas de 20 s, al final, de modo que el tren de ondas se extiende en una oscilación larga que se acorta lentamente. Medir la velocidad de grupo a muchos periodos y en muchas trayectorias da mapas del espesor de la corteza y del manto superior, independientes de los tiempos de recorrido de las ondas de cuerpo.

1.3 La Tierra resuena: modos normales

Una onda superficial de periodo muy largo tiene una longitud de onda comparable al radio terrestre; da la vuelta al globo y se encuentra consigo misma. Las ondas que interfieren constructivamente son las ondas estacionarias de toda la esfera, sus modos normales, y tras un gran terremoto son la última señal que se extingue. Su teoría, un problema de valores propios para la Tierra elástica y autogravitante, queda fuera de este capítulo; su descripción, no.

Los modos de una Tierra con simetría esférica forman dos familias. En los modos toroidales, ${}_nT_\ell$, el movimiento es horizontal, una torsión de capas esféricas unas contra otras, sin cambio de radio ni de volumen: son las ondas estacionarias de las ondas de Love. En los modos esferoidales, ${}_nS_\ell$, el movimiento tiene componentes radial y horizontal, y la Tierra cambia de forma: son las ondas estacionarias de las ondas de Rayleigh, acopladas a las ondas de cuerpo y a la gravedad. El índice $\ell$ es el grado angular, el número de líneas nodales del patrón en la superficie; $n$ cuenta los nodos del movimiento con el radio, con $n = 0$ para el modo fundamental. El modo esferoidal más grave, ${}_0S_2$, estira la Tierra alternativamente en un esferoide prolato y en uno oblato; el modo toroidal más grave, ${}_0T_2$, retuerce los hemisferios norte y sur en sentidos opuestos. Sus periodos observados son:

Cada periodo es un promedio sobre todo el planeta de sus módulos elásticos y de su densidad, ponderado según las zonas donde se desplaza el modo: los modos toroidales solo perciben el manto sólido y la corteza, y los esferoidales también el núcleo. Ajustar los periodos de centenares de modos es la única forma sismológica de pesar el interior: restringe la densidad, no solo las velocidades, y PREM se construyó para ajustarlos en conjunto con los tiempos de recorrido. La rotación y la elipticidad de la Tierra desdoblan cada modo de grado $\ell$ en $2\ell + 1$ frecuencias ligeramente distintas, el análogo del efecto Zeeman.

Historia — La Tierra después del terremoto chileno

El terremoto chileno de 1960 fue registrado por sismómetros pendulares de periodo largo y por deformímetros, instrumentos que miden el cambio de longitud de una barra o un hilo horizontal largo anclado en la roca. El análisis espectral de unos pocos días de registro mostró picos nítidos en los periodos previstos de los modos más graves: toda la Tierra había entrado en resonancia, y se cumplieron las predicciones de la teoría elástica de una Tierra esférica y autogravitante. El registro del deformímetro de Isabella, en California, fue analizado de nuevo en 1978 con métodos más modernos: mostró con claridad el desdoblamiento de ${}_0T_3$, ${}_0T_4$, ${}_0S_4$ y ${}_0S_5$ por la rotación y la elipticidad de la Tierra, y confirmó el de ${}_0S_2$ y ${}_0S_3$.

La sismología de modos normales, nacida de un solo terremoto, utiliza hoy las redes de estaciones de banda muy ancha que registran cada gran terremoto.

1.4 Problemas lineales inversos

Todos los métodos de este capítulo hasta ahora calculan datos a partir de un modelo. La tarea de la obtención de imágenes es la inversa.

Definición 1.10 (Problema directo, problema inverso)

Sea un modelo de la Tierra descrito por un vector $\mathbf m$ de $M$ parámetros y las observaciones por un vector $\mathbf d$ de $N$ datos. El problema directo consiste en calcular los datos $\mathbf d = g(\mathbf m)$ que produciría un modelo dado, mediante una teoría física. El problema inverso consiste en estimar el modelo a partir de datos que llevan errores, $\mathbf d = g(\mathbf m) + \mathbf e$, y en indicar con qué precisión queda determinada cada parte del modelo. El problema es lineal cuando $g(\mathbf m) = \mathsf G\mathbf m$ para una matriz fija $N\times M$ $\mathsf G$, el núcleo.

Ajustar una recta a tiempos de recorrido es un problema lineal inverso con $M = 2$ (ejercicio 7); una tomografía de tiempos de recorrido es una con miles de parámetros. Se presentan tres casos. Si $N > M$ y los datos restringen cada parámetro, el problema está sobredeterminado: ningún modelo ajusta cada dato exactamente, y la mejor solución se busca en sentido estadístico. Si alguna combinación de parámetros deja los datos sin cambios ($\mathsf G\mathbf m_0 = \mathbf 0$ para algún $\mathbf m_0 \ne \mathbf 0$), el problema está infradeterminado en esa dirección, y ninguna cantidad de datos de ese tipo lo corregirá. Los problemas reales son mixtos: sobredeterminados en unas partes del modelo e infradeterminados en otras.

Teorema 1.11 (Estimación por mínimos cuadrados ponderados)

Sea $\mathbf d = \mathsf G\mathbf m + \mathbf e$, donde los errores $\mathbf e$ tienen media cero y la matriz de covarianza $\mathsf C_d$ (diagonal $\sigma_i^2$ para errores independientes). Si la matriz $\mathsf N = \mathsf G^{\mathsf T}\mathsf C_d^{-1}\mathsf G$ es invertible, el modelo que minimiza el desajuste

$$\chi^2(\mathbf m) = (\mathbf d - \mathsf G\mathbf m)^{\mathsf T}\mathsf C_d^{-1} (\mathbf d - \mathsf G\mathbf m)$$

es la estimación por mínimos cuadrados ponderados

$$\hat{\mathbf m} = \mathsf N^{-1}\mathsf G^{\mathsf T}\mathsf C_d^{-1}\mathbf d ,$$

es insesgada, y su matriz de covarianza es $\mathsf C_{\hat m} = \mathsf N^{-1}$.

Los pesos hacen lo que pide la intuición: un dato con la mitad de error cuenta cuatro veces más. La matriz de covarianza da a cada parámetro su error típico $\sqrt{(\mathsf N^{-1})_{jj}}$ y las correlaciones entre ellos. Si el modelo es correcto y los errores son gaussianos, el $\chi^2$ mínimo sigue una ley de $\chi^2$ con $N - M$ grados de libertad: un valor muy superior a $N - M$ indica que los errores se subestimaron o que el modelo es demasiado simple; muy inferior, que los errores se sobreestimaron.

Cuando $\mathsf N$ es singular o casi singular, la estimación no existe o está dominada por el ruido amplificado en las direcciones mal restringidas. El remedio habitual es buscar un modelo que ajuste razonablemente los datos y que no sea innecesariamente grande.

Definición 1.12 (Mínimos cuadrados amortiguados, curva de compromiso)

Para un parámetro de amortiguamiento $\varepsilon > 0$, la estimación por mínimos cuadrados amortiguados es el modelo que minimiza

$$\Phi_\varepsilon(\mathbf m) = \chi^2(\mathbf m) + \varepsilon^2\|\mathbf m\|^2 , \qquad \mathbf m_\varepsilon = \left(\mathsf N + \varepsilon^2\mathsf I\right)^{-1} \mathsf G^{\mathsf T}\mathsf C_d^{-1}\mathbf d .$$

La curva de compromiso es la curva que se traza al variar $\varepsilon$, el par (tamaño del modelo $\|\mathbf m_\varepsilon\|$, desajuste $\chi^2(\mathbf m_\varepsilon)$): un amortiguamiento pequeño ajusta los datos de cerca con un modelo grande y rugoso; un amortiguamiento grande da un modelo pequeño y un ajuste pobre.

La fórmula se obtiene como en el Teorema 1.11: el gradiente de $\Phi_\varepsilon$ es $-2\mathsf G^{\mathsf T}\mathsf C_d^{-1}(\mathbf d - \mathsf G\mathbf m) + 2\varepsilon^2\mathbf m$, y $\mathsf N + \varepsilon^2 \mathsf I$ es siempre invertible, ya que sus autovalores son los de $\mathsf N$ (todos $\ge 0$) aumentados en $\varepsilon^2$. El amortiguamiento suele elegirse cerca del codo de la curva de compromiso, o de modo que $\chi^2$ sea aproximadamente el número de datos, el desajuste que cabe esperar solo por los errores. El precio es que la estimación deja de ser insesgada.

Definición 1.13 (Matriz de resolución del modelo)

Para una estimación lineal $\mathbf m_\varepsilon = \mathsf H_\varepsilon \mathbf d$ del modelo $\mathbf m$ a partir de los datos $\mathbf d = \mathsf G\mathbf m + \mathbf e$, la matriz de resolución del modelo es $\mathsf R = \mathsf H_\varepsilon\mathsf G$, de modo que

$$\mathbf m_\varepsilon = \mathsf R\,\mathbf m + \mathsf H_\varepsilon\mathbf e .$$

Para mínimos cuadrados amortiguados, $\mathsf R = (\mathsf N + \varepsilon^2\mathsf I)^{-1} \mathsf N$. La fila $j$ de $\mathsf R$ muestra cómo la estimación del parámetro $j$ promedia los parámetros verdaderos; $\mathsf R = \mathsf I$ significa resolución perfecta.

Si $\mathsf N$ tiene los autovectores ortonormales $\mathbf u_k$ y los autovalores $\lambda_k$, entonces $\mathsf R = \sum_k \frac{\lambda_k}{\lambda_k + \varepsilon^2} \mathbf u_k\mathbf u_k^{\mathsf T}$: las direcciones que los datos restringen bien ($\lambda_k \gg \varepsilon^2$) pasan casi sin cambios, mientras que las que restringen mal quedan suprimidas. La traza de $\mathsf R$ cuenta los parámetros que los datos resuelven efectivamente, y la covarianza de la estimación, $\mathsf H_\varepsilon\mathsf C_d\mathsf H_\varepsilon^{\mathsf T}$, se reduce a medida que disminuye la resolución: el amortiguamiento intercambia resolución por estabilidad.

1.5 Tomografía sísmica

Definición 1.14 (Tomografía sísmica, residuo de tiempo de recorrido)

El residuo de tiempo de recorrido de una llegada es el tiempo de recorrido observado menos el tiempo predicho por un modelo de referencia, $\delta t = T_{\text{obs}} - T_{\text{ref}}$, tras corregir la localización del foco y el tiempo origen, la elevación de la estación y la elipticidad de la Tierra. La tomografía sísmica es la estimación de las desviaciones tridimensionales de las velocidades sísmicas respecto de un modelo de referencia, $\delta v(\mathbf x)$, a partir de un conjunto grande de residuos (o de formas de onda), resolviendo un problema inverso.

Se escribe $s = 1/v$ para la lentitud. El tiempo de recorrido a lo largo de un rayo $\Gamma$ es $T = \int_\Gamma s\,\mathrm dl$; pero el propio rayo depende de las velocidades, lo que parece volver el problema no lineal. La respuesta de primer orden es que eso no importa.

Proposición 1.15 (Residuo de tiempo de recorrido linealizado)

Si la lentitud del modelo de referencia $s_0$ se perturba hasta $s_0 + \delta s$, el tiempo de recorrido entre dos puntos fijos cambia, a primer orden en $\delta s$, en

$$\delta t = \int_{\Gamma_0}\delta s\,\mathrm dl ,$$

siendo la integral tomada a lo largo del rayo $\Gamma_0$ del modelo de referencia; la anomalía relativa de velocidad es $\delta v/v_0 = -\delta s/s_0$ al mismo orden.

Se divide la región en $M$ celdas de perturbación uniforme de lentitud $\delta s_j$. El residuo del rayo $i$ pasa a ser $\delta t_i = \sum_j L_{ij}\,\delta s_j$, donde $L_{ij}$ es la longitud del rayo de referencia $i$ en la celda $j$: un problema inverso lineal cuyo núcleo es la matriz de longitudes de rayo, casi todas cuyas entradas son nulas. Los residuos son pequeños: en un manto a $8\,\mathrm{km}\,\mathrm{s}^{-1}$, una región de $200\,\mathrm{km}$ de anchura con una velocidad $3\,\%$ mayor adelanta un rayo vertical unos $0.7\,\mathrm{s}$ (ejercicio 9).

Método 1.16 (Tomografía de tiempos de recorrido)

1. Se toman los tiempos de llegada de una fase (por ejemplo, P) en todas las estaciones, para terremotos cuya localización es bien conocida; se calculan los residuos frente a una tabla de tiempos de recorrido, corregidos por elipticidad y por la elevación de la estación.
2. Se trazan los rayos de referencia y se divide la región en celdas (o se desarrolla $\delta s$ en funciones suaves); se rellena la matriz de longitudes de rayo.
3. Se asigna a cada residuo su error, a partir de la precisión de la lectura, y se resuelve el problema de mínimos cuadrados amortiguados para un intervalo de valores de amortiguamiento; se dibuja la curva de compromiso y se elige el amortiguamiento cerca de su codo.
4. Se comprueba la resolución: se calculan los residuos que un modelo sintético conocido produciría a través de los mismos rayos, se añade ruido y se invierten con el mismo amortiguamiento, y se comparan los resultados.
5. Se interpretan únicamente las anomalías que la prueba muestra que los rayos pueden recuperar, entendiendo sus amplitudes como cotas inferiores.

Definición 1.17 (Prueba del tablero de ajedrez)

Una prueba del tablero de ajedrez es una prueba de resolución en la que el modelo sintético es un patrón regular de celdas alternativamente rápidas y lentas: los residuos sintéticos calculados a través de la geometría real de los rayos, con ruido de magnitud similar a los errores reales, se invierten con la misma parametrización y amortiguamiento. Donde se recupera el patrón, pueden resolverse estructuras de ese tamaño; donde aparece difuminado o se pierde, no pueden.

Las figuras siguientes aplican la prueba a una sección de modelo de $800\,\mathrm{km}$ de anchura y $400\,\mathrm{km}$ de profundidad, dividida en 128 celdas de $50\,\mathrm{km}$. Veinte focos se sitúan a lo largo de su base y treinta y dos estaciones a lo largo de su parte superior, y solo se usan los rayos dentro de $45{}^{\circ}$ de la vertical: 480 rayos que todos ascienden casi verticalmente, la geometría de los rayos de terremotos lejanos que ascienden bajo una red. El modelo sintético tiene bloques de $2\times2$ celdas, $3\,\%$ más rápidos o más lentos alrededor de $8\,\mathrm{km}\,\mathrm{s}^{-1}$; los residuos llevan ruido gaussiano de $0.05\,\mathrm{s}$.

La prueba ilustra los tres límites de toda imagen tomográfica. La resolución viene fijada por los rayos: donde todos los rayos son paralelos, una anomalía se extiende a lo largo de ellos, y una celda atravesada por pocos rayos apenas queda restringida. El amortiguamiento reduce las amplitudes: las anomalías recuperadas son cotas inferiores. Además, la imagen depende del modelo de referencia y de la parametrización: solo se representan las desviaciones respecto de la referencia, y nada menor que una celda. Con estas cautelas, las imágenes globales de ondas P y S muestran las losas del volumen de primer año (en su capítulo sobre subducción) como regiones rápidas tabulares que descienden desde las fosas, algunas aplanándose sobre la discontinuidad de 670 km, otras hundiéndose en el manto inferior; y amplias regiones lentas bajo el Pacífico y África sobre el núcleo, con conductos lentos más estrechos bajo algunos puntos calientes, que son los candidatos a plumas del manto del mismo volumen.

Sobre el terreno — Una estación de banda ancha y un despliegue de fondo oceánico

Una estación permanente de banda ancha es un lugar tranquilo: el sismómetro se apoya sobre un pilar de hormigón anclado al lecho rocoso, en una cámara enterrada o en un túnel corto, bajo un escudo térmico, porque una corriente de aire o un cambio de unas pocas centésimas de grado inclina el instrumento más que un terremoto lejano. Sus tres componentes, vertical, norte y este, se digitalizan de forma continua y se envían a un centro de datos con sello de tiempo por satélite. Para obtener imágenes de una región bajo el océano, un equipo de campaña en cambio lanza sismómetros de fondo oceánico desde un barco: cada uno se hunde sujeto a un ancla, registra durante un año en el fondo marino y vuelve a la superficie cuando el barco regresa y envía una señal acústica para soltarlo. La tripulación lee entonces su deriva del reloj, corrige los registros y añade miles de llegadas a la tomografía.

1.6 Ejercicios

Ejercicio 1.1 ★

Un rayo P sale de la superficie ($R = 6371\,\mathrm{km}$, $v_0 = 5.80\,\mathrm{km}\,\mathrm{s}^{-1}$) con el parámetro del rayo $p = 500\,\mathrm{s}\,\mathrm{rad}^{-1}$. (a) Calcúlese su ángulo de salida. (b) PREM da velocidades P de $11.066\,\mathrm{km}\,\mathrm{s}^{-1}$ a $771\,\mathrm{km}$ de profundidad y $11.245\,\mathrm{km}\,\mathrm{s}^{-1}$ a $871\,\mathrm{km}$. Calcúlese $r/v$ a ambas profundidades y hállese la profundidad del punto de retorno del rayo por interpolación lineal.

Solución

Solución de Ejercicio 1.1.

(a) $\sin i_0 = pv_0/R = 500\times5.80/6371 = 0.4552$, así que $i_0 = 27.1{}^{\circ}$. (b) A $771\,\mathrm{km}$, $r = 5600\,\mathrm{km}$ y $r/v = 5600/11.066 = 506.1\,\mathrm{s}\,\mathrm{rad}^{-1}$; a $871\,\mathrm{km}$, $r = 5500\,\mathrm{km}$ y $r/v = 5500/11.245 = 489.1\,\mathrm{s}\,\mathrm{rad}^{-1}$. El rayo alcanza su punto de retorno donde $r/v = p = 500$: $z_t = 771 + 100\times(506.1 - 500)/(506.1 - 489.1) = 771 + 36 = 807\,\mathrm{km}$.

Ejercicio 1.2 ★

Los tiempos P calculados con PREM en este capítulo son $369.6\,\mathrm{s}$ a $30{}^{\circ}$ y $455.5\,\mathrm{s}$ a $40{}^{\circ}$. (a) Estímese la pendiente de la curva de tiempos de recorrido entre ellos, en $\mathrm{s}$ por grado y en $\mathrm{s}\,\mathrm{rad}^{-1}$. (b) Dedúzcase la velocidad aparente del frente de onda a lo largo de la superficie, en $\mathrm{km}\,\mathrm{s}^{-1}$. (c) ¿Cuál es la velocidad P en el punto de retorno del rayo que llega a $35{}^{\circ}$, si su punto de retorno está en el radio $5450\,\mathrm{km}$?

Solución

Solución de Ejercicio 1.2.

(a) $(455.5 - 369.6)/10 = 8.59\,\mathrm{s}$ por grado, y $8.59\times180/\pi = 492.2\,\mathrm{s}\,\mathrm{rad}^{-1}$. (b) Un grado de arco es $6371\pi/180 = 111.19\,\mathrm{km}$, así que el frente de onda barre la superficie a $111.19/8.59 = 12.94\,\mathrm{km}\,\mathrm{s}^{-1}$, mucho más rápido que cualquier roca en la superficie: el rayo llega casi vertical. (c) En el punto de retorno $v = r_t/p = 5450/492.2 = 11.07\,\mathrm{km}\,\mathrm{s}^{-1}$, una velocidad de onda P del manto inferior.

Ejercicio 1.3 ★

En el modelo de onda de Love de este capítulo la capa tiene $\beta_1 = 3.5\,\mathrm{km}\,\mathrm{s}^{-1}$, $\rho_1 = 2700\,\mathrm{kg}\,\mathrm{m}^{-3}$ y el semiespacio $\beta_2 = 4.5\,\mathrm{km}\,\mathrm{s}^{-1}$, $\rho_2 = 3300\,\mathrm{kg}\,\mathrm{m}^{-3}$. (a) Calcúlense los dos módulos de cizalla. (b) A $20\,\mathrm{s}$ la velocidad de fase es $3.79\,\mathrm{km}\,\mathrm{s}^{-1}$: calcúlese la longitud de onda y compárese con el espesor de la capa, $35\,\mathrm{km}$. (c) Con velocidades de grupo de $3.384\,\mathrm{km}\,\mathrm{s}^{-1}$ a $20\,\mathrm{s}$ y $4.333\,\mathrm{km}\,\mathrm{s}^{-1}$ a $100\,\mathrm{s}$, ¿en cuántos minutos llegan las ondas de 20 s después de las de 100 s en una estación situada a $4000\,\mathrm{km}$?

Solución

Solución de Ejercicio 1.3.

(a) $\mu_1 = \rho_1\beta_1^2 = 2700\times3500^2 = 3.31 \times 10^{10}\,\mathrm{Pa}$ ($33.1\,\mathrm{GPa}$); $\mu_2 = 3300\times4500^2 = 6.68 \times 10^{10}\,\mathrm{Pa}$ ($66.8\,\mathrm{GPa}$). (b) $\lambda = cT = 3.79\times20 = 75.8\,\mathrm{km}$, 2.17 veces el espesor de la capa: la onda penetra bien en el semiespacio, lo que explica por qué su velocidad ya supera $\beta_1$. (c) $4000/3.384 = 1182\,\mathrm{s}$ y $4000/4.333 = 923\,\mathrm{s}$: las ondas de 20 s llegan $259\,\mathrm{s}$, o $4.3\,\mathrm{min}$, después de las ondas de 100 s.

Ejercicio 1.4 ★

Los periodos observados de ${}_0S_2$ y ${}_0T_2$ son $53.89\,\mathrm{min}$ y $44.01\,\mathrm{min}$. (a) Calcúlense sus frecuencias en $\mathrm{mHz}$. (b) ¿Cuántas oscilaciones de cada uno completa la Tierra en un día? (c) Dos picos espectrales pueden separarse solo si el registro dura al menos el inverso de su diferencia de frecuencias. ¿Cuánto debe durar un registro para separar estos dos modos?

Solución

Solución de Ejercicio 1.4.

(a) ${}_0S_2$: $1/(53.89\times60) = 3.093 \times 10^{-4}\,\mathrm{Hz} = 0.3093\,\mathrm{mHz}$; ${}_0T_2$: $1/(44.01\times60) = 0.3787\,\mathrm{mHz}$. (b) $1440/53.89 = 26.7$ y $1440/44.01 = 32.7$ oscilaciones al día. (c) La diferencia es $0.0694\,\mathrm{mHz}$, cuyo inverso es $14\,400\,\mathrm{s}$: al menos $4.0\,\mathrm{h}$ de registro, en la práctica mucho más, ya que los picos se ensanchan por la atenuación y se desdoblan.

Ejercicio 1.5 ★★

En una esfera de radio $R$ y velocidad uniforme $v$ los rayos son cuerdas rectas. (a) Demuéstrese que el rayo que llega a $\Delta$ tiene el tiempo de recorrido $T = 2R\sin(\Delta/2)/v$. (b) Con $R = 6371\,\mathrm{km}$ y $v = 10.0\,\mathrm{km}\,\mathrm{s}^{-1}$, calcúlese $T$ a $30{}^{\circ}$ y $90{}^{\circ}$, y compárese con los tiempos de PREM, $369.6\,\mathrm{s}$ y $779.7\,\mathrm{s}$. (c) Explíquese por qué la esfera uniforme es demasiado rápida a $30{}^{\circ}$ y demasiado lenta a $90{}^{\circ}$.

Solución

Solución de Ejercicio 1.5.

(a) La cuerda entre dos puntos de la esfera separados por el ángulo $\Delta$ en el centro tiene la longitud $2R\sin(\Delta/2)$, recorrida a $v$. (b) A $30{}^{\circ}$: $2\times6371\times\sin15{}^{\circ}/10.0 = 329.8\,\mathrm{s}$, frente a $369.6\,\mathrm{s}$; a $90{}^{\circ}$: $2\times6371\times\sin45{}^{\circ}/10.0 = 901.0\,\mathrm{s}$, frente a $779.7\,\mathrm{s}$. (c) El rayo hasta $30{}^{\circ}$ retorna en el manto superior, donde la velocidad real es de unos 8 a $10\,\mathrm{km}\,\mathrm{s}^{-1}$, por debajo de la de $10\,\mathrm{km}\,\mathrm{s}^{-1}$ uniforme; el rayo hasta $90{}^{\circ}$ pasa la mayor parte de su trayecto en el manto inferior, a 11 a $13.7\,\mathrm{km}\,\mathrm{s}^{-1}$, y se curva hacia dentro, de modo que la Tierra real es más rápida.

Ejercicio 1.6 ★★

En un sólido de Poisson ($\lambda = \mu$, luego $\alpha^2 = 3\beta^2$) la onda de Rayleigh de velocidad $c$ satisface, con $\xi = c^2/\beta^2$, $(2 - \xi)^2 = 4\sqrt{1 - \xi/3}\,\sqrt{1 - \xi}$. (a) Elévese al cuadrado y demuéstrese que las raíces satisfacen $(\xi - 4)(3\xi^2 - 12\xi + 8) = 0$. (b) Demuéstrese que la única raíz con $0 < \xi < 1$ es $\xi = 2 - 2/\sqrt3$, y calcúlese $c/\beta$. (c) En la corteza superior de ak135, $\beta = 3.46\,\mathrm{km}\,\mathrm{s}^{-1}$: calcúlese la velocidad de Rayleigh y el retraso de la onda de Rayleigh respecto de la onda S a $1000\,\mathrm{km}$.

Solución

Solución de Ejercicio 1.6.

(a) Elevando al cuadrado, $(2 - \xi)^4 = 16(1 - \xi/3)(1 - \xi)$. Desarrollando: $16 - 32\xi + 24\xi^2 - 8\xi^3 + \xi^4 = 16 - \tfrac{64}{3}\xi + \tfrac{16}{3}\xi^2$, así que $\xi^4 - 8\xi^3 + \tfrac{56}{3}\xi^2 - \tfrac{32}{3}\xi = 0$, es decir, $\xi\,(3\xi^3 - 24\xi^2 + 56\xi - 32) = 0$, y $3\xi^3 - 24\xi^2 + 56\xi - 32 = (\xi - 4)(3\xi^2 - 12\xi + 8)$. La raíz $\xi = 0$ no es una onda. (b) Las raíces de la cuadrática son $2 \pm 2/\sqrt3$, es decir 3.155 y 0.8453; solo $\xi = 0.8453$ está en $(0, 1)$ (y $\xi = 4$ no), y satisface la ecuación sin elevar al cuadrado, siendo ambos miembros $4/3$. Así que $c/\beta = \sqrt{0.8453} = 0.9194$. (c) $c = 0.9194\times3.46 = 3.181\,\mathrm{km}\,\mathrm{s}^{-1}$; sobre $1000\,\mathrm{km}$, $314.4 - 289.0 = 25.3\,\mathrm{s}$ después de la onda S.

Ejercicio 1.7 ★★

Las llegadas de onda de cabeza en cuatro estaciones de una línea de refracción, a 50, 100, 150 y $200\,\mathrm{km}$ del disparo, son 9.0, 15.3, 21.5 y $27.9\,\mathrm{s}$, con errores de 0.1, 0.1, 0.2 y $0.2\,\mathrm{s}$. Ajústese $T = a + bx$ por mínimos cuadrados ponderados: escríbanse $\mathsf G$ y $\mathsf C_d$, fórmense las ecuaciones normales y dense $a$, $b$, sus errores estándar, la velocidad $1/b$ con su error y el mínimo $\chi^2$. ¿Es aceptable el ajuste?

Solución

Solución de Ejercicio 1.7.

$\mathsf G$ tiene las filas $(1, x_i)$ y $\mathsf C_d^{-1} = \operatorname{diag} (100, 100, 25, 25)$ ($\mathrm{s}^{-2}$). Las ecuaciones normales son

$$\begin{pmatrix} 250 & 23750 \\ 23750 & 2\,812\,500 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} 3665 \\ 418125 \end{pmatrix},$$

con determinante $1.390\,625 \times 10^{8}$. Por tanto, $b = (250\times418125 - 23750\times3665)/1.390\,625 \times 10^{8} = 0.125\,75\,\mathrm{s}\,\mathrm{km}^{-1}$ y $a = (2812500\times3665 - 23750\times418125)/1.390\,625 \times 10^{8} = 2.713\,\mathrm{s}$. La covarianza es $\mathsf N^{-1}$: $\sigma_b = \sqrt{250/1.390\,625 \times 10^{8}} = 0.001\,34\,\mathrm{s}\,\mathrm{km}^{-1}$, $\sigma_a = \sqrt{2812500/1.390\,625 \times 10^{8}} = 0.142\,\mathrm{s}$. La velocidad es $1/b = 7.95\,\mathrm{km}\,\mathrm{s}^{-1}$, con el error $\sigma_b/b^2 = 0.08\,\mathrm{km}\,\mathrm{s}^{-1}$: una velocidad del manto más superficial. Los residuos son $-0.001$, $0.011$, $-0.076$ y $0.036\,\mathrm{s}$, y $\chi^2 = 0.19$ para $4 - 2 = 2$ grados de libertad: menor de lo esperado, así que los errores indicados son, si acaso, generosos; el ajuste es aceptable.

Ejercicio 1.8 ★★

Un rayo atraviesa dos celdas, de longitudes $30\,\mathrm{km}$ y $40\,\mathrm{km}$, y tiene un residuo de $0.50\,\mathrm{s}$. (a) ¿Por qué es singular $\mathsf N = \mathsf G^{\mathsf T} \mathsf G$? (b) Entre los modelos que ajustan el residuo exactamente, el de menor norma es $\mathbf m = \mathsf G^{\mathsf T}(\mathsf G\mathsf G^{\mathsf T})^{-1}\mathbf d$: calcúlese. (c) Calcúlese la matriz $\mathsf R = \mathsf G^{\mathsf T}(\mathsf G\mathsf G^{\mathsf T})^{-1}\mathsf G$ y su traza. (d) Si la anomalía verdadera es $0.01\,\mathrm{s}\,\mathrm{km}^{-1}$ solo en la primera celda, ¿qué se recupera?

Solución

Solución de Ejercicio 1.8.

(a) $\mathsf G = (30\;\;40)$ es una única fila: $\mathsf G^{\mathsf T}\mathsf G$ es de $2\times2$ y de rango uno, por tanto singular; cualquier modelo con $30m_1 + 40m_2 = 0$ puede añadirse sin cambiar el dato. (b) $\mathsf G\mathsf G^{\mathsf T} = 900 + 1600 = 2500$, así que $\mathbf m = (30, 40)\times0.50/2500 = (0.006, 0.008)\,\mathrm{s}\,\mathrm{km}^{-1}$. (c) $\mathsf R = \frac{1}{2500} \begin{pmatrix} 900 & 1200 \\ 1200 & 1600 \end{pmatrix} = \begin{pmatrix} 0.36 & 0.48 \\ 0.48 & 0.64 \end{pmatrix}$, de traza 1: un dato resuelve una combinación de parámetros. (d) $\mathsf R\,(0.01, 0) = (0.0036, 0.0048)\,\mathrm{s}\,\mathrm{km}^{-1}$: la anomalía se reparte a lo largo del rayo, más en la celda donde el rayo es más largo.

Ejercicio 1.9 ★★

Una región de losa del manto es un $3\,\%$ más rápida que su entorno, donde $v = 9.0\,\mathrm{km}\,\mathrm{s}^{-1}$. (a) Calcúlese el residuo de un rayo que atraviesa $200\,\mathrm{km}$ de ella en vertical, mediante la Proposición 1.15, y de forma exacta. (b) Lo mismo para un rayo que atraviesa un espesor de $100\,\mathrm{km}$ de ella con un ángulo de $40{}^{\circ}$ respecto a la normal a la losa. (c) ¿Cuál es el error relativo de la linealización?

Solución

Solución de Ejercicio 1.9.

(a) Linealizando: $\delta t = -(\delta v/v)L/v = -0.03\times200/9.0 = -0.667\,\mathrm{s}$; exactamente, $L/(1.03v) - L/v = 200/9.27 - 200/9.0 = -0.647\,\mathrm{s}$. (b) El trayecto en la losa es $100/\cos40{}^{\circ} = 130.5\,\mathrm{km}$: $-0.03\times130.5/9.0 = -0.435\,\mathrm{s}$ linealizado y $-0.422\,\mathrm{s}$ exacto. (c) El valor exacto es el linealizado dividido por $1.03$: un error de $3\,\%$, del orden de $\delta v/v$ en sí, como cabe esperar de una teoría de primer orden.

Ejercicio 1.10 ★★★

En el Teorema 1.8, un sobretono de orden $n$ existe solo para periodos suficientemente cortos. (a) Demuéstrese que en su frecuencia de corte la velocidad de fase es $\beta_2$ y $kHq_1 = n\pi$. (b) Dedúzcase la frecuencia angular de corte $\omega_n = n\pi\beta_1\beta_2/\bigl(H\sqrt{\beta_2^2 - \beta_1^2}\bigr)$. (c) Calcúlese el periodo de corte del primer sobretono para el modelo del capítulo ($H = 35\,\mathrm{km}$, $\beta_1 = 3.5$, $\beta_2 = 4.5\,\mathrm{km}\,\mathrm{s}^{-1}$). (d) Demuéstrese que el modo fundamental no tiene frecuencia de corte y que su velocidad de fase tiende a $\beta_1$ para periodos cortos.

Solución

Solución de Ejercicio 1.10.

(a) El sobretono existe mientras su parte en el semiespacio decae, $c < \beta_2$; en su frecuencia de corte $c = \beta_2$, así que $q_2 = 0$ y $\tan(kHq_1) = 0$, con $kHq_1$ en $[n\pi, n\pi + \pi/2)$: $kHq_1 = n\pi$. (b) Con $k = \omega/\beta_2$ y $q_1 = \sqrt{\beta_2^2/\beta_1^2 - 1} = \sqrt{\beta_2^2 - \beta_1^2}/\beta_1$, $\omega_nH\sqrt{\beta_2^2 - \beta_1^2}/(\beta_1\beta_2) = n\pi$, que es la fórmula. (c) $\sqrt{4.5^2 - 3.5^2} = \sqrt8 = 2.828$; $\omega_1 = \pi\times3.5\times4.5/(35\times2.828) = 0.4998\,\mathrm{rad}\,\mathrm{s}^{-1}$, un periodo de $2\pi/0.4998 = 12.57\,\mathrm{s}$: a periodos más largos solo existe el modo fundamental. existe. (d) Para $n = 0$ la frecuencia de corte es cero: el modo fundamental existe a todos los periodos. A periodos cortos $kH \to \infty$ mientras $kHq_1$ se mantiene por debajo de $\pi/2$, así que $q_1 \to 0$ y $c \to \beta_1$.

Ejercicio 1.11 ★★★

Tómese el rayo del ejercicio 8, con error de datos unitario. (a) Demuéstrese que, para cualquier $\varepsilon > 0$, $(\mathsf G^{\mathsf T}\mathsf G + \varepsilon^2\mathsf I)^{-1}\mathsf G^{\mathsf T} = \mathsf G^{\mathsf T}(\mathsf G\mathsf G^{\mathsf T} + \varepsilon^2\mathsf I)^{-1}$. (b) Dedúzcase la matriz de resolución amortiguada $\mathsf R_\varepsilon = \mathsf G^{\mathsf T}\mathsf G/(\|\mathbf g\|^2 + \varepsilon^2)$, donde $\mathbf g = (30, 40)$, y su traza. (c) Calcúlese la traza y el residuo de datos para $\varepsilon^2 = 0$, 2500 y $10\,000$ (en $\mathrm{km}^{2}$). Descríbase el compromiso.

Solución

Solución de Ejercicio 1.11.

(a) $\mathsf G^{\mathsf T}(\mathsf G\mathsf G^{\mathsf T} + \varepsilon^2\mathsf I) = (\mathsf G^{\mathsf T}\mathsf G + \varepsilon^2\mathsf I)\mathsf G^{\mathsf T}$; se multiplica por la izquierda por $(\mathsf G^{\mathsf T}\mathsf G + \varepsilon^2 \mathsf I)^{-1}$ y por la derecha por $(\mathsf G\mathsf G^{\mathsf T} + \varepsilon^2\mathsf I)^{-1}$, ambas invertibles para $\varepsilon > 0$. (b) Aquí $\mathsf G\mathsf G^{\mathsf T} = \|\mathbf g\|^2 = 2500$, un número, por tanto $\mathsf R_\varepsilon = \mathsf G^{\mathsf T}\mathsf G/(2500 + \varepsilon^2)$, de traza $2500/(2500 + \varepsilon^2)$. (c) El dato predicho es $\mathbf g\cdot\mathbf m_\varepsilon = 0.50\times2500/(2500 + \varepsilon^2)$ y el residuo $0.50\,\varepsilon^2/(2500 + \varepsilon^2)$. Para $\varepsilon^2 = 0$: traza 1, residuo 0; para 2500: traza 0.5, residuo $0.25\,\mathrm{s}$; para $10\,000$: traza 0.2, residuo $0.40\,\mathrm{s}$. Un mayor amortiguamiento da lugar a un modelo más pequeño, una resolución peor y un desajuste mayor.

Ejercicio 1.12 ★★★

En la esfera uniforme del ejercicio 5: (a) demuéstrese mediante el Teorema 1.3 que $p(\Delta) = (R/v)\cos(\Delta/2)$, y compruébese con la definición de $p$ (el ángulo de salida de una cuerda). (b) Demuéstrese que el rayo que llega a $\Delta_1$ tiene su punto de retorno en $r_1 = R\cos(\Delta_1/2)$. (c) Admitiendo que $\int_0^{\Delta_1}\operatorname{arccosh}\bigl(\cos(\Delta/2) /\cos(\Delta_1/2)\bigr)\,\mathrm d\Delta = \pi\ln\bigl(1/\cos(\Delta_1/2)\bigr)$, demuéstrese que la fórmula de Herglotz–Wiechert devuelve $r_1$ y la velocidad $v$.

Solución

Solución de Ejercicio 1.12.

(a) Por el ejercicio 5, $T = (2R/v)\sin(\Delta/2)$, así que $\mathrm dT/\mathrm d\Delta = (R/v)\cos(\Delta/2)$. Una cuerda que forma el ángulo $\Delta$ en el centro es la base de un triángulo isósceles de ángulo en el vértice $\Delta$, por lo que corta al radio en la fuente con $i_0 = (\pi - \Delta)/2$, y $p = R\sin i_0/v = (R/v)\cos(\Delta/2)$: lo mismo. (b) El punto de retorno es el punto medio de la cuerda, a la distancia $R\cos(\Delta_1/2)$ del centro; equivalentemente $r_1 = p_1v$. (c) Aquí $p(\Delta)/p_1 = \cos(\Delta/2)/\cos(\Delta_1/2)$, y la fórmula da $r_1 = R\exp\bigl(-\ln(1/\cos(\Delta_1/2))\bigr) = R\cos(\Delta_1/2)$, como en (b); la velocidad es $r_1/p_1 = R\cos(\Delta_1/2)\big/ \bigl((R/v)\cos(\Delta_1/2)\bigr) = v$.

1.7 Problema: la losa bajo la estación K

Problema 1.1

Problema de fin de semana — rayos, ondas de Love y una pequeña tomografía bajo una estación ficticia

La estación K se sitúa sobre un continente, encima de una zona de subducción. Este problema emplea sus registros de tres maneras: rayos en su corteza, la dispersión de las ondas de Love a través de ella, y una pequeña tomografía del manto subyacente.

Parte I — Rayos en un gradiente. A distancias lo bastante cortas para despreciar la curvatura terrestre, la velocidad de las ondas P bajo K aumenta linealmente con la profundidad, $v(z) = v_0 + gz$, con $v_0 = 5.0\,\mathrm{km}\,\mathrm{s}^{-1}$ y $g = 0.05\,\mathrm{s}^{-1}$ ($z$ en km).

1. En capas planas, la ley de Snell mantiene constante $p = \sin i/v$. Un rayo sale K con $\sin i_0 = 0.8$: calcúlense $i_0$ y $p$.
2. ¿A qué profundidad tiene su punto de retorno, y cuál es allí la velocidad?
3. Demuéstrese que $\mathrm di/\mathrm ds = pg$ a lo largo del rayo, de modo que es un arco de circunferencia; dése su radio y la profundidad de su centro.
4. Calcúlese la distancia horizontal $X$ a la que el rayo regresa a la superficie.
5. Demuéstrese que $\mathrm dT = \mathrm di/(g\sin i)$ y calcúlese el tiempo de recorrido $T$ de forma exacta.
6. Calcúlense $X$ y $T$ para $p = 0.15\,\mathrm{s}\,\mathrm{km}^{-1}$, y compruébese que $\Delta T/\Delta X$ entre los dos rayos está comprendido entre los dos valores de $p$.

Parte II — Ondas de Love y la corteza. La corteza bajo K es una capa de velocidad S $\beta_1 = 3.5\,\mathrm{km}\,\mathrm{s}^{-1}$ y densidad $\rho_1 = 2.7\,\mathrm{g}\,\mathrm{cm}^{-3}$ sobre un manto con $\beta_2 = 4.5\,\mathrm{km}\,\mathrm{s}^{-1}$, $\rho_2 = 3.3\,\mathrm{g}\,\mathrm{cm}^{-3}$; su espesor $H$ es desconocido. Entre K y una estación vecina la velocidad de fase de las ondas de Love se mide: $3.99\,\mathrm{km}\,\mathrm{s}^{-1}$ a $25\,\mathrm{s}$ y $4.25\,\mathrm{km}\,\mathrm{s}^{-1}$ a $40\,\mathrm{s}$.

7. Calcúlense $\mu_1$ y $\mu_2$ en $\mathrm{GPa}$.
8. A $25\,\mathrm{s}$, calcúlense $k$, $q_1$ y $q_2$.
9. Resuélvase la ecuación de dispersión del modo fundamental para $H$.
10. Repítase lo mismo a $40\,\mathrm{s}$. ¿Son consistentes los dos resultados?
11. Estímese la velocidad de grupo entre los dos periodos como $\Delta\omega/\Delta k$, y compárese con las velocidades de fase.
12. Calcúlese la longitud de onda a cada periodo y explíquese por qué los dos periodos sondean la corteza y el manto superior en distinta proporción.

Parte III — Una tomografía en cuatro celdas. Bajo K, un bloque de manto de $200\,\mathrm{km}$ por $200\,\mathrm{km}$ en una sección vertical se divide en cuatro celdas cuadradas de $100\,\mathrm{km}$: 1 (superior izquierda), 2 (superior derecha), 3 (inferior izquierda) y 4 (inferior derecha); la fosa y el océano quedan a la izquierda. La velocidad de referencia es $v_0 = 8.0\,\mathrm{km}\,\mathrm{s}^{-1}$. En un experimento idealizado, seis rayos rectos la atraviesan: A a lo largo de la fila superior (celdas 1 y 2), B a lo largo de la fila inferior (3 y 4), C por la columna izquierda (1 y 3), D por la columna derecha (2 y 4), E a lo largo de la diagonal por 1 y 4, F a lo largo de la diagonal por 2 y 3. Sus residuos son:

rayo	A	B	C	D	E	F
$\delta t$ ($\mathrm{s}$)	$0.16$	$-0.48$	$-0.36$	$0.11$	$-0.15$	$-0.27$
error $\sigma$ ($\mathrm{s}$)	0.05	0.05	0.05	0.05	0.10	0.10

El modelo $\mathbf m$ se toma como las perturbaciones de lentitud de las celdas en $\mathrm{ms}\,\mathrm{km}^{-1}$ ($1\,\mathrm{ms}\,\mathrm{km}^{-1} = 10^{-3}\,\mathrm{s}\,\mathrm{km}^{-1}$), de modo que las entradas de $\mathsf G$ son las longitudes de los rayos en $\mathrm{km}$ divididas entre 1000.

13. Escríbase $\mathsf G$, una matriz $6\times4$.
14. Demuéstrese que $\mathsf N = \mathsf G^{\mathsf T}\mathsf C_d^{-1}\mathsf G$ tiene las filas $(10, 4, 4, 2)$, $(4, 10, 2, 4)$, $(4, 2, 10, 4)$, $(2, 4, 4, 10)$.
15. Calcúlese $\mathbf b = \mathsf G^{\mathsf T}\mathsf C_d^{-1}\mathbf d$.
16. Compruébese que $(1,1,1,1)$, $(1,-1,-1,1)$, $(1,1,-1,-1)$ y $(1,-1,1,-1)$ son autovectores de $\mathsf N$, con autovalores 20, 4, 8 y 8, y dedúzcase que $\mathsf N^{-1}$ tiene las filas $(11, -4, -4, 1)/80$, $(-4, 11, 1, -4)/80$, $(-4, 1, 11, -4)/80$ y $(1, -4, -4, 11)/80$.
17. Calcúlese la estimación por mínimos cuadrados ponderados $\hat{\mathbf m}$.
18. Conviértase en anomalías relativas de velocidad $\delta v/v_0$ en por ciento.
19. Dése el error estándar de cada celda en $\mathrm{ms}\,\mathrm{km}^{-1}$ y en por ciento, y la correlación entre las celdas 1 y 2, y entre las celdas 1 y 4. Explíquense sus signos.
20. Calcúlense los residuos de los seis rayos tras el ajuste y el mínimo $\chi^2$. ¿Es compatible con los errores indicados?

Parte IV — Amortiguamiento, resolución y la losa.

21. Para $\varepsilon^2 = 4$ (en las unidades de $\mathsf N$), dense los autovalores de $\mathsf R = (\mathsf N + 4\mathsf I)^{-1}\mathsf N$ y demuéstrese que $\mathsf R$ tiene la diagonal $2/3$, las entradas $1/12$ entre celdas que comparten un rayo de fila o de columna, y $0$ entre las celdas 1 y 4, y entre las celdas 2 y 3.
22. Calcúlese la estimación amortiguada de la celda 3 y su anomalía de velocidad. ¿Por qué es menor que la estimación sin amortiguar?
23. El modelo amortiguado tiene $\chi^2 = 15.8$ y norma $2.91\,\mathrm{ms}\,\mathrm{km}^{-1}$; el no amortiguado, $\chi^2 = 0.50$ y $4.24\,\mathrm{ms}\,\mathrm{km}^{-1}$. ¿Cuál se conservaría, y por qué este problema no necesita amortiguamiento?
24. Léase la fila 3 de $\mathsf R$ amortiguada: ¿con qué celdas se mezcla la estimación de la celda 3?
25. ¿Cuál celda es rápida y cuál lenta? Interprétese la celda rápida, debajo y del lado del océano de K, en términos de la zona de subducción.
26. Dése la anomalía de velocidad P de la celda de la losa en por ciento, con su error estándar.

Solución

Solución de Problema 1.1.

1. $i_0 = \arcsin 0.8 = 53.13{}^{\circ}$ y $p = 0.8/5.0 = 0.16\,\mathrm{s}\,\mathrm{km}^{-1}$.

2. El rayo alcanza su punto de retorno donde $v = 1/p = 6.25\,\mathrm{km}\,\mathrm{s}^{-1}$, es decir, en $z = (6.25 - 5.0)/0.05 = 25\,\mathrm{km}$.

3. Derivando $\sin i = pv$ a lo largo del rayo, $\cos i\,\mathrm di/ \mathrm ds = pg\,\mathrm dz/\mathrm ds = pg\cos i$, así que $\mathrm di/\mathrm ds = pg$: la dirección gira a ritmo constante, con la curvatura de un círculo de radio $1/(pg) = 1/(0.16\times0.05) = 125\,\mathrm{km}$. En el punto de retorno, el rayo es horizontal a $25\,\mathrm{km}$ de profundidad, y el centro se encuentra verticalmente a $125\,\mathrm{km}$ por encima de él: en $z = -100\,\mathrm{km}$, es decir, a $v_0/g$ por encima de la superficie, donde la velocidad extrapolada se anularía.

4. En la fuente, el radio del círculo es perpendicular al rayo, por lo que forma el ángulo $i_0$ con la horizontal, y el centro queda a $125\cos i_0$ de distancia horizontal; por simetría $X = 2\times125\times\cos i_0 = 2\times125\times0.6 = 150\,\mathrm{km}$, o $X = 2\sqrt{1 - p^2v_0^2}/(pg)$.

5. $\mathrm dT = \mathrm ds/v$ con $\mathrm ds = \mathrm di/(pg)$ y $v = \sin i/p$: $\mathrm dT = \mathrm di/(g\sin i)$. Entonces

$$T = \frac2g\int_{i_0}^{\pi/2}\frac{\mathrm di}{\sin i} = \frac2g\ln\cot\frac{i_0}2 = \frac2g\operatorname{arccosh}\frac1{pv_0} .$$

Aquí $\cot(26.57{}^{\circ}) = 2$, así que $T = 40\ln2 = 27.73\,\mathrm{s}$.

6. Para $p = 0.15$: $pv_0 = 0.75$, $\sqrt{1 - 0.5625} = 0.6614$, $X = 2\times0.6614/0.0075 = 176.4\,\mathrm{km}$ y $T = 40\operatorname{arccosh} (1.3333) = 40\times0.7954 = 31.81\,\mathrm{s}$. Entonces $\Delta T/\Delta X = 4.08/26.4 = 0.155\,\mathrm{s}\,\mathrm{km}^{-1}$, entre 0.15 y 0.16, como $\mathrm dT/\mathrm dX = p$ lo exige.

7. $\mu_1 = 2700\times3500^2 = 33.1\,\mathrm{GPa}$ y $\mu_2 = 3300\times4500^2 = 66.8\,\mathrm{GPa}$.

8. $k = 2\pi/(cT) = 2\pi/(3.99\times25) = 0.062\,99\,\mathrm{rad}\,\mathrm{km}^{-1}$; $q_1 = \sqrt{(3.99/3.5)^2 - 1} = 0.5474$; $q_2 = \sqrt{1 - (3.99/4.5)^2} = 0.4624$.

9. $\tan(kHq_1) = 66.825\times0.4624/(33.075\times0.5474) = 1.707$; el modo fundamental toma $kHq_1 = \arctan1.707 = 1.041$, así que $H = 1.041/(0.06299\times0.5474) = 30.2\,\mathrm{km}$.

10. A $40\,\mathrm{s}$: $k = 2\pi/170 = 0.03696$, $q_1 = 0.6888$, $q_2 = 0.3287$, $\tan(kHq_1) = 0.964$, $kHq_1 = 0.7671$ y $H = 0.7671/(0.03696\times0.6888) = 30.1\,\mathrm{km}$. Los dos resultados coinciden dentro de la precisión de las velocidades de fase: la corteza tiene un espesor de unos $30\,\mathrm{km}$.

11. $\omega = 0.2513$ y $0.1571\,\mathrm{rad}\,\mathrm{s}^{-1}$, $k = 0.06299$ y $0.036\,96\,\mathrm{rad}\,\mathrm{km}^{-1}$: $U \approx 0.0942/0.02603 = 3.62\,\mathrm{km}\,\mathrm{s}^{-1}$, menor que ambas velocidades de fase, como debe ser cuando $c$ crece con el periodo.

12. $\lambda = 3.99\times25 = 99.8\,\mathrm{km}$ y $4.25\times40 = 170\,\mathrm{km}$. El modo fundamental muestrea hasta una fracción de su longitud de onda: la onda de 25 s se sitúa sobre todo en la corteza, la de 40 s llega bien dentro del manto, así que su velocidad está más cerca de $\beta_2$. Dos periodos que fijan el mismo $H$ es una comprobación de que el modelo de dos capas es adecuado.

13. Con las columnas de las celdas 1 a 4 y las filas A a F, en unidades de $1000\,\mathrm{km}$ ($\sqrt2\times100/1000 = 0.1414$):

$$\mathsf G = \begin{pmatrix} 0.1 & 0.1 & 0 & 0 \\ 0 & 0 & 0.1 & 0.1 \\ 0.1 & 0 & 0.1 & 0 \\ 0 & 0.1 & 0 & 0.1 \\ 0.1414 & 0 & 0 & 0.1414 \\ 0 & 0.1414 & 0.1414 & 0 \end{pmatrix}.$$

14. $\mathsf C_d^{-1} = \operatorname{diag}(400, 400, 400, 400, 100, 100)$. Un rayo de fila o columna aporta $0.1^2\times400 = 4$, un rayo diagonal $0.1414^2\times100 = 2$. La celda 1 es atravesada por A y C (4 cada uno) y por E (2): 10. Las celdas 1 y 2 comparten el rayo A: 4; las celdas 1 y 3 comparten C: 4; las celdas 1 y 4 comparten E: 2; y así sucesivamente por simetría, lo que da las cuatro filas indicadas.

15. Un rayo de fila o columna aporta $0.1\times400\,\delta t = 40\, \delta t$ a cada una de sus celdas, y un rayo diagonal $0.1414\times100\,\delta t = 14.14\,\delta t$. Celda 1: $40(0.16 - 0.36) + 14.14\times(-0.15) = -10.12$; celda 2: $40(0.16 + 0.11) + 14.14\times(-0.27) = 6.98$; celda 3: $40(-0.48 - 0.36) + 14.14\times(-0.27) = -37.42$; celda 4: $40(-0.48 + 0.11) + 14.14\times(-0.15) = -16.92$.

16. Para $(1,1,1,1)$, cada fila de $\mathsf N$ suma 20. Para $(1,-1,-1,1)$, la fila 1 da $10 - 4 - 4 + 2 = 4$ y las demás siguen el mismo patrón. Para $(1,1,-1,-1)$, $10 + 4 - 4 - 2 = 8$, y para $(1,-1,1,-1)$, $10 - 4 + 4 - 2 = 8$. Los vectores son ortogonales con norma 2, así que $\mathsf N^{-1} = \sum_k \mathbf u_k\mathbf u_k^{\mathsf T}/(4\lambda_k)$; su elemento $(1,1)$ es $\tfrac14(\tfrac1{20} + \tfrac14 + \tfrac18 + \tfrac18) = \tfrac{11}{80}$, su elemento $(1,2)$ es $\tfrac14(\tfrac1{20} - \tfrac14 + \tfrac18 - \tfrac18) = -\tfrac4{80}$, su elemento $(1,4)$ es $\tfrac14(\tfrac1{20} + \tfrac14 - \tfrac18 - \tfrac18) = \tfrac1{80}$, y así sucesivamente.

17. $\hat{\mathbf m} = \mathsf N^{-1}\mathbf b$: $\hat m_1 = (11\times(-10.121) - 4\times6.982 - 4\times(-37.418) - 16.921)/80 = -0.081$; del mismo modo $\hat m_2 = 1.844$, $\hat m_3 = -3.706$ y $\hat m_4 = -0.931\,\mathrm{ms}\,\mathrm{km}^{-1}$.

18. $\delta v/v_0 = -v_0\,\delta s = -8.0\times10^{-3}\,\hat m$, en tanto por ciento $-0.8\,\hat m$: celda 1 $+0.07$, celda 2 $-1.48$, celda 3 $+2.96$, celda 4 $+0.75\,\%$.

19. Todo elemento diagonal de $\mathsf N^{-1}$ vale $11/80$: cada celda tiene el error estándar $\sqrt{0.1375} = 0.371\,\mathrm{ms}\,\mathrm{km}^{-1}$, es decir, $0.30\,\%$. La correlación de las celdas 1 y 2 es $-4/11 = -0.36$: comparten el rayo A, y una lentitud asignada a una se puede compensar con la otra. Las celdas 1 y 4, que solo comparten la diagonal E, de menor precisión, están débilmente y positivamente correlacionadas, $+1/11 = +0.09$.

20. Los residuos predichos son $0.176$, $-0.464$, $-0.379$, $0.091$, $-0.143$ y $-0.263\,\mathrm{s}$; los residuos tras el ajuste son $-0.016$, $-0.016$, $0.019$, $0.019$, $-0.007$ y $-0.007\,\mathrm{s}$, y $\chi^2 = 0.50$ para $6 - 4 = 2$ grados de libertad: coherente con los errores indicados, si acaso algo pequeños.

21. Los autovectores son los de $\mathsf N$ y los autovalores son $\lambda/(\lambda + 4)$: $20/24 = 5/6$, $4/8 = 1/2$, $8/12 = 2/3$ dos veces. Como en la pregunta 16, $R_{11} = \tfrac14(\tfrac56 + \tfrac12 + \tfrac23 + \tfrac23) = \tfrac23$; $R_{12} = \tfrac14(\tfrac56 - \tfrac12 + \tfrac23 - \tfrac23) = \tfrac1{12}$; $R_{13} = \tfrac14(\tfrac56 - \tfrac12 - \tfrac23 + \tfrac23) = \tfrac1{12}$; $R_{14} = \tfrac14(\tfrac56 + \tfrac12 - \tfrac23 - \tfrac23) = 0$.

22. Proyéctese $\mathbf b$ sobre los autovectores: $\mathbf b\cdot\mathbf u_k/4 = -14.37$, $0.849$, $12.80$ y $-9.40$ para los autovalores 20, 4, 8 y 8. La estimación amortiguada los divide por $\lambda + 4$ en lugar de por $\lambda$; para la celda 3 (componentes $1, -1, -1, 1$), $m_3 = -14.37/24 - 0.849/8 - 12.80/12 - 9.40/12 = -2.555\,\mathrm{ms}\,\mathrm{km}^{-1}$, una anomalía de velocidad de $+2.04\,\%$. El amortiguamiento lleva cada componente hacia cero, sobre todo las peor restringidas.

23. El modelo sin amortiguar, que ajusta los datos dentro de sus errores ($\chi^2 = 0.50$ para 2 grados de libertad), mientras que el amortiguado los desajusta notablemente ($15.8$). El amortiguamiento solo es necesario cuando $\mathsf N$ es singular o casi lo es; aquí su autovalor menor, 4, no es pequeño, la estimación es estable y sus errores son moderados.

24. La fila 3 es $(1/12, 0, 2/3, 1/12)$: la estimación amortiguada de la celda 3 es dos tercios de su valor verdadero más una doceava parte de las celdas 1 y 4, sus vecinas a lo largo del rayo C y del rayo B, y nada de la celda 2.

25. La celda 3, la celda inferior de la izquierda, es rápida en unos $3\,\%$; la celda 2 es lenta. Una losa fría y densa de litosfera oceánica que se hunde bajo el continente es más rápida que el manto circundante, así que la celda rápida, profunda y en el lado oceánico de K, es la losa; la celda lenta situada encima, al otro lado, se encuentra donde la cuña del manto está caliente.

26. A partir de las preguntas 18 y 19, la anomalía de velocidad P de la celda de la losa es $(3.0 \pm 0.3)\,\%$.