Envisat/ASAR Images for the Calibration of Wind Drag Action in the Doñana Wetlands 2D Hydrodynamic Model

Abstract

Doñana National Park wetlands are located in southwest Spain, on the right bank of the Guadalquivir River, near the Atlantic Ocean coast. The wetlands dry out completely every summer and progressively flood again throughout the fall and winter seasons. Given the flatness of Doñana’s topography, the wind drag action can induce the flooding or emergence of extensive areas, detectable in remote sensing images. Envisat/ASAR scenes acquired before and during strong and persistent wind episodes enabled the spatial delineation of the wind-induced water displacement. A two-dimensional hydrodynamic model of Doñana wetlands was built in 2006 with the aim to predict the effect of proposed hydrologic restoration actions within Doñana’s basin. In this work, on-site wind records and concurrent ASAR scenes are used for the calibration of the wind-drag modeling by assessing different formulations. Results show a good adjustment between the modeled and observed wind drag effect. Displacements of up to 2 km in the wind direction are satisfactorily reproduced by the hydrodynamic model, while including an atmospheric stability parameter led to no significant improvement of the results. Such evidence will contribute to a more accurate simulation of hypothetic or design scenarios, when no information is available for the atmospheric stability assessment.

1. Introduction

Doñana National Park (DNP) extents over 543 km² of the colmated Guadalquivir River paleo-estuary, in southwestern Spain. The area, considered the largest European habitat for migrating waterfowl, was declared a National Park in 1969, a Biosphere Reserve under the UNESCO Man and the Biosphere Programme in 1980, a Wetland of International Importance under the Ramsar Convention in 1982, a Special Protection Area under the European Union Directive on the Conservation of Wild Birds in 1988, and a UNESCO Natural World Heritage Site in 1994 [1]. About half of DNP’s extension is marshland area, which experiences annual cycles of inundation and depletion. While flooded, Doñana marshes are a vast area of shallow lakes of variable size and degree of spatial connectivity, mainly depending on the regional pluvial regime and on the time of year. The natural hydrological cycle of the marshes, which constitutes the strategic basis for the entire ecosystem functioning, has been threatened by a number of past and present human actions [2–4]. The alteration of tributary water courses (channels, deforestation, pollution) and the intensive agricultural occupation of the marshes during the 20th century have deeply altered their hydrological scheme and reduced their extension from about 1,800 to 300 km² [5]. At present, the water exchange between the marshes and the Guadalquivir River is managed by DNP’s authorities through a number of sluices located in an artificial 12 km long levee.

In order to restore the natural hydrological behavior of the marshes, Spanish authorities promoted in 1998 the project called Doñana 2005 [6]. The Institut Flumen, at the Universitat Politècnica de Catalunya and the International Center for Numerical Methods in Engineering, developed within the framework of this project a numerical model of DNP marshes in order to predict the hydrodynamic response of the wetland to the proposed restoration actions [7]. The rigorous calibration of the hydrodynamic model was assisted by field data from a network of in situ hydrometeorological gauging stations and by synoptic observations of the marshes through remote sensing images. Thus, Doñana images of the Advanced Synthetic Aperture Radar (ASAR) sensor on board on the Envisat satellite of the European Space Agency (ESA) [8], acquired since 2006 at different incidence angles and polarization configurations, enabled the calibration of different capabilities of the model, as well as the monitoring of the marshes’ flood extent evolution and the seasonal development of helophyte vegetation [9].

Given the flatness of Doñana’s topography, sustained winds can significantly modify the hydrodynamics and water extent of the marshes. Wind stress, as the vertical transfer of horizontal momentum between the atmosphere and surface, is considered a major driving force for shallow water circulation and mixing [10,11]. When sustained over a long period of time, it produces a setup of the water surface in the downwind direction. Various studies on shallow water bodies have documented the relationship between water level and local wind intensity, with short time lags of a few hours between wind speed changes and water level changes [12–14]. In the presence of low terrain slopes, as occurs in most wetlands and estuarine areas, this wind setup can eventually modify the inundation patterns [15]. If stratification occurs, persistent winds have also been identified as the cause of massive fish kills by upwellings of hypolimnetic water [16].

Current hydrodynamic modeling of wetlands and shallow reservoirs is normally taken into account using 2D depth averaged models based on either wastewater treatment or overseas experience [17–19]. Wind stress modeling over the wetland surface ranges from simple formulations using a constant wind drag coefficient [17] to more theoretically sounded complex formulations [20]. Despite wind being usually presumed as the major driver of water movement and studied phenomena such as seiche formation [17] or water turbidity [21], little discussion is found in the freshwater modeling literature about the choice of the wind stress formulation. However, wind drag modeling is still an ongoing discussion in coastal and marine engineering [22].

This paper describes the implementation of the wind stress action into the hydrodynamic model of Doñana marshes with the aid of Envisat/ASAR imagery. Two ASAR observation opportunities of an isolated water body are used to calibrate the observed wind-induced water displacement by choosing the most suitable wind stress formulation in the literature. The chosen formulation is also verified for another wind event observed in a different location of the marshes.

The paper is organized as follows. Second section describes the study area, Doñana marshes, and specifically the modeled ponds. Section 3 describes the study data, composed of radar imagery and field hydrometeorological data. Section 4 contains the methodology of flood mapping and hydrodynamic modeling. Section 5 presents and discusses the modeling results of the wind stress formulations tested at two observed events of wind-induced water displacement. Finally, Section 6 contains main conclusions and recommendations.

2. Study Area

Doñana marshes are composed of a variety of temporary shallow (<1.5 m depth) water bodies extended over an area of 298 km² on the right bank of the Guadalquivir River mouth (Figure 1). They are located at the former estuary of the river, which was filled by fine fluvial clayey sediment resulting in an extremely flat topography, with a maximum elevation difference of 2.50 m in its entire area. The main water inputs in the marshes come from direct precipitation and small and irregular northwestern tributaries, since they are disconnected from the Guadalquivir River and its tributaries by artificial structures. The climate is sub-humid Mediterranean with Atlantic influence [23], with mild winters (average temperature of 10 °C in January) and hot dry summers (24 °C in August). Mean annual precipitation is 537 mm [24], and average flooding periods occur from early autumn to late summer, when the marshes completely dry up [25,26]. On average, about 50% of the mean annual precipitation corresponds to the period from November to January, while less than 5% to the period from June to September [24]. However, the intra and inter-annual variability of the precipitation are high, resulting in flooding periods of variable magnitude and duration. Inside the marshes, small altitudinal gradients play an important role in defining the topographical elements (Figure 1), the water spatial and temporal distribution and the predominant vegetation communities.

In this work, two of the greater ponds (locally named “lucios”) of the marshes are selected for the calibration and verification of the wind stress effect on the water extent, both at the southern part of the marshes (Figure 1). First, a wind-induced water displacement event observed at Membrillo pond is selected for calibration purposes (named “Membrillo event”). The pond’s size is about 11.5 km² and has the main longitudinal axis oriented to the N-S direction. Using the chosen wind stress formula, one more event is simulated at the Ánsares pond (named “Ánsares event”), which has an area of about 7.5 km² and the main longitudinal axis at the WSW-ENE direction.

Most of the Membrillo and Ánsares ponds are formed by clayey bare soil, and virtually no vegetation develops in them at any time of the year [27]. Helophyte communities have colonized the pond’s slightly higher terrain, located at the western and northern ends of Ánsares and Membrillo, respectively. These communities are dominated by Scirpus maritimus (Castañuela) and Scirpus litorales (Bayunco). The Castañuela is an herb of height ranging between 0.6 and 1.0 m. The Bayunco is reed-like, and often taller than 1.0 m. Both species start emerging from the water surface towards the end of February. They experience rapid growth throughout the spring season and then dry out during the summer. Their brown stems stick out of the water surface when the marshes flood in autumn, but they progressively decay and sink by the spring.

3. Study Data

3.1. Envisat/ASAR Imagery

Three images of the ASAR sensor were used in this study for mapping the inundation in the DNP marshes. The ASAR sensor was installed on board of the Envisat satellite of the ESA [8] and was operational since the launch of the satellite in 2002, until its failure in 2012. ASAR imaged the Earth’s surface by emitting and receiving electromagnetic radiation at C-band (5.34 GHz) and could be operated in different measurement modes and at seven predetermined incidence angles or swaths, designated as IS1 to IS7. The appropriateness of the ASAR data for flood mapping and herbaceous wetlands monitoring has been widely reported. Some examples can be found in [28–34].

The images used in this study were acquired in Alternated Polarization and Image acquisition modes, and were received respectively in the form of Alternated Polarization Ellipsoid Geocoded (ASA_APG_1P) and Image Mode Precision (ASA_IMP_1P) products. Both products are multi-look, ground-range projected digital images, with a nominal resolution (range × azimuth) of 30 m × 30 m. In the Alternating Polarization mode, the ASAR sensor acquired data in two of the three polarization channels as follows: HH, VV and HV. In the HH and VV channels, also referred to as like-polar channels, the sensor antenna used respectively horizontally or vertically polarized radiation, and the same polarization was used for transmitting and receiving. HV is the cross-polar channel, where the transmitted and received pulses have different polarization. In the Image acquisition mode only one like-polar channel could be used, either HH or VV.

Table 1. Envisat/ASAR data.

Acquisition Date	ASAR Swath	Incidence Angle Range (°)	Acquisition Mode	Polarization Configuration
19 October 2006	IS1	15.0–22.9	Alternated polarization	HH/HV
1 March 2006	IS4	31.0–36.3	Image	HH
4 March 2006	IS2	19.2–26.7	Image	HH

summarizes the acquisition date, swath, incidence angle range, mode and polarization configuration for the images used in this study.

3.2. In Situ Hydrometeorological Measurements

A network of six hydrometeorological gauging stations was progressively installed in Doñana marshes between September 2004 and September 2007. These stations were operational until October 2011 and acquired 10 min records of water level, water temperature, conductivity, dissolved oxygen, precipitation, wind direction and velocity, air temperature and humidity, soil temperature, incident and reflected solar radiation and net all-wave radiation [35]. The stations used in this study are shown in Figure 1. The following field data were available for each case of study:

• - For the Membrillo event, measuring station D05 was still not operational. Wind speed and direction, air temperature and precipitation concurrent with the ASAR scenes were recorded at station D04. Water temperature was available at station D02.
• - For the Ánsares event, wind speed, wind direction and precipitation concurrent with the ASAR scene were recorded at station D05. Water level was available at station D02.

Horizontal wind velocity and direction were measured at 2 m height by 0513-5 Young helicoid anemometers with precisions of 0.3 m/s in the wind velocity and 3° in the wind direction. Air temperature was measured with two Vaisala HMP45C sensors with 0.2 °C accuracy. Water level was measured with GE Druck PTX 1830 pressure sensors located 10 cm over the ground, with accuracy of 1.6 mm for water depths less than 1 m. Water temperature was measured by duplicated Campbell Scientific PT-100 probes (0.1 °C accuracy) at 0, 0.5 and 1 m above the ground. The surface water temperature, when needed, was obtained by selecting the highest underwater record of temperature with the aid of the water level measurements. All measurements were taken at 1 min time intervals (except for the anemometers, at 2 s time intervals) and recorded at 10 min averages by a Campbell Scientific’s CR10 datalogger. Data were periodically downloaded by modem.

The wind drag effect was first observed in oscillations of the water level consistent with those of the wind speed and direction. Figure 2 shows hydrometeorological data from stations D02 and D04 (as named in Figure 1) from 20 February to 6 March 2006, coincident with the Membrillo event. The overall increasing tendency of water level responds to direct precipitation, since the marshes at that time were divided in several water bodies and the measurement points were disconnected from the northwestern tributaries. Short-term variations of the water level respond to the wind forcing, as they coincide with the wind velocity variations. Water level response to the sea breeze during calm days is of the order of 5 cm (e.g., 28 February–2 March in Figure 2), while wind events over 10 m/s can modify up to 10 cm the height of the water column (e.g., 26 February, 4 March in Figure 2). Further analyses of the temporal and spectral correlation of measured wind and water level at Doñana marshes are presented in [14].

The water level response to wind stress depends on location and wind direction. In this work, the latter is not supposed to change significantly across the marshes due to the practical absence of topographical obstacles. Locations shown in Figure 2 respond in an opposite way to the same wind forcing, depending on their position in their respective water body. For the same location, opposite responses are also found for N-NE and S-SW wind events (e.g., 28 February and 1 March, respectively).

4. Methodology

4.1. Flood Mapping from the ASAR Images

The ASAR scenes in Table 1 were received from the ESA as radar brightness and calibrated to backscattering coefficient according to [36]. The images were then precisely co-registered to the existing DTM. Approximately 80 ground control points were selected to perform the co-registration. Due to the terrain flatness, the ASAR images approximated fairly the projected geometry of the DTM, so a first degree polynomial and the nearest neighborhood methods were considered appropriate for the co-registration warping and resampling, respectively. The total root mean square error yielded values below 0.6 pixels for the co-registration of the three images.

4.1.1. ASAR Image Interpretation

Figure 3 shows the Membrillo pond on the calibrated ASAR images from 1 March and 4 March 2006. The Ánsares pond is depicted in the ASAR image from 19 October 2006 in Figure 4. The backscattering from these ponds was interpreted and related to their physical conditions at the time of the image acquisition, based on the concurrent wind records and the ASAR backscattering characterization of Doñana wetlands undertaken in [9]. This interpretation is explained as follows.

The ASAR image from 1 March 2006, in Figure 3a, was acquired under low wind conditions, as it can be appreciated in Figure 2. In this situation, open water surfaces are smooth and the radar beam is reflected specularly forward, away from the sensor. Hence, no return is received by the SAR antenna and the open water areas appear dark on the radar image. The brightest areas on the scene correspond to emerged land. Flooded surfaces with emerging sparse helophyte vegetation yield intermediate returns caused by the radar beam being partially reflected forward by the water surface and partially reflected back by direct and double bounce backscattering from the vertical stems [9].

The March 4th image, in Figure 3b, was acquired after a sustained, strong-wind episode, as revealed by Figure 2. The wind-roughened water surface induces high backscattering on like-polarized SAR channels at low incidence angles [9], such as the IS2 swath of the March 4th HH-polarized image. As a consequence, open water areas in the Doñana ponds and the Guadalquivir River are bright on this scene, although not as much as the emerged land. A noticeable dark area appears on the upwind side (southern end) of several ponds. This phenomenon had been previously observed on a number of low incidence angle ASAR images acquired on windy days, and the nature of the dark area was verified by ground truth data. When the water mass is pushed by the wind, a portion of the pond’s flooded area emerges at the upwind end. If the emerging surface consists of bare soil, this exhibits a smooth surface, formed by a puzzle of unconnected extremely shallow puddles, captured in Figure 5. The shallowness and reduced fetch of these paddles impede the formation of the waves responsible for the radar backscattering. Hence, bare soil areas emerging as a consequence of the wind drag action also cause forward specular reflection of the radar beam and appear dark on the ASAR images.

The area of flooded helophytes at the pond’s northern end in Figure 3b also shows lower backscattering than from the wind-roughened open water. As for the just unflooded bare soil, this phenomenon had already been observed on windy days. It was verified by field work that the flood surface among emerging vegetation remains smoother than in open water areas, and therefore appears darker than the latter on SAR images. This fact is attributable to the wind energy attenuation effect of the vegetation.

Figure 4a,b depicts the Ánsares pond in HV and HH polarization respectively, on the ASAR image from 19 October 2006. This scene was acquired after an episode of sustained southwest winds. The wind-roughened open water surface in Doñana ponds induces little depolarization of the radar signal, so backscatter from open water areas is low in the cross-polar channel and subsequently appears dark on Figure 4a. The HV return is most powerful on the emerged areas, were the volume scattering from the vegetation bushes causes depolarization in the radar signal. As for the March 1st image in Figure 3a, the HV power backscattered by flooded helophytes takes intermediate values, due to the combined volume backscattering from the vegetation and the specular reflection on the water surface.

In the HH channel, the emerged pastureland is again the cover type showing the largest backscattering. As already mentioned, the water surface roughness does induce significant radar returns in like-polar channels, and subsequently this target appears bright in Figure 4b. Although the open water shows slightly lower backscattering than the emerged land, both cover types are more easily separable in the HV channel, as can be observed by the simple visual contrast between them in Figure 4a.

Dark elongated features appear on the eastern end of the Ánsares pond, in Figure 4b. They are believed to correspond to patches of smoother water surface. This sort of dark patch with capricious geometries is often observed on wind-roughened water surfaces at low incidence angle SAR images. Some authors have related these patches to wind-induced upwellings in thermally stratified coastal water bodies [37]. However, the physical explanation to these features has not yet been studied in Doñana.

4.1.2. Image Filtering

Prior to the delineation of the ponds’ flooded area, the ASAR images were filtered out to smooth backscattering fluctuations owing to the speckle noise. The applied filtering method takes advantage of the tight relationship found between Doñana’s cover types and the terrain topography [26,27,38,39]. As a consequence of this relation, neighboring pixels at the same elevation are most likely to be of the same class and therefore constitute a stationary neighborhood. This peculiarity of Doñana’s landscape was utilized as follows:

At every image location P, the algorithm selects those connected pixels which are ±25 mm apart from the pixel to be filtered, P. This set of pixels forms the P’s irregular filtering neighborhood. The coefficient of variation (CV) is then used to assess the neighborhood stationarity. If compliant with the CV expected from the corresponding ASAR product ENL, pixel P’s filtered value is computed as the filtering region mean. Otherwise, the neighborhood is iteratively split in a direction orthogonal to the region’s main gradient, until the sub-part containing P meets the stationarity requirements. Further details on the filtering procedure and performance examples can be found in [38]. Figure 6 shows the filtering result for the scenes from 1 and 4 March 2006.

It is worth mentioning that the topography-based neighborhood sampling method is irrelevant for those pixels located in the middle of an open water patch. However, the flooded area perimeter still approximates the terrain contours, so topography-based neighborhoods are more likely to adjust to them than square filtering windows. As a result, the flood boundaries in DTM-filtered images are notably sharper than those obtained with fixed filtering window geometries. This comparison is discussed and sample images included in [38,40,41].

4.1.3. Flood Mapping

Regions representative of the target’s main land cover classes (ROIs) were selected on the three ASAR scenes. The selection was undertaken in the backscattering coefficient-terrain elevation space. The elevation data was introduced in order to increase the class separability. Figure 7 illustrates the ROIs selection for the March 4th image: points in Figure 7a’s plot depict the image pixels by their filtered backscattering coefficient and their terrain elevation. Points in the plot with the highest elevation and backscattering values correspond almost certainly to emerged land pixels. Ring A in Figure 7a indicates the selection of the emerged land ROIs. Ring B encircles a large set of pixels with intermediate elevation and radar returns which correspond to wind-roughened water areas. Most of the low backscattering points, in Ring C, correspond to just unflooded pixels at the southern end of the Membrillo pond. Pixel clusters A, B and C are depicted on the corresponding ASAR image in Figure 7b.

An iterative region growing procedure was then applied to extend the ROIs on the filtered image space. Firstly, the image pixels adjacent to the classes initial ROIs were selected. Next, the Mahalanobis distance [42] between every adjacent pixel and each class is computed. The elevation data was not used in the distance computation since it would somehow push the perimeter of the classes to be horizontal, and detection of the flood perimeter deviation from the horizontal plane was key for the wind drag action calibration. If the minimum Mahalanobis distance between an adjacent pixel and a sample class happens for the contiguous class, then the pixel is assigned to that class. When all assignments have been performed, the just classified pixels are adjoined to the sample regions. New adjacent pixels are determined and the whole process is repeated. The process ends when no new assignments are made in two consecutive iterations. Further details of the region growing procedure and sample results are included in [38].

Some flooded helophyte pixels were selected as just unflooded bare soil (region F in Figure 7b) at the ROIs selection stage. The reason is that backscattering and elevation values in flooded helophyte and recently unflooded bare soil areas overlap considerably in the backscattering-elevation space, and it was not possible to select two separate ROIs for them. As a consequence, the darkest flooded helophyte areas were classified as unflooded bare soil (blue class in Figure 7b), after the region growing procedure. Ancillary spatial data on Doñana marshes land cover types had to be applied to disentangle this confusion. The existing vegetation maps [43,44] indicate the presence of helophytes around region F. Such presence was corroborated by several ground truth campaigns. The unflooded bare soil class (blue in Figure 7b) resulting from the region growing procedure was intersected with the vegetation map. Those blue class areas at the north of the Membrillo pond coinciding with helophyte areas according to the vegetation map were reclassified as flooded helophytes.

It has to be said that the water/land discrimination accuracy would be greatly enhanced by the use of optical data. However, the wind-induced water displacement can change rapidly and capriciously with the wind speed and direction, so images taken only hours apart can show flooded areas shifts of hundreds of meters. Hence, an optical image appropriate to assist the water displacement delineation should have been acquired simultaneously to the ASAR observation, and such data were not available.

The above-described classification procedure was also applied to the ASAR images from 1 March 2006 and 19 October 2006 in order to determine the flood perimeter in the Membrillo and Ánsares ponds, respectively. For the October image, only the HV channel was used, since this polarization provides better discrimination than HH between wind-roughened open water and emerged land, as can be observed in Figure 4a. The March 1st scene was acquired under low wind conditions. As a consequence, the vast open water surface in the Membrillo pond was smooth and appears dark on the ASAR image. All this dark area was mapped as flooded by the classification algorithm. However, if waterlogged soil areas, as those depicted in Figure 5, were adjacent to the Membrillo open water, they would have not been discriminated from the open water and would be included as flooded in the resulting map.

Unfortunately, no ground truth data were collected coinciding with the ASAR acquisitions presented in this manuscript, so the flood mapping accuracy could not be assessed. However, the same flood mapping procedure was applied to other similar ASAR scenes of Doñana for which ground truth data were available. As reported in [38], the classification accuracy for those images was around 92%. This value provides an estimate of this study’s flood mapping accuracy.

The flood perimeters for the Membrillo and Ánsares main water bodies obtained by the above presented methodology are shown in conjunction with the hydrodynamic modeling results in Section 5.

4.2. Hydrodynamic Modeling

Within the frame of the Doñana 2005 restoration project, a two-dimension hydrodynamic model of Doñana marshes was developed [7,45]. Given the boundary conditions (topography, inflow hydrographs, precipitation, wind speed and direction, situation of the sluices communicating the marshes with the Guadalquivir river and water levels in the river) and the loss functions (mainly evapotranspiration), the hydrodynamic model simulates the water level and velocity evolution in the marshes. This model was built on the basis of the Iber software [46], which is a numerical tool for the simulation of unsteady open-channel turbulent flow and sediment transport in two dimensions. The water movement is computed by solving the 2D Shallow Water Equations or Saint Venant equations.

4.2.1. Numerical Scheme

The equations solved by the hydrodynamic model are the 2D Saint Venant equations (or two-dimensional shallow water equations), which can be written in conservative form as:

$$\frac{\partial }{\partial t}\mathbf{U}+\mathbf{\nabla }\mathbf{F}=\mathbf{H}$$

(1)

where U_2D is the conserved variables vector, F_2D is the flux tensor and H_2D is the source term. In this work the wind stresses have been included as part of the source term:

$$\mathbf{U}=\left(\begin{array}{c}h\\ \mathit{hu}\\ \mathit{hv}\end{array}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathbf{F}=\left[\begin{array}{cc}\mathit{hu}& \mathit{hv}\\ {\mathit{hu}}^2+\frac{{\mathit{gh}}^2}{2}& \mathit{huv}\\ \mathit{huv}& {\mathit{uv}}^2+\frac{{\mathit{gh}}^2}{2}\end{array}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathbf{H}=\left(\begin{array}{c}i\\ \mathit{gh}\left(S_{\mathit{ox}}-S_{\mathit{fx}}\right)-{\tau }_x/\rho \\ \mathit{gh}\left(S_{\mathit{oy}}-S_{\mathit{fy}}\right)-{\tau }_y/\rho \end{array}\right)$$

(2)

Here, h is the depth, u and v the two horizontal depth averaged velocity components, g gravity, S_o the bed slope, S_f the friction slope and τ the friction on the water surface induced by wind.

A high resolution finite volume scheme, based on the Godunov scheme and the Roe Approximate Riemann solver [47], and using the Finite Volume Method is used. Special care is taken in the discretization of geometric source terms following the ideas of improved treatment developed by Vazquez–Cendon [48].

The resulting scheme is explicit, and is the base of Iber, a numerical model for simulating turbulent free surface unsteady flow and environmental processes in rivers and estuaries [49]. The ranges of application of Iber cover river hydrodynamics, dam-break simulation, flood zones evaluation, sediment transport calculation and wave flow in estuaries. The model can be run over an irregular mesh of triangles, quadrilaterals or a combination of both of them [50]. The wind effect is considered in the source term H in Equation (2) through τ_x and τ_y, which are the decomposition in the x and y directions of the wind stress over the water surface τ.

4.2.2. Wind Stress Formulations

The wind stress over the surface is expressed in terms of the bulk drag relation:

$$\tau =\rho u_{*}^2=\rho C_D{U}_z^2$$

(3)

where ρ is the mass density of air (about 1.225 kg·m⁻³) and μ_* is the friction velocity, which is replaced by the measured wind velocity at a reference height U_z and a dimensionless drag coefficient C_D corresponding to that height z. The wind velocity profile (∂U/∂z) depends on surface roughness and atmospheric stability, whose influence is usually included in C_D. A variety of wind drag formulations have appeared in the past literature, most of them referred to a reference height of 10 m. In this work, we test five different schemes in order to achieve a satisfactory balance between field data requirements and quality of the modeling results.

Table 2. Characteristics of the wind drag formulations of the simulations performed at the Membrillo event.

Simulation	Wind Profile	Drag Formulation	Atmospheric Stability	Low Wind Speed Correction
S1	Exponential	Constant	Neutral	No
S2	Exponential	Van Dorn	Neutral	No
S3	Logarithmic	Charnock	Non-neutral	No
S4	Logarithmic	Charnock	Neutral	No
S5	Logarithmic (exponential for low wind)	Charnock, (Wuest and Lorke for low wind)	Neutral	Yes

summarizes the main characteristics of the formulation used in each simulation, named S1 to S5.

The simplest option for the wind drag modeling is performed in simulation S1 as a constant value of C_D = 2.6 × 10⁻³ given by Munk [51] for wind speeds measured at 10 m. Simulation S2 uses Van Dorn [52] wind drag expression:

$$C_D=1.2\times {10}^{-3}\hspace{0.17em}\hspace{0.17em}\text{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{U}_{10}\le 5.6\hspace{0.17em}\text{m}/\text{s}$$

(4)

$$C_D=1.2\times {10}^{-3}+2.25\times {10}^{-3}\left(1-\frac{5.6}{U_{10}}\right)\hspace{0.17em}\hspace{0.17em}\text{for}\hspace{0.17em}\hspace{0.17em}{U}_{10}>5.6\hspace{0.17em}\text{m}/\text{s}$$

(5)

In both cases, the wind speed at 10 m height is obtained from the power law wind profile:

$$\frac{U_{z1}}{U_{z2}}={\left(\frac{z_1}{z_2}\right)}^m$$

(6)

where m = 0.11 is the empirically derived coefficient obtained by Hsu [53] for neutral atmospheric stability over the sea.

Taking into account the more realistic logarithmic wind profile above the surface and non-neutral atmospheric stability conditions, C_D can be written as [54],

$$C_D=k^2{\left[\text{ln}\left(\frac{z}{z_0}\right)-{\psi }_m\left(\frac{z}{L}\right)\right]}^{-2}$$

(7)

where k = 0.41 is von Kármán’s constant. The surface roughness length z₀ is defined according to Charnock [55],

$$z_0=\alpha \frac{u_{*}}{g}$$

(8)

where α = 0.0185 is the Charnock constant [54] and g is the gravitational constant. ψ_m(z/L) is the Monin–Obukhov similarity function for momentum flux, which is neglected under neutral stability conditions (simulations S4 and S5 in Table 2). Under non-neutral stability conditions (simulation S3), this function is defined according to Hsu [56]:

$${\psi }_m\left(\frac{z}{L}\right)=-5\frac{z}{L}\hspace{0.17em}\hspace{0.17em}\text{for} \text{stable} \text{conditions} (z/L\ge 0)$$

(9)

$${\psi }_m\left(\frac{z}{L}\right)=1.0496{\left(-\frac{z}{L}\right)}^{0.4591}\text{for} \text{unstable} \text{conditions} (z/L<0)$$

(10)

According to Donelan [57], the stability parameter (z/L) is related to the bulk Richardson number Ri_B, such that for unstable conditions (i.e., T_w > T_a)

$$\frac{z}{L}=7.6{\mathit{Ri}}_B$$

(11)

and for unstable conditions (i.e., T_a > T_w)

$$\frac{z}{L}=6.0{\mathit{Ri}}_B$$

(12)

T_a and T_w being the air and water temperatures, respectively, and [54,56]

$${\mathit{Ri}}_B=\frac{\mathit{gz}\left(T_a-T_w\right)}{\left(T_a+273.2\right)U_z^2}$$

(13)

Equations (3), (7) and (8) define an implicit scheme which is solved for C_D. The solution is shown in Figure 8 for the range of observed values of U₂ and ψ_m(z/L) during the period of study.

Furthermore, there is experimental evidence of the increase of C_D over water surfaces for weak winds under 5 m/s [58–60]. As about 60% of U₁₀ estimations by Equation (6) during the modeling period fall under such wind conditions, simulation S5 performs the empirical expression proposed by Wuest and Lorke [61]:

$$C_D\approx 0.0044U_{10}^{-1.15}\hspace{0.17em}\hspace{0.17em}\text{for} U_{10}<5\hspace{0.17em}\text{m}/\text{s}$$

(14)

4.2.3. Model Setup

The numerical model was based on a digital terrain model (DTM) of Doñana marshes built from a point elevation data set acquired on a LIDAR survey flight in September 2002. It has 2.00 m² and 0.15 m planimetric and elevation resolution, respectively [62]. The numerical computations were carried out on two irregular meshes of triangular elements (TIN), one for each of the two modeling areas of Ánsares and Membrillo drawn in Figure 1, obtained from the DTM.

Regarding the initial conditions, water volume for the Membrillo case of study was obtained via GIS techniques from the flood maps derived from the ASAR scenes together with the DTM of the marshes. For such purpose, the terrain elevation was subtracted from the water surface to obtain the observed water volume. In the Ánsares event, where water depth data were available, the water volume was chosen so as to coincide with the observed initial water depth at the measuring location. At initial time, a horizontal water surface is not expected to coincide with ASAR observations. For that reason, simulations started 4.5 days before the first ASAR observation of the Membrillo event, and 1.5 days before the unique ASAR observation of the Ánsares event.

Bottom friction of each case, implemented as a manning value, was chosen according to field observations by the authors. For the Membrillo event, a spatial distinction was made in the modeling area, setting a value of 0.125 at the scarcely vegetated northern third of the pond and a value of 0.030 at the non-vegetated rest of the pond. The Ánsares case was performed as a completely non-vegetated area with a manning value of 0.023.

Wind stress formulations are implemented in the numerical model without spatial variations or dependence on the presence of vegetation. Further details of the model setup of each case of study are shown in

Table 3. Characteristics of the model setup for each case of study.

Case of Study	Computational Mesh (TIN)				Initial Volume (m3)	Initial Time	Final Time
	Num. of Elements	Num. of Nodes	Min. Element Size (m)	Max. Element Size (m)
Membrillo	52,348	27,381	12.5	200	1,956,337	25 February 2006 0000 h	4 March 2006 1,200 h
Ánsares	7,352	3,804	13.6	398	550,000/600,000	18 October 2006 0000 h	20 October 2006 0000 h

.

4.2.4. Model Performance

As the wind stress exerted over the water surface is theoretically dependent on atmospheric conditions, a greater number of observed events would be necessary to perform a desirable calibration-verification procedure. For this reason, wind stress formulations are tested but calibration of their parameters is not performed.

In order to quantify the degree of agreement between simulated and ASAR imagery based flood maps, the fitting index in Equation (15) has been used [63–65]:

$$F=\frac{A}{B+C-A}$$

(15)

Parameter B in Equation (15) represents the observed flooded area on the ASAR image. C corresponds to the simulated flooded area at the time of the image acquisition. A is the overlapping between areas B and C.

5. Results and Discussion

5.1. Comparison of Wind Stress Formulations, Membrillo Event

Figure 9a shows the time series of the bulk wind drag coefficient C_D of the simulations performed in the Membrillo case of study, whose wind conditions were shown in Figure 2. To make them comparable, S1 and S2 coefficients have been adapted to a reference height of 2 m using the power law relation for neutral atmospheric stability Equation (6). The constant value of C_D in simulation S1 is noticeably higher than the rest, being from 1.74 to 2.54 times the average value of the other computed C_D coefficients during the period of study. The low wind speed correction used in simulation S5 greatly modifies the time evolution of C_D with respect to those time-dependant C_D of simulations S2, S3 and S4. However, as the highest values of C_D in S5 correspond to the lowest wind intensities, this difference is barely transferred to the effective wind stress exerted over the water surface, as shown in Figure 9b. According to the resulting wind drag coefficient, the highest wind stress is computed in simulation S1. Its intensity is only reached by the Van Dorn formulation (S2) in the case of strong winds, as can be seen on days 25, 26 February and 4 March. During most of the simulation period, wind stress computed in simulation S2 remains at a medium magnitude between S1 and the Charnock-based formulations of simulations S3, S4 and S5.

In order to better distinguish the effect of the different formulations performed, Figure 9c shows the accumulated wind stress over the water surface. At the end of the study period, simulations S1 and S2 have performed an accumulated stress of 1.89 and 1.32 times the wind stress obtained by the Charnock formulation for neutral stability (S4), respectively.

Modifications of the Charnock formulation computed in simulations S3 and S5 are intended to account for the physically sounded effect of the atmospheric stability and for previously reported experimental evidences of higher wind stresses at lower wind speeds, respectively. However, their accumulated effect over the 7.5 days of simulation only differs from the wind stress computed in simulation S4 by ±5%, as shown in Figure 9c. From a practical point of view, the influence of such modifications in this work seems to be negligible, taking into account their higher field data and computational requirements.

5.2. Modeling Results

5.2.1. Selection of the Wind Stress Formulation (Membrillo Event)

In order to assess the model capability to reproduce the wind-driven water movement and to choose the most suitable wind stress formulation, the flood extent of each simulation coincident with the ASAR scenes were compared with the flood contours obtained from the ASAR images. F parameters obtained for these comparisons are shown in

Table 4. F parameters obtained for the water extent of each simulation at the ASAR observation times.

	ASAR 01Mar06	ASAR 04Mar06	F med
S1	0.78	0.74	0.76
S2	0.78	0.77	0.77
S3	0.78	0.76	0.77
S4	0.79	0.76	0.77
S5	0.79	0.76	0.77
F med	0.78	0.76	0.77

.

As expected from the results presented in the previous section, little differences could be found between the Charnock-based formulations S3, S4 and S5, with a maximum model performance difference of 1% (F parameters in Table 4). Little but negligible differences were also found between them and the results of simulation S2, which despite having lower data requirements performs 1% better the flood extent of the second ASAR scene. Simulation S1 produced similar results for the first ASAR scene (1 March 2006), but overestimated the wind drag action at the time coincident with the second ASAR image (4 March 2006, the analogous situation of Figure 10d). Overall mean F values are 0.77 for simulations S2 to S5 and 0.76 for simulation S1. Tested formulations performed flood extent on 1 March 2006 between 1% and 4% better than on 4 March 2006, mainly due to the excess of water displacement to the vegetated northern part of the pond. This points to the necessity of introducing wind stress spatial variations or an excess of resistance to flood (manning value) if modeling more vegetated areas during the growing season.

By balancing field data requirements, computational efforts and modeling results, wind stress formulation of simulation S2 was chosen for further analyses and recommendations. Using the same data requirements as the formulation of S1, S2 performs an improvement up to 3% on the flood area estimation (Table 4). Charnock-based formulations S3 to S5 require more field data and computational time though they do not perform more than 1% improvement of the water extent estimation.

5.2.2. Analysis of the Membrillo Event

Figure 10 shows four time steps of the chosen model configuration: initial condition of horizontal water surface (Figure 9a); maximum flooded extension towards the south (Figure 10b); simulation time step coincident with the first ASAR image (Figure 10c); and simulation time step coincident with the second ASAR image (Figure 10d).

The initial time of all simulations is settled 4.5 days before the acquisition of the first ASAR image, a period of time expected to be long enough to minimize the influence of the initial condition of horizontal water surface. About 1.5 days before the first image, the wind starts to blow consistently from the north, with maximum intensity about 23 h before the time of the image acquisition (Figure 2). As shown in Figure 10b, the maximum simulated flooded extension towards the south, on 28 February 2006 1,500 h, approximates the flood contour of 1 March 2006. However, during the last 17 h before the ASAR image is taken, the wind remains calm (mostly under 1.5 m/s, Figure 2), and the water body tends to recover the horizontal surface. The lack of coincidence of the southern flooded area at the very same time of the first ASAR image, in Figure 10c, may respond to the different sensitivity of the techniques applied to determine the water extent. In the absence of vegetation, the algorithms applied to the ASAR image may be classified as a water-saturated sediment that was flooded a few hours before, whereas the numerical model identifies this area as dry through a wet-dry threshold of 1 mm of water height over a smooth computational mesh.

After the time of occurrence of the first image, a S-SW wind event starts to develop with maximum intensity approximately at the time of acquisition of the second ASAR image, as is represented in Figure 10d. Both flooded areas compared in this figure coincide to a great extent. However, the numerically simulated water area seems to be moved towards the west with respect to the flood contour obtained from the ASAR image. The cause of this difference might be twofold: on the one hand, small altimetric errors in the marshes’ DTM or in the computational TIN can change significantly the simulated water extent in this area, as the slopes of the terrain are minima; on the other hand, the wind measurement station D04 is located about 14 km from the modeling area, and slight rotations of the wind direction can occur inside the marshes despite their flatness, which are not accounted for in the numerical model.

5.2.3. Analysis of the Ánsares Event

The time series of measured wind velocity, wind direction and water depth during the Ánsares event on 18–19 October 2006 are shown in Figure 11, together with the estimated water depth at the measurement point by two simulations with different initial water volume. According to the measured water depth at the initial time, the event is firstly modeled with 5.5 × 10⁵ m³ (green line in Figure 11). In this configuration, the average error for estimated water depth is −8.8 × 10⁻³ m and the mean absolute error is 10.8 × 10⁻³ m. The model satisfactorily reproduces the observed water depth variations on 18 October, reproducing a decrease of up to 6 cm about 3 h after the wind event reaches its maximum velocity of 9.4 m/s. However, initial water volume seems to be insufficient, since water depth on 19 October with lower wind velocities tends to underestimate the measurements by 15 × 10⁻³ m. Moderate winds measured at the beginning of the simulation could explain a non-horizontal water surface and a higher water volume. Thus, a second simulation with an initial volume of 6.0 × 10⁵ m³ is also presented (blue line in Figure 11). In this case, average error for estimated water depth is 5.3 × 10⁻³ m and mean absolute error is 7.6 × 10⁻³ m. Although the results of 18 October are less satisfactory (mean error of 11.7 × 10⁻³ m), the modeled water depth on 19 October coincides to a great extent with measurements (mean error of −1.1 × 10⁻³ m). At the time of radar image acquisition, the error of water depth estimation is reduced to −3.3 × 10⁻³ m.

Figure 12 shows the water level and the water extent obtained at three different time steps of the Ánsares event simulation using an initial water volume of 6 × 10⁵ m³, together with the location of the D02 gauging station and the water perimeter obtained from the ASAR scene. Figure 12a shows the horizontal water surface at the beginning of the simulation (18 October 2006, 0000 h). The maximum water displacement towards the east due to the wind event on 18 October is presented in Figure 12b. At that time, the water body performs about 1 km displacement from the initial position and develops a mean surface slope of 0.1%, reaching the eastern water perimeter estimated by remote sensing techniques. Similar to the Membrillo event studied previously, the ASAR image is taken when the water body is recovering its horizontal surface (Figure 12c). The lack of coincidence of the eastern perimeter at the time of the ASAR image is also explained by the different sensitivity of the techniques applied to determine the water extent. This assumption is reinforced by results of simulated water depth in Figure 11.

Taking into account that both Figure 12b and Figure 12c are jointly compared to the water extent obtained by the ASAR scene, simulated flood extents at those time steps were superimposed to evaluate the model performance. F values obtained were 0.87 for both initial water volumes of 5.5 and 6 × 10⁵ m³.

6. Conclusions

This paper describes the implementation of the wind stress action into the two-dimensional hydrodynamic model of Doñana marshes with the aid of Envisat/ASAR imagery. Given the flatness of the Doñana’s terrain, remote sensing images provided valuable spatial data on the wind-induced water bodies’ displacement. This information could barely be achieved by point field measurements.

Five wind stress formulations in the literature were introduced in the marshes’ hydrodynamic model. On-site wind records were then used to simulate the water wind drag effect in an isolated water body, the Membrillo pond, by using the five formulations. Two ASAR observations of the Membrillo pond flooded area during the same wind episode were used to assess the modeling results.

The Van Dorn’s expression [52] and the power law wind profile for neutral stability conditions were selected as the most appropriate wind stress modeling technique. This formulation requires less field data and computational effort than the Charnock-based expressions, and yielded satisfactory results with an overall F coefficient value of 0.77 [63–65]. The chosen formulation was verified by simulating a different wind event at the Ánsares pond, for which on-site wind records and ASAR imagery were also available.

Wind-induced water displacements up to 2 km were satisfactorily reproduced by the selected formulation, and water level variations were modeled with an average absolute error under 11 mm at 10 min calculation time steps. The maximum water displacements took place in time intervals of 12 to 18 h, in response to moderate and rather frequent wind velocities, below 15 m/s. In the presented study cases, the inclusion of atmospheric stability considerations and the increase of the wind drag coefficient at low wind speeds led to no significant improvement of the modeling results, with an increment of the F value under 1%.

Future lines of research will aim at mapping vegetation biomass by using polarimetric SAR imagery. These vegetation data will be implemented in the hydrodynamic model to spatially reproduce bedload resistance to flow, surface wind shear-stress reduction and evapotranspiration water losses. Further efforts will focus on modeling and calibrating water quality parameters of ecological relevance, such as wind-induced sediment resuspension and turbidity.