A predictive evolution model of rocky coasts

Resumen

El retroceso de las costas acantiladas es un fenómeno muy extendido sobre los litorales rocosos expuestos a la incidencia combinada de los procesos marinos y meteorológicos que se dan en la franja costera. Este fenómeno se revela violenta y ocasionalmente como movimientos gravitacionales del terreno que pueden causar pérdidas materiales y/o humanas. Aunque el conocimiento de estos riesgos de erosión resulta de vital importancia para la correcta gestión de la costa, su predicción es complicada. Los modelos de predicción publicados son escasos y con importantes inconvenientes: extrapolación (extienden información de registros históricos); empíricos (sobre registros históricos estudian la respuesta al cambio de un parámetro); estocásticos (determinan la cadencia y la magnitud de los eventos futuros extrapolando las distribuciones de probabilidad extraídas de catálogos históricos); proceso-respuesta, de estabilidad y propagación del error inexplorada; en EDP, computacionalmente costosos y poco exactos. En este trabajo se desarrolla un modelo combinado del tipo proceso-respuesta basado en incorporar un balance de fuerzas de los mecanismos que actúan sobre el proceso erosivo en el frente del acantilado. El modelo simula la evolución espaciotemporal de un perfil-2D del acantilado, formado por materiales heterogéneos, acoplando la dinámica marina con la evolución del terreno en cada periodo de marea. Integra en cada periodo una función de erosión, dependiente de la pendiente de la zona afectada, que se desplaza sobre la onda de marea. Se ha estudiado el error de discretización del modelo y su propagación en el tiempo a partir de las soluciones exactas para los 2 primeros periodos de marea para diferentes aproximaciones numéricas en la integración y de la pendiente. Los resultados obtenidos han permitido justificar las discretizaciones que minimizan el error y los métodos de aproximación más adecuados.

Abstract

The recession of coastal cliffs is a widespread phenomenon on the rocky shores that are exposed to the combined incidence of marine and meteorological processes that occur in the shoreline. This phenomenon is revealed violently and occasionally as gravitational movements of the ground and can cause material or human losses. Their prediction is difficult; however it is basic for the proper coastal management and the clear understanding of erosion risks. There are several models of the coastal cliff recession processes. From the stochastic type models based on historical erosion events, to the theoretical models as described by Eikonal or Boussinesq equations. In this work an intermediate solution is adopted (process-response model), such as models based on simplified balance of forces as trigger mechanisms involved in the erosion of the rocks. This model fits the marine dynamics: sea level changes, tidal range and wave transformation; along with the evolution of the land: the slope of the surf zone, erosion and rock falls. So far these models have been limited to small slope cliffs, since its numerical stability and the propagation of errors are unknowns. For these reasons, the development of a generic model that reproduces the spatial and temporal evolution of a cliff-2D profile (platform, beach and slope) consisting on semi-consolidated heterogeneous materials, is presented. In addition, the computational implementation, the study of different numerical resolution techniques and the produced errors, are also exposed and analysed.

Palabras clave

Acantilado ; Erosión ; Peligro geológico ; Recesión costera ; Simulación

Keywords

Cliff ; Erosion ; Geological hazard ; Coastal recession ; Computer model

1. Introducción

Uno de los mayores retos en las ciencias ambientales es desarrollar modelos de sistemas dinámicos naturales que permitan simular su estado futuro cuando se combinan los procesos humanos con los propios. Sin embargo, para la mayoría de los sistemas naturales este reto no tiene solución trivial. Los grandes avances logrados durante las 2 últimas décadas en el estudio de sistemas dinámicos no lineales ha reforzado la necesidad de desarrollar modelos que permitan la simulación o la reproducción realista de la naturaleza no lineal del comportamiento de ciertos sistemas medioambientales.

Las distintas aproximaciones al conocimiento de la franja costera como sistema dinámico reconocen un comportamiento complejo que es el resultado de las interacciones entre un gran número de procesos naturales y humanos que actúan sobre un amplio rango de escalas espaciales y temporales [1] . Estas relaciones se producen a través de mecanismos de retroalimentación y dan lugar a comportamientos dinámicos fuertemente no lineales. En este sistema se ha descrito el acoplamiento entre fenómenos erosivos fuertemente dependientes de los procesos marinos, como características de oleaje y marea o transporte de sedimentos, así como de las condiciones del medio terrestre, como litología, grado de fracturación o meteorización [2] . A estos fenómenos hay que añadir la concentración de usos y actividades propia de este ámbito territorial, lo que introduce los efectos de las intervenciones antrópicas directas o indirectas en su evolución. Como consecuencia, las interrelaciones que se producen en este espacio de transición mar-tierra son complejas y llegan a producir situaciones de riesgo, bien por la propia dinámica del sistema natural o bien derivadas de la ocupación y alteración de los procesos naturales.

La naturaleza no lineal de la morfodinámica costera es bien conocida para playas arenosas y estuarios [3] . Sin embargo, el 80% del perímetro costero terrestre está formado por acantilados rocosos [4] , y este tipo de costas constituyen el 60% de la franja costera peninsular [5] . Estas zonas presentan una destacada incidencia de movimientos del terreno de rango extraordinario [6] que hacen especialmente vulnerables algunas áreas donde hay poblaciones, complejos turísticos, industrias, etc. Recuérdense, por ejemplo, los sucesos más recientes tristemente acaecidos en la playa de los Gigantes (1/11/2009, isla de Tenerife, 2 fallecidos) o en la playa María Luisa (21/08/2009, Algarve, Portugal, 5 fallecidos). A pesar de la importancia socioeconómica del medio litoral acantilado, los modelos para reproducir o pronosticar fielmente su comportamiento son muy escasos.

En este artículo se trata el diseño y el estudio empírico del error de un modelo matemático de erosión de acantilados costeros. El desarrollo del artículo se centra en el análisis empírico de los errores cometidos en la resolución del modelo, con el fin de mejorar y depurar las inestabilidades numéricas presentes en modelos anteriormente publicados y que hasta la fecha no han sido realizados por otros autores. Esto permitirá el uso del modelo con mayor confianza en los resultados obtenidos, disminuyendo las restricciones de uso y mejorando, además, la eficiencia del cálculo. Este modelo está basado en una formulación simplificada de los procesos que combinan sus efectos en la franja costera y permite simular la evolución temporal del perfil vertical de la costa acantilada. La figura 1 a muestra los principales elementos de este tipo de costas afectadas por la recesión costera. El alcance temporal que permite la simulación con este modelo es entre uno y 100 años, lo que se conoce como mesoescala.

Figura 1. Sección esquemática de la franja costera modelada. a) Bloque diagrama que representa los elementos geomorfológicos y los procesos erosivos que dan lugar a la recesión costera. b) Sistema de referencia, discretización espacial del perfil, y funciones más relevantes incorporadas en el modelo.

En primer lugar se clasifican los modelos de recesión de costas rocosas más relevantes en la actualidad. Posteriormente, y a partir de los modelos teóricos y empíricos descritos en la literatura, se fundamenta y explica el modelo aquí desarrollado. Se presentan varias aproximaciones numéricas para la resolución del modelo, y sobre ellas se evalúa el error cometido, a partir de una serie de soluciones exactas, para diferentes momentos de las simulaciones, lo que permite estudiar cómo se propaga el error a lo largo del tiempo. No llega a tratarse de un estudio de la estabilidad y aproximación teórica del modelo, lo cual queda fuera del propósito de este artículo. Finalmente se exponen las conclusiones referentes a los resultados obtenidos y se comentan algunos de los módulos que se están incorporando al modelo para incrementar su capacidad y su potencial.

2. Modelos previos de recesión de costas rocosas

La morfodinámica de la recesión de un acantilado, aunque se produzca continuamente una erosión gradual en el frente o al pie del mismo, puede ser un proceso intermitente, pues su colapso ocurre repentina e infrecuentemente [7] . Así, las variaciones en el medio geológico, en el medio marino y en la climatología condicionan de manera determinante la compleja dinámica de la recesión costera. Por tanto, para desarrollar un modelo predictivo es fundamental tener en cuenta la mayor cantidad de información posible, las relaciones de retroalimentación existentes y el acoplamiento a distintas escalas que existe entre los procesos que controlan dicha morfodinámica.

Los primeros modelos publicados, como los de Kamphuis J. W. [8] y Sunamura T. [9] , tomaron la tasa de erosión como un promedio que cuantifica la distancia de terreno perdida por la recesión costera a lo largo del tiempo (generalmente un año). Su valor es resultado del balance entre la resultante de las fuerzas marinas y atmosféricas incidentes y la capacidad resistente del material. Aunque se constataba que en su evolución contribuían numerosos procesos acoplados a diferentes escalas, el posible efecto causado quedaba englobado en una u otra fuerza del balance. De esta forma, a pesar del elevado número de procesos que interactúan en la recesión de un acantilado, gran parte de los modelos consideran simplificadamente que la erosión, o cambio de posición por unidad de tiempo, en un punto x del mismo depende de: a) la fuerza de impacto (hidrodinámica e hidrostática [9] ) del oleaje (F_w); b) el tiempo (t) en que esta fuerza actúa sobre el punto y que depende del ciclo de marea, y c) las propiedades resistivas (F_R) del macizo rocoso. En los primeros modelos se propone que el desplazamiento de un punto de la franja costera a causa de la erosión marina a lo largo del tiempo t es:

$\delta (x,t)=f(F_{W}(x),F_{R}(x),t)\leftrightarrow {\frac {dx}{dt}}\approx {\frac {F_{W}(x,t)}{F_{R}(x,t)}}$

( 1)

Se han realizado varios intentos para asociar F_w y F_R con parámetros físicos medibles. Entre otros, los experimentos de Kamphuis [8] , [10]  and [11] proporcionaron una expresión de f(·) en la que la resistencia del material se representa mediante un término a calibrar con experimentos. Posteriormente, el modelo de Sunamura [9] , deducido igualmente a partir de experimentos sobre una pared de hormigón en un tanque de oleaje artificial, propone que F_R sea directamente proporcional a la resistencia a compresión (σ_c) de la roca expuesta, y formula además la tasa temporal de erosión instantánea δ(x,t) proporcional al logaritmo de la razón entre ambas fuerzas, F_R y F_W , aplicadas en la posición x . La identificación del parámetro más apropiado para representar la resistencia de la roca a la erosión F_R ha sido investigada en numerosos estudios sobre el deterioro de acantilados. La resistencia a compresión σ_c es el más frecuentemente adoptado [12] . También se han propuesto, aunque con menos frecuencia, la resistencia a tracción [13] , la cohesión [14] , la resistencia al corte [10] o el módulo de Young [15] , entre otros. Entre estos parámetros, la resistencia a la compresión está directamente ligada con la mayoría de ellos y se trata de un valor ampliamente utilizado, por lo que en los modelos que incorporan F_R se suele tomar σ_c[12] como característica geomecánica del macizo intacto.

A partir de estas aproximaciones experimentales han surgido varias líneas de desarrollo de modelos (tabla 1 ) para la estimación/predicción de tasas de recesión en acantilados. La primera se basada en los catálogos históricos sobre los cuales se practicaban una serie de regresiones (lineales) y extrapolaciones (simples) de las tendencias pasadas para cuantificar las tasas en el futuro [16] , [17]  and [18] . Sin embargo, estos modelos no son capaces de proporcionar una distribución precisa de las recesiones futuras [18] , habida cuenta de su limitada extensión temporal y, por ende, su incapacidad para reflejar la variabilidad e incertidumbre de todos los procesos involucrados. Sin embargo, algunos trabajos recogen estos inconvenientes y tratan de crear un marco estructural para reflejar sobre los catálogos históricos los cambios en los factores clave que controlan la recesión costera. Por ejemplo, Lee [19] aplica un procedimiento de aproximación que, en un marco probabilístico, está basado en un árbol de decisiones o sucesos. En dicho árbol propone una función de distribución de las recesiones, partiendo de su distribución en los catálogos históricos existentes para cada lugar.

Tabla 1. Resumen de los métodos de predicción de recesión de acantilados costeros

Tipo de modelo	Subtipo	Referencias
Extrapolación	Regresión lineal	[16] , [17]  and [18]
	Árboles de sucesos	[19]
Empíricos	Proyección de datos históricos	[20]
	Regla de Bruun	[15] , [21]  and [22]
Probabilísticos	Distribución doble	[23]
	Distribuciones múltiples	[24]  and [41]
	Estocástico	[25]
Proceso-respuesta	Transporte de sedimentos	[28]  and [29]
	Modelo combinado	[9] , [13] , [26] , [27] , [40]  and [48]

En una segunda línea se pueden englobar los modelos empíricos. Estos son modelos muy simples que tratan de aportar una indicación de los posibles cambios en la tendencia de recesión futura como respuesta al cambio en algún factor determinante (p. ej., aumento en el nivel del mar, altura del acantilado…). Dentro de esta clase de modelos existen 2 grandes tipos: los que son simplemente una proyección de datos históricos [20] y los basados en la regla de Bruun [9] , [15] , [21]  and [22] . Estos modelos no tienen en cuenta la incertidumbre y la variabilidad que envuelven a los procesos que generan la recesión costera. Los datos de tasas de recesión históricos suelen ser bastante limitados, y por tanto se revelan como una imagen parcial de los eventos acaecidos. Sin embargo, la incierta relación entre la recesión pasada y la futura se puede acomodar adoptando un modelo probabilístico y tomando conciencia del comportamiento actual. Bajo esta condición aparece otro tipo de modelos al que pertenecen los métodos probabilísticos y los modelos de proceso-respuesta.

La tercera línea está representada por los modelos puramente probabilísticos, como por ejemplo el denominado cliffplan[23] . En estos modelos se incorpora una versión simplista de los procesos implicados, para así poder realizar múltiples simulaciones con un coste computacional relativamente bajo. El resultado de nuevo es una función de distribución de la distancia de recesión en un instante de la simulación. La amplia variabilidad introducida en los múltiples muestreos aleatorios de las funciones de probabilidad de las variables del modelo exige que sus parámetros sean calibrados utilizando inventarios históricos. Dentro de este grupo se pueden destacar varios subtipos de modelos (tabla 1 ). Los modelos de distribución doble utilizan una función de distribución para representar el intervalo entre eventos y otra para la magnitud del evento [23] . Existen otros modelos, como el de Furlan C. [24] , que utilizan distribuciones múltiples donde todos los parámetros influyentes en la recesión costera se extraen de distribuciones de probabilidad (p. ej., resistencia de la roca, condiciones de oleaje…). Por último, a este grupo pertenecen los modelos estocásticos, que se diferencian de los demás en que las probabilidades de un suceso cambian con el tiempo, es decir, son modelos dinámicos [25] .

La cuarta línea desarrolla los modelos de proceso-respuesta como el desarrollado por Sunamura [9] , el desarrollado por Trenhaile [13] o como el denominado scape[26]  and [27] . Este tipo de modelos se basan en el conocimiento y la descripción matemática de las interacciones entre el medio marino, el litoral y el propio acantilado, aunque están todavía en su fase inicial de desarrollo. Para incorporar la incertidumbre y la variabilidad, los modelos se suelen enmarcar en un procedimiento de análisis probabilista [19] que permite abordar la complejidad del sistema utilizando modelos simples de cada proceso. En este tipo de modelos, como en los anteriores, los registros históricos se utilizan para su calibración. Otro problema que presenta este tipo de modelos es que la mayor parte de las interacciones físicas se estudian en laboratorio, lo que limita su conocimiento y exactitud. A este grupo pertenecen los modelos centrados en el transporte de sedimentos [28]  and [29] y los modelos combinados [13] , [26] , [27]  and [30] .

La formulación de los modelos anteriores es relativamente simple comparada con la de las ecuaciones de Boussinesq [31] para ondas cortas irregulares y multidireccionales en aguas someras [32] . Esta opción es de un enorme coste computacional, al tener que calcular cada ola y, sobre esta, el potencial erosivo y de arrastre y su incidencia sobre la roca, para el tiempo final de simulación que puede ser del orden de siglos. De igual forma ocurre con los modelos matemáticos basados en una ecuación no lineal de tipo Eikonal [33] para simular el avance del frente de la erosión [34] . En estos se simula la propagación de una onda en un medio resistivo como un análogo al problema planteado sobre el cambio del perfil del acantilado. El frente de la onda se asemeja a la superficie de la interfase roca/aire en el pie del acantilado que va cambiando de posición a causa de la erosión. Las soluciones explícitas de este modelo, presentadas por Belov [34] , muestran cómo la intensidad de la erosión decrece exponencialmente con la altura. Este resultado es solo aplicable para acantilados desarrollados en materiales homogéneos no fracturados y en los tramos superiores del notch , si bien su forma, aun en rocas con tales características, es mucho más compleja y depende fuertemente de su estado evolutivo.

3. Desarrollo del modelo numérico

A partir del marco conceptual de funcionamiento del sistema acantilado en la franja costera, el modelo numérico desarrollado utiliza las ecuaciones simplificadas del oleaje [35] para evaluar la interacción entre la dinámica marina (mareas y oleaje), la morfología costera y la erosión de la roca en el perfil de un acantilado. El modelo no considera el fenómeno de deposición de sedimentos transportados por las corrientes y las olas, ni la erosión química, biomecánica y pluviomecánica. El modelo se basa en la hipótesis común en latitudes medias de que la acción mecánica del oleaje es el mecanismo erosivo marino más importante en acantilados abiertos sometidos a eventos tormentosos [9]  and [36] . Este proceso mecánico actúa sobre una estrecha franja, limitada por la fluctuación superior e inferior del nivel medio del mar a causa del ritmo de las mareas.

Para su formulación, se considera el sistema de referencia ZY global, con origen sobre el nivel medio del mar en el instante inicial ( fig. 1 b), manteniéndose en esta posición aunque este nivel cambie con el paso del tiempo. Se representa por y(z,t) el frente de la interfase roca/agua/aire ( fig. 1 b) cuya evolución se modela en forma del desplazamiento de este para cada punto del perfil del acantilado con cota z en cada instante de tiempo t . Durante su evolución, el perfil y(z,t) puede desplazarse hacia el medio terrestre, produciéndose la regresión de la línea de costa, o hacia el medio marino, ganándose terreno al mar (acreción). En consecuencia, la cresta del acantilado se desplaza en una cantidad R(t) que cambia a lo largo del tiempo. Es en este punto elevado donde habitualmente se practican las medidas (terrestres o aéreas) de control y seguimiento (los puntos

en la fig. 1 b) de la recesión [37] . Sobre el medio terrestre se han supuesto unas características litológicas homogéneas o con una estratificación subhorizontal (fig. 1 b). Estas se incorporan en el modelo a partir del valor de la resistencia a compresión σ_c(z) del material que se encuentre a cota z sobre el nivel medio del mar, como valor representativo de las características geomecánicas del material.

Considerando un punto de altura z situado bajo el nivel de pleamar (altura máxima afectada por la erosión en un ciclo de marea), el desplazamiento producido en el perfil por la acción erosiva durante un periodo de tiempo T se puede expresar como:

$\delta (z,T)={\frac {y(z,t+T)-y(z,t)}{T}}$

( 2)

Kamphuis [7] y Sunamura [9] suponen la expresión (2) proporcional a la razón de fuerzas incidentes (capacidad erosiva del oleaje) frente a las resistentes (capacidad resistiva del macizo): δ(z,T) ∝ F_W/ F_R . Como razón de las fuerzas en este trabajo se ha utilizado una expresión modificada de Kamphuis [7] en la que se ha explicitado el factor de proporcionalidad a la resistencia del macizo como σ_c(z) en F_R e introducido el efecto de la energía disipada por el oleaje, al aumentar el recorrido que ha de hacer la ola rompiente sobre la plataforma costera, en función de su pendiente β en F_W , resultando:

$\delta (z,T)\propto {\frac {H_{b}^{13/4}T_{b}^{3/2}}{{\sigma }_{c}(z)}}tan\beta$

( 3)

El término de la pendiente de la plataforma tan β en la expresión (3) puede relacionarse con el comportamiento del perfil y(z,t) , siendo aquella el inverso de su pendiente local: (∂_zy(z,t))^–1 .

Las fuerzas erosivas F_W generadas por el mar dependen de numerosos parámetros. De algunos se conoce o se puede cuantificar su efecto, aunque sea de manera aproximada. Estos son la altura de la ola en la rompiente H_b(t) , el periodo del oleaje T_b(t) , la variación de altura debido a las mareas w_t(t) , la amplitud de la marea A_m(t) y la variación del nivel medio del mar w_rl(t) provocado por el cambio climático o por movimientos isostáticos o tectónicos [35] . Aunque dependen también de otros muchos factores no cuantificables, como del contenido en burbujas, de los sedimentos en suspensión, de la eficiencia del impacto para provocar un efecto pistón en las fisuras y diaclasas de la roca, etc. Su incidencia, en mayor o en menor grado, de forma conjunta se traduce en una importante incertidumbre sobre la distribución de F_W en profundidad. Este problema se resuelve en el modelo agregando una función de forma p_w(z,t) que modifica proporcionalmente a F_W y que se obtiene a partir de datos experimentales [38]  and [39] . Dicha función de forma p_w(z,t) está relacionada con la longitud de onda de las olas, de su periodo T , del momento de ruptura y de su orientación θ_b(t) .

$\delta (z,T)={\frac {H_{b}^{13/4}T_{b}^{3/2}}{K{\sigma }_{c}(z)}}{\left({\frac {\partial y(z,t)}{\partial x}}\right)}^{-1}p_{w}(z,t)$

( 4)

Se llega finalmente a la expresión (4) para obtener el desplazamiento del perfil para un periodo de tiempo T . En ella se ha introducido una constante K que, obtenida mediante calibración, recoge el efecto de las posibles variaciones en la competencia del material, debido a heterogeneidades litológicas y estructurales sobre F_R . Y también se ha introducido un ajuste sobre ciertas constantes hidrodinámicas no del todo conocidas.

Como δ(z,T) es la recesión tras un tiempo T , transformándola a recesión instantánea (T →0), según Sunamura [9] , y suponiendo que el perfil y(z,t) se desplaza en dicha cantidad en los puntos a cota z afectados por la erosión:

${\begin{array}{cc}{\underset {\Delta t\rightarrow 0}{lim}}\delta (z,\Delta t)={\underset {\Delta t\rightarrow 0}{lim}}{\frac {y(z,t+\Delta t)-y(z,t)}{\Delta t}}={\frac {\partial y(z,t)}{\partial t}}&\forall {\mbox{   }}z\in [z_{min},z_{max}]\end{array}}$

( 5)

En consecuencia, supuesta y(z,t) suficientemente regular para que este límite exista y teniendo en cuenta (4), es posible formular el problema de calcular la evolución del perfil de un acantilado como:

hallar y (z , t ) ∈ C¹ ([z_min , z_max ]x [0, t_fin ]) que verifique:

${\begin{array}{c}{\frac {\partial y(z,t)}{\partial t}}={\frac {H_{b}^{13/4}T_{b}^{3/2}}{K{\sigma }_{c}(z)}}{\left({\frac {\partial y(z,t)}{\partial z}}\right)}^{-1}p_{w}(z,t);\forall z\in [z_{min},z_{max}],t\in ]0,t_{fin}]\\y(z,0)=pf(z);\forall z\in [z_{min},z_{max}];c.i.\\y(z_{min},t)=y_{min},y(z_{max},t)=y_{max},t\in [0,t_{fin}];c.c.\end{array}}$

( 6)

para unas condiciones iniciales (c.i.) del perfil del acantilado conocido pf(z) y unas condiciones de contorno (c.c.) y_min e y_max que permanecen constantes a lo largo del ciclo erosivo sobre el notch .

La formulación del problema (6) corresponde a un régimen de oleaje constante y un nivel del mar fijo. Para su generalización a regímenes variables, con trenes de olas de diferentes altitudes o periodos sobre un nivel del mar que fluctúa con la marea w_t(t) , bastará con contemplar la posible variación con el tiempo de los términos que se refieren a F_W . Para ello bastará con suponer que la erosión sobre un punto z sometido a diferentes tipos de olas es aditiva, esto es: H_b  = H_b(t) y T_b  = T_b(t); además de desplazar la función de ponderación p_w(z,t) sobre la marea de la forma p_w(w_t(t)-z) .

El modelo propuesto en (6) posee algunos aspectos comunes a los modelos de Walkden y Hall [40] , Lee [41] y Trenhaile [13] , aunque también presenta importantes diferencias. Este modelo generaliza las formulaciones anteriores de la evolución del perfil y(z,t) que se presentaban en su versión discreta. Además se ha introducido una propiedad mecánica del material reduciendo la variabilidad paramétrica sobre la que calibrar las simulaciones. El modelo aquí propuesto permite, tal y como se recoge la resistencia del material, trabajar no solo con medios geológicos heterogéneos estratificados para los que σ_c  = σ_c(z) , sino con medios con heterogeneidades estructurales si σ_c  = σ_c(z,y) . Ningún modelo de los encontrados hasta el momento en la bibliografía recoge esta posibilidad.

3.1. Comportamiento temporal asintótico del modelo

A partir de (6) es sencillo evaluar el comportamiento en el perfil de un acantilado a largo plazo en función del efecto predominante, si los efectos marinos se promedian en un intervalo temporal suficientemente largo t^∞ . Sin embargo, la falta de estacionariedad a corto y medio plazo en los procesos involucrados (mareas, oleaje, clima) afecta al comportamiento de ∂_ty(z,t) y lo hace fuertemente variable a lo largo del tiempo, lo que motiva el desarrollo de un modelo numérico de (6) si se desea estudiar su evolución a mesoescala.

Por un lado, a igualdad de condiciones de exposición a los fenómenos marinos, es posible comparar 2 acantilados con diferencias en las resistencias. En este escenario, la recesión a largo plazo puede expresarse, según (6), como:

${\frac {\partial y}{\partial t}}={\frac {1}{K}}{\frac {F_{W}(H_{b},T_{b},t^{\infty })}{F_{R}({\sigma }_{c},z)}}=$$K_{W}F\ast _{R}^{-1}$

( 7)

donde K_W  = F*_W/K , F*_W es la fuerza de la acción del mar promedio a largo plazo y F*_R es la resistencia del macizo, supuestas estacionarias e independientes del tiempo. En (7) se refleja cómo los acantilados más «duros» poseen una tasa de recesión menor.

Por otro lado, si ambos acantilados se encuentran afectados por una dinámica marina diferente, a igualdad de propiedades geomecánicas en el macizo rocoso, reorganizando la ecuación (7) puede escribirse que:

${\frac {\partial y}{\partial t}}={\frac {1}{K}}{\frac {F_{W}(H_{b},T_{b},t^{\infty })}{F_{R}({\sigma }_{c},z)}}=$$K\ast _{R}F\ast _{W}$

( 8)

lo que se interpreta como que, a largo plazo, los acantilados que sufrirán tasas de recesión mayores serán aquellos que, para características geomecánicas semejantes, estén sometidos a temporales marinos más severos. De aquí la importancia de establecer correctamente el régimen de tormentas que pueden afectar un tramo de litoral para la correcta gestión territorial.

3.2. Discretización espacial y temporal del modelo

Para la resolución numérica del problema evolutivo planteado en (6) y teniendo en cuenta (4) para un Δt suficientemente pequeño (desplazamiento instantáneo):

$y(z,t+\Delta t)=y(z,t)+{\frac {H_{b}^{13/4}(t)T_{b}^{3/2}(t)}{K{\sigma }_{c}(z)}}{\left({\frac {\partial y(z,t)}{\partial z}}\right)}^{-1}p_{w}(w_{t}(t)-$$z)\Delta t$

( 9)

Si para un Δt  ≤  T el comportamiento de H_b y T_b es constante y tampoco hay variaciones significativas en la geometría del perfil en el rango de escalas temporales consideradas para resolver la evolución de y(z,t) , entonces lo único que afecta durante un ciclo de periodo T es la variación del potencial erosivo que actúa desplazándose sobre la onda de marea ( fig. 2 a). De esta forma, durante el tiempo que actúa un ciclo de marea, el efecto en el estado y(z,t) para llegar a y(z,t+T) es resultado de la acumulación de los cambios en la posición de p_w(t) tras cada Δt :

$y(z,t+T)=y(z,t)+{\frac {H_{b}^{13/4}(t)T_{b}^{3/2}(t)}{K{\sigma }_{c}(z)}}{\left({\frac {\partial y(z,t)}{\partial z}}\right)}^{-1}\sum _{k=0}^{T/\Delta t}p_{w}(w_{t}(t+$$k\Delta t)-z)\Delta t$

( 10)

Figura 2. Modelo de onda senoidal de marea con periodo T  = 12,46 h. a) Discretización temporal Δt , sobre la que se calcula el efecto acumulado de la función pw(z,t) que afecta a un determinado z(y). b) Función pw(z,t) sobre la onda de marea en unidades (escala de grises) que representan tasa de erosión entre pendiente.

Y cuando se toma Δt  →  0, la ecuación (9) se transforma en:

$y(z,t+T)=y(z,t)+{\left(-1\right)}^{\nu }{\frac {H_{b}^{13/4}(t)T_{b}^{3/2}(t)}{K{\sigma }_{c}(z)}}{\left({\frac {\partial y(z,t)}{\partial z}}\right)}^{-1}\int _{0}^{T}p_{w}(w_{t}(t)-$$z)dt$

( 11)

ecuación semidiscreta de (6) que permite calcular la nueva posición del frente erosionado después que ha pasado un ciclo de marea de periodo T . El cambio de signo que introduce ν en (11) depende de la dirección del oleaje respecto al frente acantilado, si la erosión evoluciona de derecha a izquierda ν  = 0, y si es al contrario ν  = 1.

Con el fin de evaluar numéricamente el desplazamiento y su comportamiento frente a la elección del esquema de cálculo respecto al paso de tiempo Δt ( fig. 2 a), la integral sobre cada ciclo de marea se ha aproximado mediante diferentes fórmulas de integración numérica [42] , utilizándose Newton-Cotes (rectángulo, trapecio y Simpson) y fórmulas de cuadratura gaussiana (de 4 y 7 puntos) para cada subintervalo de tamaño Δt ( tabla 2 ) en el que se ha dividido el intervalo [0,T ]:

$\int _{0}^{T}p_{w}(w_{t}(t)-z(y))dt=\sum _{i=0}^{(T/\Delta t)-1}\int _{t+i\Delta t}^{t+(i+1)\Delta t}p_{w}(w_{t}(\tau )-$$z(y))d\tau =\sum _{i=0}^{(T/\Delta t)-1}Ip_{w}^{i}(z,t)$

( 12)

Tabla 2. Fórmulas de integración numérica de Newton-Cotes y de cuadratura gaussiana para aproximar (12)

Fórmula de integración numérica para ${\textstyle Ip_{w}^{i}(z,t)}$		Orden del error
Rectángulo (T-I)	${\textstyle \Delta tp_{w}(w_{t}(t+i\Delta t)-}$$z)$	O (Δt2 )
Trapecio (T-II)	${\textstyle {\frac {\Delta t}{2}}[p_{w}(w_{t}(t+}$$i\Delta t)-z]+p_{w}(w_{t}(t+(i+1)\Delta t)-z)]$	O (Δt3 )
Simpson (T-III)	${\textstyle {\frac {\Delta t}{6}}[p_{w}(w_{t}(t+}$$i\Delta t)-z)+4p_{w}(w_{t}(t+{\frac {(i+1)\Delta t}{2}})-$$z)+p_{w}(w_{t}(t+(i+1)\Delta t)-z)]$	O (Δt5 )
Cuadratura (n puntos)	${\textstyle {\frac {\Delta t}{2}}\sum _{j=1}^{n}c_{nj}p_{w}(w_{t}({\theta }_{j})-}$$z)$ ; θj  = (rnj  + 1)Δt  + 2(t  + iΔt )	O (Δtm )

En la fórmula de cuadratura c_nj y r_nj son los coeficientes y raíces, respectivamente, de cada tipo de integración de Gauss, siendo m el orden del error cometido en la aproximación. Si el número de puntos n es par, m  = n  + 1; pero si n es impar, m  = n  +  2. Se han seleccionado los valores de n  = 4 que corresponde con una aproximación tipo IV cuyo orden de error es O(Δt⁵) , y de n  = 7 que corresponde con una aproximación tipo V o Weddle, cuyo orden de error es O(Δt⁹) . Los coeficientes y raíces utilizados para la generación de estas fórmulas de integración gaussiana son los usuales [42] .

La discretización espacial del acantilado se ha llevado a cabo en segmentos de altura Δz ( fig. 1 b), lo que permite obtener un conjunto de N   puntos ${\textstyle {\left\{\left(y_{i},z_{i}\right)\right\}}_{i=1}^{N}}$ , tales que y_i  = y (z_i , t ), con z_{i +1}  − z_i  = Δz , cuyos extremos son z₁  = z_min y z_N  = z_max . Sobre esta división del perfil es posible calcular aproximadamente la pendiente local β evaluando numéricamente la derivada [42] de y(z_i,t) respecto a z en cada instante t de la simulación. Se han implementado diferentes técnicas de derivación numérica ( tabla 3 ) para su aproximación, lo que se utilizará para estudiar la sensibilidad de la solución frente al uso de uno u otro método.

Tabla 3. Fórmulas de derivación numérica para aproximar localmente ∂_zy(z_i,t)

Fórmula de derivación numérica para ∂zy(zi,t) .		Orden del error
2 puntos	${\textstyle {\frac {1}{\Delta z}}\left[y(z_{i+1},t)-\right.}$$\left.y(z_{i},t)\right]$	O (Δz )
3 puntos	${\textstyle {\frac {1}{\Delta z}}\left[y(z_{i+1},t)-\right.}$$\left.y(z_{i-1},t)\right]$	O (Δz2 )
3 puntos	${\textstyle {\frac {1}{2\Delta z}}\left[-3y(z_{i},t)+\right.}$$\left.4y(z_{i+/-1},t)-y(z_{i+/-2},t)\right]$	O (Δz2 )

A partir de las discretizaciones temporal y espacial anteriormente expuestas, la resolución numérica del problema (6) mediante (11) consiste en determinar la evolución del perfil discretizado en espacio y tiempo y_i,k  = y(z_i,t_k) , según:

$y_{i,k+1}=y_{i,k}+{\left(-1\right)}^{\nu }{\frac {H_{b}^{13/4}(t_{k})T_{b}^{3/2}(t_{k})}{K{\sigma }_{c}(z_{i})}}{\left({\partial }_{z}y_{i,k}\right)}^{-1}\sum _{j=0}^{(T/\Delta t)-1}Ip_{w,i,k}^{j}$

( 13)

desde i  = 1,…,N , y t_k  = kT para valores de k  = 0,…,k_max–1 , y para una discretización temporal de la onda de marea con paso Δt . Al ser T  <<  t_fin , puede tomarse el valor de k_max como el del entero más próximo al cociente t_fin/T sin que el error cometido afecte notablemente al alcance temporal t_fin de la simulación: k_maxT/t_fin ≈ 1.

3.3. Modelización del desprendimiento rocoso

A partir de (11) puede determinarse el volumen de material extraído por la abrasión marina del frente, por metro lineal de costa, entre 2 instantes de tiempo separados δt. Teniendo en cuenta el desplazamiento del perfil en el notch , único tramo del acantilado afectado por la erosión, tras un número de ciclos Δk de marea de periodo T:

$V_{e}(\delta t)=\int _{w_{t}^{B}-D_{b}}^{w_{t}^{P}}\left|y(z,t)-\right.$$\left.y(z,t+\delta t)\right|dz$

( 14)

donde ${\textstyle w_{t}^{P}}$ es el máximo nivel de la marea alta o pleamar y ${\textstyle w_{t}^{B}}$ es el mínimo nivel de la marea baja o bajamar en los Δk ciclos de marea, siendo D_b la profundidad a la que se produce la ruptura de la ola y a partir de la cual se desprecian los efectos de las fuerzas de arrastre del oleaje en profundidad. Teniendo en cuenta (13) y suponiendo que el intervalo de tiempo δt es proporcional a un número de ciclos de marea δt  =  ΔkT , la versión discreta de (14) puede escribirse como:

$V_{e}(\Delta kT)\approx \Delta z\sum _{i=Nmin}^{Nmax}\left|y_{i,k}-\right.$$\left.y_{i,k+\Delta k}\right|$

( 15)

siendo N_min y N_max los índices extremos de los puntos afectados (z_min y z_max ), pertenecientes al notch , por la acción marina en los Δk ciclos de marea transcurridos desde el ciclo k -ésimo.

Cuando se produce la desestabilización del acantilado en el momento t_r , el desprendimiento se encuentra limitado en su extensión lateral en la dirección OY por la superficie de ruptura cuya traza en el plano YOZ es la recta Π(z,t_r). Esta ruptura parte de un punto z_MAX-y(t_r) en el que se inicia y que se encuentra en la zona de acumulación de esfuerzos que corresponde con la posición más erosionada y profunda del notch[43]  and [44] . El buzamiento del plano de ruptura, inverso de la pendiente de Π(z,t_r) , se introduce como parámetro al modelo, al que se le incorpora un grado de incertidumbre sobre el que fluctuar. El instante en el que se produce la inestabilidad del frente del acantilado se considera descrito por un modelo estocástico, a partir de la distribución de probabilidad de los intervalos temporales δt_r entre 2 instantes consecutivos en los que se produce una ruptura. Las funciones de distribución tipo gamma, uniforme o pareto son las más habitualmente utilizadas [24]  and [25] .

Una vez producida la ruptura, el volumen de material desprendido por metro lineal de costa respecto al perfil del acantilado y hasta el punto y(z_out,t_r) , donde aflora en superficie la traza de la ruptura, es:

$V_{d}(t_{r})=\int _{z_{MAX-y}(t_{r})}^{z_{out}(t_{r})}\left|\Pi (z,t_{r})-\right.$$\left.y(z,t_{r})\right|dz$

( 16)

Para las escalas de tiempo manejadas en los pronósticos de recesión, el instante de la ruptura t_r puede aproximarse en unidades de tiempo por ciclo de marea: t_r  = k_rT   . Teniendo en cuenta la discretización espaciotemporal del perfil ${\textstyle y(z_{i},t_{r})\cong y(z_{i},k_{r}T)=}$$y_{i,k_{r}}$ , la versión discreta de (16) permite aproximar el volumen desprendido como:

$V_{d}(t_{r})\approx V_{d}(k_{r}T)=\Delta z\sum _{i=N{}'min}^{N{}'max}\left|\Pi (z_{i},k_{r}T)-\right.$$\left.y_{i,k_{r}}\right|$

( 17)

siendo N’_min y N’_max los índices extremos de los puntos z_MAX-y(t_r) y z_out(t_r) , respectivamente, sobre el perfil discretizado.

A partir de (14) y (16) puede determinarse el volumen total de material detrítico, por unidad lineal de costa, aportado entre 2 eventos de recesión costera:

$V_{total}(t_{r})=V_{e}(t_{r}-t_{r-1})+V_{d}(t_{r})$

( 18)

Es decir, la suma del volumen de roca desgastado por la acción marina entre la ruptura actual t_r y la inmediatamente anterior t_r-1 , y del volumen de la roca suprayacente al notch que se ha desprendido del frente del acantilado sobre el plano de ruptura en el momento actual. El volumen así calculado permite estimar la cantidad de sedimentos aportados por el acantilado al transporte marino que se produce por el oleaje y las corrientes, y que a su vez modifica la morfo-hidrodinámica de playas, puertos y otros acantilados colindantes. Este volumen puede ser totalmente arrastrado mar adentro o tan solo parcialmente, creando un repié de detritus o pequeño talud de acreción formado con materiales de resistencias menores que la roca sana del acantilado. En cualquiera de las 2 situaciones la geometría del acantilado se modifica. En el caso de que todo el material sea arrastrado instantáneamente tras la ruptura, la nueva configuración geométrica del acantilado se establece desplazando todos los puntos del frente afectados entre N’_min y N’_max sobre su cota z_i hasta la proyección sobre la traza del plano de ruptura Π(z_i,t_r) . Si se forma un pequeño talud frente al pie del acantilado, este se construye como un plano inclinado cuya pendiente viene dada por el ángulo de fricción interna del material desprendido y va desde la plataforma costera hasta el frente del acantilado reconfigurado geométricamente. Al material que forma este talud se le ha de asignar una resistencia σ_c menor que la del acantilado por su alto grado de fracturación y deterioro.

Todo el proceso de ruptura, desprendimiento y reconfiguración geométrica ocurre con unas duraciones que son instantáneas en comparación con los pasos de tiempo manejados en la integración temporal durante la cual se supone que el material caído ha sido retirado o se ha depositado al pie.

4. Estructura del código

Para obtener computacionalmente el desplazamiento del perfil mediante (13) y calcular los volúmenes de costa perdidos (18) durante el proceso de recesión, se ha desarrollado un código compuesto por 3 módulos organizados secuencialmente en su ejecución: primero, el de lectura de datos y preproceso; segundo, el de proceso computacional, y tercero, el de posproceso geomecánico y escritura de resultados. Este diseño permite que el cálculo se integre fácilmente en un esquema de análisis de sensibilidad paramétrico, y además es posible evaluar el riesgo y la incertidumbre en las variables más relevantes que intervienen en la simulación a largo plazo, como la resistencia del macizo rocoso (7) y la fuerza del oleaje (8), entre otras. La estructura general del código queda reflejada de manera simplificada en la figura 3 .

Figura 3. Diagrama de flujo del modelo desarrollado.

El bloque inicial se lleva a cabo en 2 etapas: primero lectura y luego adaptación de esta información para el cálculo. Los datos marinos —que pueden ser obtenidos de los institutos oceanográficos— son, si se considera una onda de marea, su amplitud A_m y su periodo T . Las condiciones del oleaje incidente se introducen mediante los valores de la altura de la ola en aguas profundas H₀ y su periodo T_b . La información geomecánica se incorpora mediante una columna litológica de resistencias a compresión σ_c(z) , referida la cota z al nivel medio del mar, además de la elección de la función de densidad de probabilidad paramétrica de los tiempos de cadencia de las rupturas. También se informa del plano de ruptura más probable según su buzamiento aparente respecto al eje OY ( fig. 1 b) del modelo. Un conjunto de puntos referidos al nivel medio del mar como origen son los datos necesarios para la geometría del perfil. También se han de introducir los datos correspondientes a las configuraciones del esquema numérico a utilizar, esto es, el tipo de integración (tabla 2 ) o derivación (tabla 3 ), pasos de tiempo Δt en el cálculo de la integral, tamaño de la discretización espacial Δz y alcance temporal de la simulación t_fin . La distribución de puntos así obtenida corresponde con la versión discreta de la curva y_i,0  = y(z_i,t₀) en el instante t₀  = 0 (condición inicial).

En el segundo bloque se realiza el cálculo, propiamente dicho, de la evolución del perfil según (13) sobre cada ciclo de marea. Los procesos implementados, como el cálculo del oleaje en la franja litoral y la rompiente [45] , la erosión producida (13), la probabilidad de producirse una ruptura y el tipo de ruptura, son calculados por pasos de tiempo T a lo largo del tiempo total de simulación [t₀,t_fin ]. El resultado final de una simulación son los efectos acumulados a lo largo del tiempo de dichos procesos sobre el acantilado.

Por último, en el tercer bloque se lleva a cabo un posproceso de los cálculos efectuados y se escriben los resultados. El posprocesado consiste en que en cada paso T se evalúan las condiciones de ruptura del frente erosionado, bien analizando el estado de equilibrio geomecánico mediante el factor de seguridad de ruptura en cuña [46] , o bien evaluando si la probabilidad del tiempo transcurrido tras la última ruptura es suficientemente alta para que haya una nueva ruptura. Si se produce una ruptura el posprocesado, suprime el material caído desde el plano de ruptura hasta el frente y reorganiza geométricamente el perfil del acantilado para entrar en un nuevo ciclo de marea. Así mismo, en caso de ruptura, evalúa los volúmenes que intervienen en (18).

5. Análisis del modelo

Para estudiar el comportamiento, precisar los mejores parámetros de discretización y establecer el método numérico más eficiente en la resolución del modelo de cálculo propuesto (13), se ha asumido que la oscilación de la marea w_t(t) es, de la forma más simple, una onda de amplitud A_m y periodo T constantes:

$w_{t}(t)=A_{m}(t)sen(2\pi t/T)$

( 19)

Dado que el objetivo es analizar su estabilidad frente a las discretizaciones espacial y temporal, se ha obtenido una solución exacta del modelo semidiscreto (11) utilizando (19) e incorporando p_w(z,t) explícitamente, lo que permite evaluar el error cometido en la aproximación numérica (13). La interpretación de este análisis es útil para controlar casos más complejos en los que no se disponga de solución exacta o para simulaciones a muy largo plazo en las que se desconoce el error, y es preciso tener un control de la precisión y de la exactitud de la solución.

5.1. Elección de la función de forma p_w(z,t)

Los estudios experimentales de Skafel y Skafel y Bishop [38]  and [39] permiten aproximar la función p_w(z,t) analizando las variaciones producidas por el oleaje en un modelo reducido de acantilado asentado sobre un perfil inicial de equilibrio z  = z(y)  = 0,18y^2/3 y sometido a una serie de olas seudoaleatorias en regímenes de vuelco y derrame. Sin embargo, dado que los experimentos se realizaron sobre rocas blandas, el rango de aplicación en función de los valores de resistencia a compresión σ_c de los materiales se debe limitar entre rocas muy blandas (1-5 MPa) y las moderadamente duras (50-100 MPa). Tomando sus resultados como tasas de erosión con respecto a la línea de costa (en la fig. 4 a se reproducen los datos de erosión medidos en el eje Y ), la distancia y con respecto a esta línea se ha convertido en profundidad de agua z , según el perfil inicial de equilibrio z(y) . A continuación se han normalizado las profundidades z respecto a la profundidad de rompiente D_b. Para calcular esta profundidad se ha supuesto que el valor crítico entre la altura de la ola rompiente H_b y la profundidad D_b a la que esta se produce es 0,78, si se consideran despreciables los efectos del periodo de onda, de la pendiente del fondo y del parámetro de semejanza del oleaje [47] . Para determinar H_b se ha utilizado la aproximación de [45] con la que, a partir de la altura de ola en aguas profundas H₀ y su periodo T₀ , se obtiene:

$D_{b}={\frac {1}{2}}g^{0.2}{\left(T_{0}H_{0}^{2}\right)}^{0.4}$

( 20)

donde g es la aceleración de la gravedad. Al efectuar Skafel y Bishop [38]  and [39] sus experimentos a diferentes profundidades, se aprecia en las erosiones resultantes una clara influencia de la pendiente del perfil z(y) que forma el fondo. Este efecto se ha corregido dividiendo las tasas de erosión por la pendiente del perfil según la distancia y a la línea de costa: z ′ = 0,12y^–1/3 . Los datos experimentales convertidos de las tasas de erosión (fig. 4 b) para cada tipo de oleaje son muy semejantes y presentan un máximo común en la relación z/D_b  = 0,25 entre las profundidades, lo que significa que para esas profundidades se produce una amplificación de la erosión sobre la roca. Este aumento se compensa con el paso del tiempo a causa del efecto que la pendiente local presenta en la ecuación (4), tendiendo a ser nulo a medida que la erosión progresa en la base del notch .

Figura 4. a) Datos experimentales en tasas de erosión (modificado de Skafel y Bishop [16] ). b) Calibración de la función del potencial erosivo pw(z,t) (línea continua) sobre datos experimentales reescalados para cada tipo de ola (puntos).

La semejanza entre los 2 experimentos ha permitido asignar una expresión matemática a la función de ponderación p_w(z,t) mediante una interpolación conjunta de los valores corregidos para cada ratio z/D_b . La función que se ha encontrado más apropiada para la interpolación en el intervalo [0,1] ha sido un polinomio racional de Chebyshev de orden 5/6:

$p_{w}(x,t)={\frac {p_{0}^{1}T_{0}(x)+p_{1}^{1}T_{1}(x)+p_{2}^{1}T_{2}(x)+p_{3}^{1}T_{3}(x)+p_{4}^{1}T_{4}(x)+p_{5}^{1}T_{5}(x)}{p_{0}^{2}T_{0}(x)+p_{1}^{2}T_{1}(x)+p_{2}^{2}T_{2}(x)+p_{3}^{2}T_{3}(x)+p_{4}^{2}T_{4}(x)+p_{5}^{2}T_{5}(x)+p_{6}^{2}T_{6}(x)}}$

( 21)

siendo T₀(x)  = 1 ; T₁(x)  = x ; T_n+1(x)  = 2xT_n(x)-T_n-1(x) ; los polinomios de Chebyshev de primer tipo son definidos entre [–1,1], y por tanto para poder utilizarlos de manera general es necesario su ajuste x  = (z/D_b-(0,5111604219506374))/0,5110522089614126 . Los 12 coeficientes utilizados para la construcción de este polinomio se resumen en la tabla 4 .

Tabla 4. Valores de los coeficientes ${\textstyle p_{j}^{i}}$ estimados del polinomio racional de Chebyshev (21) para los datos experimentales de Skafel y Bishop [39]

${\textstyle p_{j}^{i}}$	i = 1	i = 2
j = 0	−1,069838130676276	1.0
j = 1	−0,01221236978255141	−1,441825011828062
j = 2	0,4610744748787188	1,303412347447682
j = 3	0,07844716807628178	−0,9727210591773651
j = 4	0,652841133940385	0,5803220956908709
j = 5	−0,08926720346045524	−0,2877147080827893
j = 6	–	0,1702945908127132

El desarrollo y trasladado de p_w(z,t) sobre la onda (19) permite evaluar el potencial erosivo (escala de grises en la fig. 2 b) que se aplica en cada instante de tiempo entre [0,T ].

5.2. Solución exacta

Para este estudio se ha obtenido la solución exacta del modelo semidiscreto (11) introduciendo (19) y (21) para un frente rocoso inclinado 45° y 5 m de altura, con 2 m sobre el nivel del mar y 3 m sumergidos. Sobre este incide un oleaje estático y uniforme, de derecha a izquierda (ν  = 1), cuya altura en aguas profundas es de H₀  = 1,5 m y periodo T_b  = 1 s. La amplitud de la marea es constante A_m  = 1,46 m, y con periodo T  = 44850 s. Con el objetivo de amplificar los errores se ha supuesto una roca extremadamente blanda cuya resistencia a compresión es σ_c  = 1 MPa y para el valor de K  = 1, con el propósito de apreciar notablemente el efecto erosivo. En estas condiciones, mediante (11) es posible obtener el perfil erosionado tras T segundos calculando:

$y(z,T)=z-{\frac {H_{b}^{13/4}T_{b}^{3/2}}{K{\sigma }_{c}}}tan(\pi /4)\int _{0}^{T}p_{w}({\frac {A_{m}sen(2\pi t/T)-z}{D_{b}}})dt$

( 22)

Y, análogamente, tras 2T segundos o 2 ondas de marea:

$y(z,2T)=y(z,T)-{\frac {H_{b}^{13/4}T_{b}^{3/2}}{K{\sigma }_{c}}}{\left({\frac {\partial y(z,T)}{\partial z}}\right)}^{-1}\int _{0}^{T}p_{w}({\frac {A_{m}sen(2\pi t/T)-z}{D_{b}}})dt$

( 23)

En (23) hay que calcular el término de la pendiente local, dado que el perfil se diferencia de los 45° de partida al haber sido erosionado según (22) en el primer ciclo de marea. Los efectos de (22) y (23) se obtienen sobre el mismo intervalo de z  ∈ [–3,2], tal como se puede observar en la figura 5 a y b, respectivamente, trazadas con línea continua, tras el perfil de 45° inicial. A medida que el oleaje incidente ataca el frente tras cada onda de marea, el perfil se va desplazando en sentido negativo del eje OY .

Figura 5. Comparación entre la solución exacta y la simulación: a) tras una onda de marea y b) tras la segunda onda de marea, con un Δt  = 4,5 s. Los parámetros del modelo y discretización se detallan en el texto.

El motivo de escoger una pendiente de 45° para el frente del acantilado ha sido analizar el comportamiento de las simulaciones en las situaciones más desfavorables. Esto es porque al modelo se le ha introducido un límite de 50° para acotar los valores muy altos de la pendiente en (11) que se traducirían en erosiones anormales, como el máximo de erosión que se encuentra en z  = –1,75 m, tras la segunda marea (fig. 5 b). Por tanto, la utilización de 45° sitúa el caso muy cerca de su límite de validez, generando así una situación extrema para el estudio del error.

El modelo de erosión presentado representa un sistema retroalimentado en el que la historia de su evolución afecta significativamente a la solución en un instante. Sin embargo, gracias a que se toma una pendiente constante y conocida en todo el frente, por las expresiones (22) y (23) es posible estudiar separadamente y en 2 etapas el error en la solución numérica debido a las discretizaciones temporal y espacial. Así puede entenderse cómo se transmite el error entre 2 ciclos de marea. Para ello, en primer lugar, se estudian sobre (22) los efectos de la discretización temporal en la integral, de forma independiente a la espacial, ya que no afecta en el cálculo de las aproximaciones de la pendiente local por ser un valor exacto conocido (tan 45°) durante el primer ciclo de marea. Algunos de los primeros resultados en este sentido (fig. 5 a) muestran cómo la elección del paso de tiempo resulta crítica para la calidad de la solución aproximada. A continuación, y evaluada la sensibilidad de (22) a Δt , se estudia el efecto de la discretización espacial en (23), en la que interviene, por primera vez en el caso, alguna de las aproximaciones de la tabla 3 y cuyos efectos se aprecian en la figura 5 b, para un Δt  = 4,5 s, en diferentes Δz.

5.3. Efecto de la aproximación temporal (Δt)

Los errores cometidos en la aproximación temporal están sujetos al orden de precisión utilizado en las distintas técnicas de integración numérica (tabla 2 ). Para estudiar este tipo de error se ha realizado una serie de 1.100 simulaciones para diferentes valores de Δt , repartidos entre los 10⁴  s y los 4,5 s, y tomando los diferentes métodos de integración temporal de la tabla 2 , sobre un acantilado rocoso de perfil inicial inclinado con pendiente de 45°, con las mismas características enunciadas en el epígrafe anterior.

Los errores absolutos obtenidos para cada Δt tomado ( fig. 6 a) han sido calculados como la máxima diferencia entre la posición del frente de erosión obtenido de manera exacta y el obtenido de manera aproximada según el tipo de integración (tabla 2 ), para todos los z  ∈ [–3,2]. Los errores absolutos tienden a un valor umbral a medida que Δt disminuye, y por tanto se aumenta el número de pasos de tiempo necesarios para evaluar la erosión tras un periodo T . Al mismo tiempo que se utiliza una integración de mayor orden ( tabla 2 ), se obtiene una mejora notable en el valor de los errores absolutos para cada Δt . El rápido incremento en la amplitud del comportamiento oscilante frente a la solución límite surge antes (Δt más pequeños) en las integraciones de bajo orden con respecto a las de mayor orden, en un Δt  ≈900 s. Para los métodos de cuadratura gaussiana, el patrón oscilante aumenta su amplitud en el caso de la integración de 4 puntos (integral tipo IV) frente al de 7 puntos (integral tipo V), y con respecto a los demás tipos de integración utilizados su error absoluto siempre es más bajo.

Figura 6. Comparación entre la solución exacta y la aproximada con varios tipos de integración numérica. a) Errores absolutos. b) Errores relativos. Los parámetros del cálculo se detallan en el texto.

Los errores relativos (fig. 6 b) han sido calculados como la diferencia máxima en valor absoluto de la posición del frente de erosión en la zona de afección del oleaje obtenido de manera exacta, con respecto a la obtenida de manera aproximada, entre el valor exacto. Estos errores se comportan de manera similar a los absolutos. De la misma manera, es apreciable una tendencia de los valores del error relativo oscilatoria sobre un valor límite según aumenta el valor del paso de tiempo Δt . Las oscilaciones más pequeñas corresponden a los tipos de integración de mayor orden de precisión.

5.4. Efecto de la aproximación espacial (Δz)

Tras el primer ciclo de marea, para el cálculo de (23) entra en juego la aproximación numérica de la derivada local de la función y(z,t) , lo que da lugar a que se incorporen en la solución los errores debidos a la discretización espacial. Para su estudio, se parte del mismo caso que se utilizó en los 2 epígrafes anteriores, con el que se ha evaluado de forma aproximada la solución de (23), según un esquema de integración ( tabla 2 ) y derivación numérica (tabla 3 ), sobre un conjunto de 9.360 pares (Δt,Δz). Tras evaluar el valor de Δt en los que la solución se mantiene con una calidad relativamente aceptable, se han tomado 1.040 valores de Δt entre 4,5 y 900 s, valor a partir del cual se considera que las oscilaciones de este error comienzan a ser importantes. Por otro lado, se han escogido 9 valores de Δz comprendidos entre 0,01 y 0,45 m con los que analizar el comportamiento del error en 2T segundos y entre los que se encuentra el valor de la discretización propuesta por Walkden y Hall [40] (Δz  = 0,05 m). A partir de la solución numérica obtenida de (23) para cada par (Δt,Δz) se ha calculado el error absoluto (su logaritmo decimal se representa en las Figura 7 , Figura 8 , Figura 9 , Figura 10  and Figura 11 ) como la máxima diferencia en valor absoluto de la posición del frente de erosión en la zona de afección del oleaje obtenido de manera exacta, con respecto a la posición obtenida con el modelo.

Figura 7. Representación del logaritmo de los errores relativos cometidos en función de Δz y de Δt , para una simulación de 24,92 h e integración tipo I. La isolínea marcada corresponde a un error absoluto de 0,616 m.

Figura 8. Representación del logaritmo de los errores relativos cometidos en función de Δz y de Δt , para una simulación de 24,92 h e integración tipo II. La isolínea marcada corresponde a un error absoluto de 0,616 m.

Figura 9. Representación del logaritmo de los errores relativos cometidos en función de Δz y de Δt , para una simulación de 24,92 h e integración tipo III. La isolínea marcada corresponde a un error absoluto de 0,616 m.

Figura 10. Representación del logaritmo de los errores relativos cometidos en función de Δz y de Δt , para una simulación de 24,92 h e integración tipo IV. La isolínea marcada corresponde a un error absoluto de 0,616 m.

Figura 11. Representación del logaritmo de los errores relativos cometidos en función de Δz y de Δt , para una simulación de 24,92 h e integración tipo V. La isolínea marcada corresponde a un error absoluto de 0,616 m.

Del comportamiento que se aprecia en la representación de los errores obtenidos (Figura 7 , Figura 8 , Figura 9 , Figura 10  and Figura 11 ) es posible afirmar que, según aumenta el orden de precisión del método de integración, para ambos tipos de aproximación de la derivada, la estabilidad de la solución es mucho más alta, y por tanto es posible obtener una buena solución con mayores pasos de tiempo o Δt . Por ende, a medida que un método de integración de mayor orden es utilizado para el cálculo de la integral, es posible reducir el número de Δt necesarios para su obtención. También se puede observar una tendencia en las soluciones de O(Δz) , donde se tiene un óptimo en el valor del error relativo para Δz  = 0,05 m, lo que justifica el valor sugerido por Walkden y Hall [40] , e Δz  = 0,40 m; sin embargo, las soluciones para O(Δz²) presentan un óptimo en Δz  = 0,15 m. Por consiguiente, estos 3 valores de Δz son buenas opciones a priori para su uso en el modelo, según el esquema de derivación, y permiten reducir costes computacionales si se escoge un esquema de mayor orden en tiempo.

6. Conclusiones

A la hora de planificar y evaluar las posibles medidas de protección frente a peligros geológicos, es esencial disponer de la capacidad de pronosticar comportamientos y tendencias futuras de, en este caso, sistemas complejos litorales. Las revisiones efectuadas sobre los métodos para la obtención de predicciones verosímiles de las tasas de recesión han destacado la importancia de comprender el comportamiento morfodinámico de un acantilado [6] , [19] , [41] , [48] , [49]  and [50] para crear un modelo conceptual adecuado. El esquema morfodinámico desarrollado pone de manifiesto complicadas relaciones y retroalimentaciones que dan lugar a un comportamiento dinámico altamente complejo —y posiblemente caótico— en el sistema acantilado-playa, e implica que la recesión futura es fuertemente dependiente de qué ocurre en la actualidad y lo ocurrido en el pasado. Esta dependencia de las condiciones actuales y del pasado se traduce en un efecto memoria sobre todo el sistema. Si un sistema como este, además de poseer un efecto memoria, es no lineal, presenta un elevado número de grados de libertad y se encuentra sometido a excitaciones lentas, es muy probable que sea de comportamiento crítico autoorganizado [51] , mecanismo por el cual emerge la complejidad en la naturaleza. Por ello, la extrapolación de tasas de recesión históricas o incrementar el alcance de modelos deterministas, con el objeto de hacer pronósticos a largo plazo, puede resultar muy poco realista.

Con el propósito de simular la evolución temporal de una línea de costa rocosa se ha desarrollado un modelo de la erosión, como desplazamiento del frente interfase roca – agua/aire, que incorpora los principales procesos que desencadenan un cambio en el equilibrio del frente rocoso en la franja costera. Entre otros, los más importantes que se tienen en cuenta son la acción del oleaje, el efecto de las mareas, las propiedades geomecánicas de los macizos rocosos y la pendiente de la plataforma costera. En el modelo desarrollado se ha incorporado un conjunto de expresiones matemáticas relativamente simples de estos procesos, extraídas de la bibliografía consultada, y una ventajosa representación de sus retroalimentaciones a lo largo del tiempo. Este hecho lo diferencia de otros modelos mucho más complejos [16] , [17] , [23] , [24] , [25] , [34]  and [41] , que tratan los procesos de forma detallada y exhaustiva y, por tanto, demandan una amplia y precisa base de datos, lo que los hace muy sensibles a las condiciones iniciales, y que requieren un elevado tiempo de ejecución para proporcionar un resultado que, si bien es muy preciso, es poco exacto. El hecho de adoptar un modelo como el aquí propuesto proporciona la ventaja adicional de que, por precisión, exactitud y bajo coste, permite explorar la sensibilidad a las incertidumbres globalmente hasta la mesoescala.

Los resultados del modelo en las 2 primeras etapas de simulación (2 ciclos de marea) permiten afirmar que su comportamiento es mucho más robusto que los modelos más recientemente publicados, como los de Walkden y Dickson [26] o Walkden y Hall [27]  and [40] , para lo cual se ha generalizado el procedimiento de obtención de la función de forma p_w(z,t) a partir de estudios experimentales, siendo de mayor continuidad que otras obtenidas a trozos [40] , lo cual facilita su tratamiento numérico al describirse de manera continua el potencial erosivo en profundidad.

En los modelos [26] , [27]  and [40] previos no se estudian en ningún caso los errores cometidos en la resolución numérica del modelo presente. El análisis aquí realizado de los errores cometidos para las distintas aproximaciones implementadas, tanto para la integración como para la derivación numérica, permite concluir que el uso de estas herramientas matemáticas —incluso las de más bajo orden— da lugar a unos resultados comparables a los obtenidos de manera exacta. Por lo expuesto con anterioridad y a la vista de los resultados que aparecen en las Figura 7 , Figura 8 , Figura 9 , Figura 10  and Figura 11 , se puede concluir que el modelo presenta un mayor rango de estabilidad en la solución para métodos de mayor orden en la aproximación temporal Δt . Los resultados obtenidos en el análisis de los efectos de la aproximación espacial Δz , donde el valor de 0,05 m es un óptimo, se ven reafirmados por los datos presentados por dichos autores. Esto aporta un nuevo enfoque al problema, puesto que los modelos presentes en la bibliografía [26] , [27]  and [40] basan su cálculo en la aproximación de la integral con el menor orden (aproximación tipo trapecio). Es más, dichos modelos resuelven el problema con un paso de tiempo de una hora, o incluso con un golpe de marea, lo que puede provocar grandes errores en la solución del modelo y, por ende, formas de erosión poco realistas con abruptos cambios de dirección e irregularidades extremas no constatados hasta el momento y que podrían ser origen de falsas interpretaciones en la morfodinámica del sistema acantilado a largo plazo. Este cambio de concepto nos permite introducir un cálculo más preciso, sin que de ello se deriven mayores costes computacionales.

Además, dado que la estabilidad del modelo está fuertemente condicionada a la pendiente local de la plataforma litoral y puede tomar valores muy elevados en los tramos cuasi verticales entre el pie y la cresta del acantilado, con el modelo presentado y gracias a las aproximaciones utilizadas en la implementación y al estudio del error se ha logrado ampliar el valor límite de la pendiente de 10 a 50°. Esto mejora notablemente la limitación a 10° que han propuesto modelos más recientes [26]  and [27] y amplía considerablemente el rango de pendientes en el perfil del acantilado sobre el que es aplicable el modelo para que sean valores más próximos a los que se constatan en los casos reales [4] , [7] , [9] , [13] , [19] , [35]  and [50] . Para inclinaciones mayores de 50° la solución exacta se hace muy inestable tras un par de ciclos de marea. Esta inestabilidad se traduce en que el perfil se hace zigzagueante y con tasas de erosión extremadamente elevadas e irreales que habrá que limitar, para lo cual en versiones futuras se incorporarán otros métodos más sofisticados que suavicen este hecho.

Finalmente, se ha mostrado que es posible la mejora en la formulación y eficiencia de los modelos de erosión costera existentes. Los resultados aquí obtenidos sugieren la continuidad del desarrollo del modelo, con el claro objetivo de extenderlo a una escala temporal más amplia, junto con la implementación de algunos de los módulos representados en el esquema conceptual de su morfodinámica. Dichos módulos, en fase de desarrollo, son los considerados más influyentes en el fenómeno físico de la erosión [4] , [7] , [9] , [13] , [19] , [35]  and [50] , contemplan el cambio en el nivel del mar, la evaluación del factor de seguridad geomecánico del macizo rocoso y los posibles eventos tectónicos. Por otra parte, la calibración y la validación del modelo frente a datos reales son necesarias para demostrar que puede ser utilizado con confianza, y para ello se están testando diferentes localizaciones piloto aportadas por los proyectos que cofinancian este trabajo.

Agradecimientos

Este trabajo ha sido parcialmente financiado por el Instituto Geológico y Minero de España a través de los proyectos: Metodologías de estudio de la peligrosidad en acantilados costeros (MEPAC I), Metodología de estimación de la peligrosidad por caídas de bloques en acantilados costeros (MEPAC II), y Modelización y experimentación sobre peligrosidad por avenidas, movimientos del terreno y volcanismo (MODEX). Así mismo, la Universidad Politécnica de Madrid sufraga parte del desarrollo de este trabajo a través del Programa de Becas de Formación de Profesorado Universitario. Igualmente, los autores agradecen las aportaciones y las interesantes sugerencias de los 2 revisores anónimos que han contribuido indudablemente a la mejora del manuscrito.

Bibliografía

1. [1] M.J.F. Stive, S.G.F. Aarninkhof, L. Hamm; Modeling shoreface profile evolution; Mar. Geol., 126 (1) (1995), pp. 235–248
2. [2] A.S. Trenhaile; Modeling the role of weathering in shore platform development; Geomorphology, 94 (2008), pp. 24–39
3. [3] V.C. Lackhan; Advances in coastal modelling; Elsevier, The Netherlands (2003)
4. [4] E. Bird; Coastal geomorphology: an introduction; (2nd ed.)John Wiley & Sons, New York (2008)
5. [5] Eurosion; Living with coastal erosion in Europe: sediment and space for sustainability; Results from the Eurosion study, Office for Official Publications of the European Communities, Luxembourg (2004)
6. [6] J.D. Quinn, N.J. Rosser, W. Murphy, J.A. Lawrence; Identifying the behavioural characteristics of clay using intensive monitoring and geotechnical numerical modelling; Geomorphology, 120 (2010), pp. 107–122
7. [7] A.S. Trenhaile; The Geomorphology of Rock Coast; Clarendon Press, Oxford (1987)
8. [8] J.W. Kamphuis; Influence of sand or gravel on the erosion of cohesive sediment; J. Hydraul. Res., 28 (1) (1990), pp. 43–53
9. [9] T. Sunamura; The Geomorphology of Rocky Coasts; Wiley, Chichester (1992)
10. [10] J.W. Kamphuis; Recession rate of glacial till bluffs; ASCE J. Wtrwy., Port, Coast., and Oc. Engrg., 113 (1) (1987), pp. 60–73
11. [11] J.W. Kamphuis; Alongshore sediment transport rate; ASCE J. Wtrwy., Port, Coast., and Oc. Engrg., 117 (6) (1991), pp. 624–641
12. [12] P. Budetta, G. Galietta, A. Santo; A methodology for the study of the relation between coastal cliff erosion and the mechanical strength of soils and rock masses; Eng. Geol., 56 (2000), pp. 243–256
13. [13] A.S. Trenhaile; Modeling the erosion of cohesive clay coasts; Mar. Geol., 56 (1) (2009), pp. 59–72
14. [14] P.R. Wilcock, D.S. Miller, R.H. Shea, R.T. Kerkin; Frequency of effective wave activity and the recession of coastal bluffs: Calvert Cliffs, Maryland; J. Coastal. Res., 14 (1) (1998), pp. 256–268
15. [15] A. Mano, S. Suzuki; Erosion characteristics of sea cliffs on the Fukushima coast; Coast. Eng. J., 41 (1) (1999), pp. 43–63
16. [16] M. Crowell, B.C. Douglas, S.P. Leatherman; On forecasting future U.S. shoreline positions: a test of algorithms; J. Coastal Res., 13 (4) (1997), pp. 1245–1255
17. [17] S.M.N. Amin, R.G.D. Davidson-Arnott; A statistical analysis of controls on shoreline erosion rates, Lake Ontario; J. Coastal Res., 13 (4) (1997), pp. 1093–1101
18. [18] P. Milheiro-Oliveira, I.C. Meadowcroft; A methodology for modelling and prediction of coastal cliff recession; Proc. 4th Int. Conf. on Coastal Dynamics, ASCE, New York (2001)
19. [19] E.M. Lee; Coastal cliff recession risk: a simple judgement-based model; Q. J. Eng. Geol. Hydroge., 38 (2005), pp. 89–104
20. [20] S.P. Leatherman; Modelling shore response to sea-level rise on sedimentary coasts; Prog. Phys. Geo., 14 (1990), pp. 447–464
21. [21] A.S. Trenhaile; The effect of Holocene changes in relative sea level on the morphology of rocky coasts; Geomorphology, 114 (1-2) (2009), pp. 30–41
22. [22] M.J. Bray, Hooke J; Prediction of soft-cliff retreat with accelerating sea-level rise; J. Coastal Res., 13 (2) (1997), pp. 453–467
23. [23] J.W. Hall, I.C. Meadowcroft, E.M. Lee, P.H.A.J.M. Van Gelder; Stochastic simulation of episodic soft coastal cliff recession; Coast. Eng., 46 (3) (2002), pp. 159–174
24. [24] Furlan C; Hierarchical random effect models for coastal erosion of cliffs in the Holderness coast; Stat. Meth. Appl., 17 (2008), pp. 335–350
25. [25] P. Milheiro-Oliveira; Bayesian statistical methods for modelling and prediction of major landslides in coastal cliffs; Coast. Eng. J., 49 (1) (2007), pp. 45–61
26. [26] M.J.A. Walkden, M. Dickson; Equilibrium erosion on soft rock shores with shallow or absent beach under increased sea level rise; Mar. Geol., 251 (2008), pp. 75–84
27. [27] M.J.A. Walkden, J.W. Hall; A mesoscale predictive model of the evolution and management of a soft-rock coast; J. Coastal Res., 27 (3) (2011), pp. 529–543
28. [28] A. Bayram, M. Larson, H. Hanson; A new formula for the total longshore sediment transport rate; Coast. Eng., 54 (2007), pp. 700–710
29. [29] Z. YongHui, L. JinYou, L. HongZhi, W. JiaSheng, F. BeiLin, Y. ShiMing; Research on cohesive sediment erosion by flow: An overview; Sci. China Ser. E: Technol. Sci., 51 (11) (2008), pp. 2001–2012
30. [30] K. Appeaning Addo, M. Walkden, J.P. Mills; Detection, measurement and prediction of shoreline recession in Accra, Ghana, ISPRS J. Photogramm; Remote Sens., 63 (2008), pp. 543–558
31. [31] Boussinesq J; Théorie des ondes et des remous qui se propagent le long d’un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond; J. Math. Pure Appl. Deuxième Série, 17 (1872), pp. 55–108
32. [32] D.H. Peregrine; Long waves on a beach; J. Fluid Mech., 27 (1967), pp. 815–827
33. [33] C.C. Mei; The Applied Dynamics of Ocean Surface Waves; Advanced Series on Ocean Engineering, 1, World Scientific, Singapore (1989)
34. [34] A.P. Belov, P. Davies, A.T. Williams; Mathematical modeling of basal coastal cliff erosion in uniform strata: a theoretical approach; J. Geol., 107 (1998), pp. 99–109
35. [35] USACE; Shore Protection Manual; (4th ed)Department of the Army U. S. Corps of Engineers, Washington DC (1984)
36. [36] R.J. Dawson, M.E. Dickson, R.J. Nicholls, J.W. Hall, M.J.A. Walkden, P.K. Stansby, et al.; Integrated analysis of risks of coastal flooding and cliff erosion under scenarios of long term change; Climatic Change, 95 (1-2) (2009), pp. 249–288
37. [37] J.D. Quinn, L.K. Philip, W. Murphy; Understanding the recession of the Holderness Coast, east Yorkshire, UK: a new presentation of temporal and spatial patterns; Q. J. Eng. Geol. Hydroge., 42 (2009), pp. 165–178
38. [38] M.G. Skafel; Laboratory measurement of nearshore velocities and erosion of cohesive sediment (till) shorelines; Coast. Eng., 24 (1995), pp. 343–349
39. [39] M.G. Skafel, C.T. Bishop; Flume experiments on the erosion of till shores by waves; Coast. Eng., 23 (1994), pp. 329–348
40. [40] M.J.A. Walkden, J.W. Hall; A predictive mesoscale model of the erosion and profile development on soft rock shores; Coast. Eng., 52 (2005), pp. 535–563
41. [41] E.M. Lee, J.W. Hall, I.C. Meadowcroft; Coastal cliff recession: the use of probabilistic prediction methods; Geomorphology, 40 (2001), pp. 253–269
42. [42] J.D. Faires, R. Burden; Métodos Numéricos [traducido por Pedro José Paul Escolano]; Thomson, Madrid (2004)
43. [43] T. Kogure, Y. Matsukura; Critical notch depths for failure of coastal limestone cliffs: case study at Kuro-shima Island, Okinawa, Japan; Earth Surf. Proc. Land, 35 (2010), pp. 1044–1056
44. [44] G. Wolters, G. Müller; Effect of cliff shape on internal stresses and rock slope stability; J. Coastal Res., 24 (1) (2008), pp. 43–50
45. [45] P.D. Komar, M.K. Gaughan; Airy wave theory and Brecker height prediction; Proceedings 13th Coastal Engineering Conference (1972), pp. 405–418
46. [46] T. Kogure, H. Aoki, A. Maekado, T. Hirose, Y. Matsukura; Effect of the development of notches and tension cracks on instability of limestone coastal cliffs in the Ryukyus, Japan; Geomorphology, 80 (2006), pp. 236–244
47. [47] J.H. Balsillie, W.F. Tanner; Red flags on the beach, part II; J. Coastal Res., 16 (2000), pp. iii–x
48. [48] A.S. Trenhaile; Modeling cohesive clay coast evolution and response to climate change; Mar. Geol., 277 (2010), pp. 11–20
49. [49] B.D. Collins, N. Sitar; Processes of coastal bluff erosion in weakly lithified sands, Pacifica, California, USA; Geomorphology, 97 (2008), pp. 483–501
50. [50] R. Davidson-Arnott; An Introduction to Coastal Processes and Geomorphology; Cambridge University Press, New York (2010)
51. [51] P. Bak, C. Tang, K. Wiesenfeld; Self-organized criticality: an explanation of 1/f noise; Phys. Rev. Lett., 59 (1987), pp. 381–384