El trabajo se focaliza en un segmento costero de 1800 m en la playa de l’Auir en el extremo norte del municipio de Gandia, al sur de la provincia de Valencia (Figura 1). Esta playa se encuentra al sur de la desembocadura del río Xúquer (Cullera) que actúa como la principal fuente de aprovisionamiento sedimentario de todo este tramo (Sanjaume y Pardo-Pascual, 2004). Este hecho ha determinado que, históricamente, todo este tramo tuviera un claro superávit sedimentario que favoreció el desarrollo de una amplia restinga y la formación de varias alineaciones dunares, alguna de las cuales ha sido históricamente desmantelada para ser ocupadas por cultivos (Sanjaume y Pardo-Pascual, 2004). Actualmente, la alineación dunar existente se extiende 7,3 km de los que aquí se estudiarán los últimos 1,8 km, entre la desembocadura del río Vaca, en el límite municipal con Xeraco, y la playa Norte de Gandia, cuyo borde meridional queda flanqueado por el puerto de este municipio. Los diques del puerto de Gandia tienen un papel clave en la evolución de este tramo dado que atrapan y paralizan el sedimento que las corrientes longitudinales costeras arrastran debido a que todo el sector presenta una clara deriva litoral hacia el sur (Sanjaume y Pardo-Pascual, 2004). Este hecho provoca que La playa de l’Auir se encuentre en un tramo claramente acumulativo y esté formada por arena fina y dorada, con partículas de entre 0,2 y 0,3 mm de tamaño (Cabezas-Rabadán et al., 2019b). Debido a la riqueza ecológica y ecosistémica, aparecen catalogados distintos hábitats de interés comunitario. Este espacio está incluido en la Red Natura 2000 como Lugar de Interés Comunitario (LIC), declarado hoy en día como Zona de Especial Conservación (ZEC) ‘ES5233038 Dunes de la Safor’, y en el catálogo de Zonas Húmedas de la Comunidad Valenciana (https://mediambient.gva.es).

Las dunas se extienden a lo largo de la totalidad del frente costero de la zona estudiada, presentando una anchura predominante de 60-70 metros en la mayoría de los tramos. El perfil longitudinal es regular y continuo, salvo en los accesos peatonales, donde el cordón dunar se reduce. En el perfil transversal, las dunas se alinean paralelas a la orilla en un patrón típico de sistemas dunares costeros progradantes, desde las formaciones más recientes hasta las más antiguas. Este sistema dunar continúa hacia el norte hasta la playa de Tavernes de la Valldigna, donde su anchura disminuye hasta desaparecer en algunos puntos, mientras que hacia el sur la ocupación humana lo ha hecho desaparecer.

En este tramo costero, la marea astronómica rara vez supera los 20 cm de altura (Cabezas-Rabadán et al., 2019b), y el oleaje suele ser moderado, con olas cuya altura significante (Hs) no sobrepasa habitualmente los 1,5 m (Figura 2). Sin embargo, en condiciones meteorológicas adversas, el nivel del mar puede experimentar aumentos significativos. Este es el caso de episodios como el temporal Gloria (20 de enero de 2020), en el que la boya oceanográfica del puerto de Valencia registró una altura de ola significante de 8,67 m (Pardo-Pascual et al., 2022). Los oleajes más intensos provienen principalmente del ENE y, en menor medida, del NE (Figura 3A). Esta dinámica, combinada con la orientación de esta costa, ha favorecido históricamente el transporte sedimentario longitudinal hacia el sur (Pardo-Pascual et al., 2021). Por su parte, el régimen de vientos es bidireccional, alternando principalmente entre los vientos del W (los más mayoritarios) y del WNW y aquellos procedentes del ESE (Figura 3B).

Durante el periodo analizado, diferentes temporales costeros afectaron a la zona de estudio. Los eventos más destacados fueron identificados de acuerdo con la altura de ola significante, aplicando el método de picos sobre un umbral (Peaks over threshold, POT) empleado habitualmente para la caracterización de temporales costeros (e.g., Dorsch et al., 2008; Ferreira y Soares, 1998). Así, se estableció un umbral de 1,42 m basado en una probabilidad de ocurrencia del 5%, considerando aquellos eventos con una duración mínima de 12 horas y, al menos, cinco días de separación entre ellos. Con el objetivo de poder comparar la magnitud y el potencial erosivo de los eventos, el índice de Dolan y Davis (1994) fue calculado para cuantificar la intensidad relativa de las tormentas combinando la máxima altura de ola significante (Hs_max) y su duración (D) mediante la fórmula: I=Hsmax2⋅D, donde I se expresa en m²h. La duración se calculó sumando las horas en las que la altura de ola superó el umbral, y Hs_max correspondió al valor máximo registrado durante cada temporal (Tabla 1).

Para este estudio, se ha partido de un conjunto de nueve ficheros de nubes de puntos 3D (Tabla 2), tres de ellos obtenidos mediante tecnología LiDAR aeroportado (2008, 2009 y 2015) y seis mediante técnicas de fotogrametría aplicadas sobre fotografías aéreas adquiridas en diferentes años (2017-2022) por el Instituto Cartográfico Valenciano (ICV, https://icv.gva.es/es/). Las nubes de puntos correspondientes a los años 2009 y 2015, pertenecientes al Plan Nacional de Ortofotografía aérea (PNOA), se adquirieron a través del Centro de Descargas del Instituto Geográfico Nacional (IGN). Por su parte, los modelos Digitales de Superficie fueron derivados en el seno del grupo de Cartografía GeoAmbiental y Teledetección (CGAT-UPV) mediante fotogrametría a partir de fotografías aéreas, siguiendo la metodología descrita por Almonacid-Caballer et al. (2025) que permiten obtener precisiones altimétricas comparables a las de los vuelos LiDAR (0,15 - 0,30 m). El proceso parte de fotografías aéreas de 25 cm/píxel, con un solapamiento longitudinal mayor al 60% y del 30% en sentido transversal. Asimismo, se dispone de puntos de control mediante la identificación de puntos invariables en el tiempo sobre las ortofotografías anuales del ICV. La cota de esos puntos de control se obtiene de un MDS generado a partir de los datos LiDAR de 2015. El proceso fotogramétrico posterior se realiza con la herramienta Metashape Photoscan, incluyendo las etapas de enmascaramiento de las zonas acuáticas, la orientación de las imágenes, su georreferenciación mediante los puntos de control, la generación de los mapas de profundidad y la triangulación ulterior para obtener cada MDS anual. Tal y como demuestran Almonacid-Caballer et al. (2025), estos modelos aún pueden ser refinados para mejorar su intercomparabilidad. Para ello, el resultado obtenido se compara frente a posiciones invariantes, como son puntos del terreno situados sobre los ejes de las carreteras de la zona y se cuantifica el sesgo (residuo) de cada posición. Tomando todos los residuos medidos se genera un modelo de residuos con el que se corrige el modelo fotogramétrico original, lográndose así una modelización tridimensional más precisa.

La medida de error de cada uno de los MDS (definida como el error medio cuadrático en la Tabla 2) se ha definido comparando la altura de puntos invariantes (no usados para la modelización) medidos en fuentes documentales diferentes. En el caso de los modelos procedentes de LiDAR, la comparación se ha hecho frente a datos tomados mediante técnicas DGNSS. En el caso de los modelos procedentes de técnicas fotogramétricas se ha comparado frente al MDS procedente de los datos LiDAR de 2015.

Para calcular los MDS a partir de los datos LiDAR, se aplicó en primer lugar un proceso de filtrado para eliminar puntos con valores de elevación anómalos. A partir de los ficheros filtrados se calcularon MDS con un tamaño de píxel 1×1 m. Estos procesos se realizaron a partir del software FUSION (McGaughey, 2023) integrado en QGIS (3.34.11). Las zonas de estudio, frente de duna y zona de playa emergida se caracterizan por presentar escasa vegetación y de poca altura por lo que no es necesario aplicar ninguna operación adicional a la nube de puntos 3D para extraer los puntos pertenecientes al suelo.

Para analizar los cambios que ha habido en la zona del frente de duna se siguieron dos metodologías: una volumétrica y una morfométrica. La primera de ellas consistió en calcular diferencias de volúmenes de arena entre los periodos 2009-2008 (2009 menos 2008), 2015-2009, 2017-2015, 2019-2018, 2020-2019, 2021-2020, 2022-2021 y 2022-2008. Para ello, en primer lugar, fue necesario delimitar la zona del frente dunar (Figura 1). Su borde más alejado a la costa se obtuvo a partir de la digitalización de la línea del frente de duna del año 2008, año en el que se observó esta línea más tierra adentro. Esta línea coincidiría con la divisoria del frente de duna de este año. Por su parte, el borde del polígono más cercano a la orilla se digitalizó siguiendo los postes de protección de la duna existentes en la zona de estudio (Figura 4A). En algunos casos, se observó que la duna había sobrepasado estos postes. En estos casos, se modificó esta línea midiendo puntos adicionales que se tomaron como puntos pertenecientes a la base de la duna con un equipo GNSS Leica Zeno FLX 100. Para el cálculo de los volúmenes, se calcularon las diferencias entre los diferentes MDS en cada periodo de análisis. Estas imágenes se convirtieron a formato punto y mediante una selección por atributos se seleccionaron los puntos con diferencias mayores y menores a 0. Estas operaciones se realizaron con el programa ArcGIS 10.2 (ESRI, 2013). Para cada conjunto de puntos, se sumaron los valores y, teniendo en cuenta el tamaño de píxel (1×1 m), el resultado de la suma fue el valor del volumen. Es importante indicar que el cálculo de volumen no ha tenido en consideración las variaciones que quedan dentro del límite de detección (Tabla 2) ya que no es posible asegurar que se haya producido cambios. Dicho límite ha sido definido por la precisión altimétrica expresada mediante el error medio cuadrático, RMSE (coincidente con la desviación típica en los modelos gestionados) y, dado que se comparan dos fechas concretas, siempre se toma el valor de error más alto de las dos, para que la estimación del cambio volumétrico sea más robusta. En el caso de los MDS procedentes de fotogrametría, el error está obtenido frente al MDS derivado de LiDAR de 2015 (que también tiene su propio error), por lo que se ha hecho una nueva estimación de error como la suma cuadrática del RMSE de este modelo (LiDAR) y el de cada MDS fotogramétrico. Posteriormente se calcularon las diferencias entre los MDS de todos los periodos analizados y se generó un mapa con esta información. Al mismo tiempo, también se calcularon diferencias de volúmenes en la zona de playa emergida para los mismos periodos de análisis que el caso anterior. Para ello, fue necesario definir un nuevo polígono. La delimitación del borde más alejado a la costa coincidió con la línea de la zona del frente dunar, definida a partir de los postes de protección y los datos medidos en campo. Para definir el borde más cercano al mar se superpusieron todos los MDS de la serie de datos y se digitalizó la línea más cercana al mar en la que había información de elevación en todos los MDS del periodo de análisis. En este caso no se consideró el MDS de 2008 dado que no se disponía de información de la parte norte de playa.

En cuanto a la segunda metodología utilizada para analizar los cambios en el frente de duna, se generaron 18 perfiles homogéneamente distribuidos (aproximadamente cada 100 m) a lo largo de la playa (Figura 1) sobre los que se realizó una caracterización morfométrica. Estos perfiles tuvieron como puntos de inicio las intersecciones con un camino paralelo a la costa (tierra adentro fuera de la zona dunar) estable durante el periodo analizado, y como punto final el límite de la zona del frente de duna más cercano a la orilla. En algunos casos, esta distancia de 100 m de separación fue diferente para evitar zonas de paso y áreas en las que el frente de duna no presentaba una estructura bien definida. Para cada perfil y periodo de análisis se definió la elevación del frente de duna con diferentes estadísticos (promedio, máximo, mínimo y la desviación estándar), a la vez que se midió la distancia desde el inicio de cada perfil hasta el punto de máxima cota (Figura 4B). Para cada periodo definido por un par de fechas (MDSs) a analizar, y para cada perfil, se calculó la diferencia de estos dos parámetros: cambios de cota de la cresta de la duna, y desplazamiento de posición de la cresta del frente dunar.

A pesar de que la acción de los temporales costeros puede verse reflejada a corto plazo y de forma directa sobre el pie de la duna (Van Gent et al., 2008), en este trabajo no se ha realizado una caracterización empleando el pie de la duna. La subjetividad e inexactitud en la identificación de parámetros como el pie de duna ha sido señalada previamente en la literatura (Wernette et al., 2018), y también se ha evidenciado en l’Auir debido a la presencia de un perfil sin cambios abruptos de pendientes.

Con respecto a la influencia de las características oceanográficas sobre el comportamiento dunar, se analizó la relación entre el oleaje y la morfología de la duna. Para ello, una vez definidos los parámetros morfométricos y volúmenes asociados a cada MDS y sus variaciones a lo largo del tiempo se exploró su relación con los episodios de oleaje más significativos registrados a lo largo del periodo de estudio (ver Tabla 1).

La disponibilidad de datos tridimensionales a lo largo del tiempo ha permitido analizar la variación de volumen de sedimento en la zona del frente de duna de la playa de l’Auir. Así, entre los años 2008 y 2022 ha habido un incremento de 25095 m³ de arena (Tabla 3, Figura 5a). Cabe resaltar un incremento de arena en la mayoría de los periodos analizados exceptuando los años 2009-2008, con un descenso de –919 m³ (Figura 5b); entre los años 2021 y 2020, de –172 m³ y; en el último periodo analizado, 2021-2022, con un valor de –81 m³ (Figura 5h). En el periodo de análisis se puede destacar el incremento de volumen entre los años 2009 y 2015 (12 005 m³) (Figura 5c), lo que representa casi la mitad del volumen total acumulado. También fueron significativos los incrementos de volumen entre 2015 y 2017 con 6262 m³ (Figura 5d) y entre los años 2019 y 2020 con un total de 3300 m³ (Figura 5f).

Al analizar los cambios de volumen de arena en la zona de playa, se ha observado también, aunque con menor magnitud, un aumento de volumen de 2646 m³ en el periodo analizado (Tabla 3). En esta zona de análisis, se han alternado periodos de pérdida y de ganancia de arena. Así, los periodos 2015-2009, 2018-2017 y 2021-2019 presentaron valores positivos de volúmenes, siendo la acumulación más significativa la correspondiente al primero de estos periodos con 4883 m³. Este incremento máximo coincidió con el observado en la zona del frente de duna. Los periodos 2017-2015, 2019-2018 y 2022-2021 fueron periodos en los que se observaron valores negativos de volúmenes. En estos casos hay que destacar el periodo 2019-2018 con un valor de –1715 m³. Cuando se consideran ambas zonas de análisis se observa que, en 4 de los periodos analizados hay una coincidencia de comportamiento del cambio positivo. En contraste, los periodos 2017-2015, 2019-2018 y 2021-2020 presentaron volúmenes con signo opuesto en ambas zonas.

La identificación de la cresta del frente dunar en diferentes perfiles perpendiculares a la línea de costa ha permitido caracterizar su cota y su posición a lo largo del tiempo, así como cuantificar sus cambios entre los diferentes periodos analizados. En cuanto a la cota de la cresta del frente dunar se puede observar cómo, en promedio, ésta presenta valores cercanos a los 4 m (Tabla 4). El año 2008 presentó el menor promedio de cotas (3,9 m), así como la cota mínima por perfil (3,4 m) y una variabilidad entre perfiles relativamente baja (desviación estándar de 0,2 m). Por su parte, el año 2009, el frente dunar presenta valores similares tanto de cota promedio (4,0 m) como de variabilidad entre perfiles (0,2 m). En contraste, 2015 es el año con el promedio de cota más elevado (4,3 m), las cotas máximas y mínimas más altas (4,8 m y 3,8 m respectivamente) y la mayor diversidad entre perfiles (0,3 m). En 2017 el frente dunar presenta tanto una cota promedio (4,1 m) como valores máximos y mínimos absolutos (4,5 m y 3,6 m, respectivamente) ligeramente inferiores. A partir de este momento las elevaciones de la cresta de la duna parecen estabilizarse, con escasas variaciones tanto en sus valores promedio (manteniéndose entre los 4,2 m y 4,2 m) como con lo que respecta a los máximos y mínimos absolutos.

Además de los valores de elevación del frente dunar, en este apartado se van a describir los resultados obtenidos con relación a las diferencias entre las elevaciones máximas del frente dunar en los perfiles digitalizados considerando los periodos de análisis: 2022-2008, 2009-2008, 2015-2009, 2017-2015, 2018-2017, 2019-2018, 2021-2019, 2021-2020, 2022-2021. Además, se van a comentar los resultados de los incrementos de distancias entre los puntos iniciales de un conjunto de perfiles al punto de elevación máxima del frente de dunar para cada periodo analizado.

En el caso de las diferencias de elevaciones máximas, se obtuvo una diferencia media de 0,3 m en el periodo 2008-2022 (Tabla 5). Entre los periodos analizados, el valor máximo fue 0,3 m en el periodo 2015-2009, lo que coincide con el mismo periodo de mayor volumen de arena descrito en los párrafos anteriores. Todas las diferencias fueron positivas exceptuando los periodos 2017-2015 y 2022-2021 con valores de –0,2 m y –0,1 m, respectivamente. También se puede destacar que el promedio de las diferencias fue cercano a 0 entre los años 2018 y 2021, lo que indica cierta estabilidad en las cotas máximas del frente de dunas durante este periodo. Durante este intervalo de tiempo también se observó que las desviaciones típicas de las diferencias de cotas fueron las menores de todo el periodo analizado.

Por lo que se refiere al análisis de la posición del frente de duna, hay que destacar un desplazamiento medio hacia el mar de 1,8 m en el periodo analizado (Tabla 6). Este incremento no fue uniforme en toda la serie de datos analizados. Así, los periodos con mayor desplazamiento de la línea del frente de duna tuvieron lugar durante los periodos 2009-2008 y 2015-2009 con valores de 0,6 m y 1,4 m, respectivamente. Los periodos 2019-2020 y 2020-2021 presentaron también incrementos positivos, aunque de menor magnitud con valores de 0,3 m y 0,2 m, respectivamente. También se observaron valores negativos en la serie de datos analizados destacando los periodos 2018-2017 y 2022-2021, ambos con valores de –0,3 m. Con relación al análisis de los valores de desviaciones típicas se puede resaltar que los primeros tres periodos de análisis (2009-2008, 2015-2009, 2017-2015) fueron los que presentaron valores más altos frente al resto de años analizados lo que puede indicar entre qué periodos los cambios tuvieron una variación de mayor calado.

Con respecto a la influencia de las características oceanográficas sobre el comportamiento dunar, se ha analizado la relación entre el oleaje (ver Figura 2, Tabla 1) y la morfología de la duna. Así, se puede observar que la intensidad media anual de las tormentas registradas a lo largo del periodo de estudio describe un patrón que guarda cierta correspondencia con la cota media de la cresta del frente dunar (Figura 6). De este modo, en 2015 tiene lugar el crecimiento más abrupto de la cresta de la duna registrado a lo largo de la serie, el cual coincide con un episodio de baja intensidad en los temporales (265,9 m²h). Este crecimiento contrasta con el descenso registrado en 2017 (alcanzando los 4,11 m), asociado a un aumento de intensidad en los temporales (532,1 m²h). En 2018 tiene lugar una nueva elevación de la cresta, que alcanza los 4,2 m, asociada nuevamente a un periodo de baja intensidad (312,7 m²h). Por el contrario, los valores de intensidad relativamente elevada registrados a partir de 2019 (por encima de los 600 m²h y con un máximo de 748,7 m²h en 2020) no se traducen en descensos significativos de la elevación de la cresta, aspecto que será comentado en la discusión.