Evidence of Hydroperiod Shortening in a Preserved System of Temporary Ponds

Based on field data simultaneous with Landsat overpasses from six different dates, we developed a robust linear model to predict subpixel fractions of water cover. The model was applied to a time series of 174 Landsat TM and ETM+ images to reconstruct the flooding regime of a system of small temporary ponds and to study their spatio-temporal changes in a 23-year period. We tried to differentiate natural fluctuations from trends in hydrologic variables (i.e., hydroperiod shortening) that may threaten the preservation of the system. Although medium-resolution remote sensing data have rarely been applied to the monitoring of small-sized wetlands, this study evidences its utility to understand the hydrology of temporary ponds at a local scale. We show that the temporary ponds in Doñana National Park constitute a large and heterogeneous system with high intra and inter-annual variability. We also evidence that the conservation value of this ecosystem is threatened by the observed tendency to shorter annual hydroperiods in recent years, probably due to aquifer exploitation. This system of temporary ponds deserves special attention for the high density and heterogeneity of natural ponds, not common in Europe. For this reason, management decisions to avoid its destruction or degradation are critical.


Introduction
Conservation strategies should take into account the spatial distribution of habitats because it conditions the distribution and dynamics of the species associated with them [1][2][3][4].However, the spatial distribution of habitats that change over time (fluctuating habitats), such as temporary ponds, is not constant.Thereby, conservation programs focused in such fluctuating ecosystems should also incorporate the monitoring of temporal changes in their distribution and extent.Temporary ponds are shallow water bodies with a recurrent annual dry phase that remain flooded long enough to allow the development of (semi) aquatic vegetation and animal communities [5][6][7].They present wide variability in filling onset and duration, depending on rainfall input and pattern [8], evapotranspiration and groundwater inputs.Temporary ponds are the main breeding habitat of many aquatic invertebrate [5,6,9] and amphibian species [1,7,10].Due to their conservation value, temporary ponds are a priority habitat under the European Union Habitats Directive, code 3170 [11].However, pond ecosystems are threatened worldwide by their drastic reduction in number [12].For this reason, long-term studies recording the temporal changes in the number, distribution, extent and duration of flooding (hydroperiod) of temporary pond ecosystems are necessary, although they have been scarce until now.Moreover, one of the most important assessments would be the discrimination of natural fluctuations from trends of habitat degradation or disappearance, which may threaten the population stability of species associated with them.
Remote sensing is a valuable tool for monitoring ecosystems and, in particular, of wetlands [13].Satellite images have been used to map large permanent wetlands [see 14 for a review] and, occasionally, to reconstruct their temporal dynamics over the past decades [15][16][17].However, its application on seasonally flooded systems has been scarce [18][19][20][21], and, in particular, for small-sized wetlands [22] due to the difficulty of identifying waterbodies if below pixel size, e.g., [23].High spatial resolution remote sensing may be an alternative for the delineation of small-sized wetlands.However, its application has been scarce [19,24] due to the high cost of such imagery.Landsat images may provide adequate data for long-term reconstruction of wetland dynamics due to its temporal resolution (16 days repeat cycle) and the continuity of the missions (operational for more than 30 years).However, its spatial resolution (30 m pixel) may compromise the identification of small sized ponds.A retrospective study of the presence of water in temporary ponds can constitute a high quality data set to analyse the relationship between rainfall and pond dynamics.The relation between rainfall and water presence, as derived from remote sensing data, has been previously evaluated in mediumsized and large temporary wetlands, e.g., [15,16], but such an approach has never been attempted for small temporary ponds.
Doñana National Park is one of the most important wetlands in southern Europe.It is located in southwestern Spain, eastwards from the mouth of the Guadalquivir River.Apart from the Guadalquivier river marshes, it preserves a large network of small temporary ponds [25,26].These temporary ponds are a critical habitat for many species of aquatic flora and fauna: macrophytes [27], invertebrates [28][29][30], and amphibians [10,31,32].Although Doñana temporary ponds are mainly fed by rainfall, its hydrology is also directly dependent on groundwater.The aquifer system in the Doñana region consists of two units: a relatively thick unconfined aquifer overlying a lower and more permeable aquifer [33].It is mainly the unconfined aquifer, which has a shallow water-table and several flow systems, that feeds temporary ponds [34].Groundwater feeding is complex in some ponds due to temporal variation in the connection to aquifer flow systems [35].The water-table follows the hydraulic head gradient, which roughly corresponds to topographic elevation.So, the water-table is close to the surface everywhere except beneath the dunes ridges [33].Doñana temporary ponds have no direct connection to the sea except through airborne salt deposition [35].Serrano and Zunzunegui [36], in a study conducted in six ponds, indicated that some of them may be under threat, since their hydrologic regime may have been damaged by aquifer water extraction by a nearby tourist resort [see also 37].Scientifically-based management of Doñana National Park requires information on historical flooding patterns, both spatially and temporally, and their relationship with natural fluctuations and anthropogenic modifications [38].
In this study, we investigate the applicability of Landsat satellite images to assess the spatio-temporal dynamics of a large system of temporary ponds (more than 800 waterbodies).We assess (i) the usefulness of Landsat imagery for identifying the presence of water in small temporary ponds; (ii) the intra-annual pattern of flooding and drying out; (iii) the spatial variation, i.e., differences in flooded area or duration of flooding, among ecosections within the National Park; (iv) the variations in water presence during a 23-year period, mainly to discriminate natural fluctuations from trends of habitat degradation (i.e., reduction in flooded area or annual hydroperiod shortening); and (v) the lagged response of temporary ponds extent to rainfall.

Study Area
The study was conducted in the Doñana Biological Reserve (6,794 ha), within Doñana National Park [see 39 for a geo-morphological description], in southwestern Spain (Figure 1a).Two main ecosystems are differentiated within Doñana Biological Reserve: the aeolian sands and the marshes.The object of the present study was the system of temporary ponds located on the aeolian sands.In this ecosystem, dominant vegetation is Mediterranean scrub (Halimio halimifolii-Stauracanthetum genistoides and Erico scopariae-Ulicetum australis as defined by Rivas-Martínez et al. [40]) with small patches of pine (Pinus pinea L.) and juniper forests (Juniperus phoenicea L. subsp.turbinata (Guss.)Nyman).We followed the ecosection zonation proposed by Montes et al. [41], who differentiated six ecosections (areas with similar ecosystems) within the aeolian sands (see Supplementary Material Table 1 for details).We further divided the largest ecosection into its northern and southern part, thus identifying seven zones in the aeolian sands of the study area.Those ecosections showed different temporary pond density and different geomorphologic, stratigraphic and hydrodynamic characteristics (Figure 1b and, for details on characteristics, Supplementary Material).
A total of 883 temporary ponds and two permanent ponds were delineated in this study area at a time of a large flooding event in April 2004 using a hyperspectral remote sensing image with high spatial resolution (5 m pixel size) [see 42 for details on cartographic methods].A detailed description of these temporary ponds can be found elsewhere [43].Filling dates and flooding duration of temporary ponds are highly variable among years depending on rainfall pattern and quantity.The wet season may extend from September to May and the dry season from June to August.In this study, we consider a hydrological year to be the period from September 1st to August 31st, which corresponds to the traditional farming year in many Mediterranean countries.Mean

Pre-Processing of Time-Series Landsat Images
We used 174 cloud-free Landsat TM and Landsat ETM+ images of the study area (scene = path 202, row 34) taken between November 1984 and July 2007 (see Supplementary Material Table 2) and processed by the GIS and Remote Sensing Lab of the Doñana Biological Station (LAST-EBD).The number and dates of the images for each year differed because of cloud cover or acquisition failures.This set constitutes the largest remote sensing time-series with comparable sensors for the area.
Pre-processing of Landsat data was carried out by the LAST-EBD.A total of 62 images, including those with ground-truth data, were co-registered with 80-100 control points (GCPs) to a Landsat ETM+ reference image acquired on July 18th 2002 (RMS < 1 pixel).Due to time constraints, the rest of the images were georeferenced using a short-cut approach that consisted on using nine GCPs to georeference the image (RMS < 1 pixel).The resampling method was cubic convolution.Images were radiometrically corrected and transformed into reflectance values using Pons and Solé-Sugrañes [44] method implemented in MIRAMON [45].Images were then normalized using a set of pseudo-invariant areas expected to have low seasonal changes in reflectance and covering all the range of reflectance values [see 46 or 47 for a detailed description of the image processing procedure].In images georeferenced with nine GCPs, radiometric correction and normalization was applied only to Band 5 (λ = 1.55-1.75µm) due to time constraints.We chose Band 5, an infrared spectral band, because we expected an accurate discrimination between water and land classes due to differences in absorption of radiation in the near-infrared part of the spectrum [48].

Ground-Truth Data: In Situ Delineation of Flooded Area
During the two periods 2005-2006 and 2006-2007, we systematically surveyed different portions of the study area in six different dates, coincident with an overpass of Landsat 5 or Landsat 7 satellite (day of overpass ± one day).On each survey date, we visited all the temporary ponds in the portion of the area surveyed (see Supplementary Material Table 3 and Supplementary Material Figure 1).Surveyed ponds were located in the LOW SANDS, NORTHERN SANDS or SOUTHERN SANDS ecosections.The water presence was recorded and the limit of surface water was registered with a differential Leica GS20 GPS receiver.GPS data were post-processed to achieve sub-meter locational accuracy.Field visits were guided with the map of ponds with a resolution of 5 m (see Figure 1).

Model Building
Using the limits of temporary ponds, as delineated in situ, and the Landsat image of the closest date, we fitted a Generalized Linear Model (GLM) [49] using as response variable the fraction of a 30 m pixel covered by water, with binomial errors and a logit link.The response variable was in the form of a bound vector (number of 3 m flooded pixels/number of 3 m dry pixels, within each 30 m pixel, as extracted from field data projected on a 3 m grid, see Figure 2).As field-assessed completely dry pixels outnumbered flooded pixels and this could be a severe bias in model building, we balanced the data in order to have the same number of flooded and dry pixels in each sampling date.We balanced the data by means of randomly selecting a sample of dry pixels equal in number to those partially or totally covered by water (balanced data set).We tested as potential predictors the pixel normalized reflectance in each spectral band (B) as well as the value obtained when correcting for its reflectance in the absence of water [Corrected band: CB = (B − B DRY ) / B DRY ] (Table 1).Soil reflectance in the absence of water (B DRY ) was extracted from a Landsat image taken in August 2006, when all temporary ponds were dry.Predictor selection was based on percentage of deviance reduction using a manual stepwise procedure.

Model Validation
Model validation was assessed by calculating the Spearman correlation coefficient between the observed fraction of water cover in a given pixel and the ratio predicted by the model.We assumed that the predicted probability of flooding is a surrogate of the fraction of the area covered by water.We used a cross-validation procedure to test the model.In each turn, we split the data set into a construction set (with data from five of the six sampling dates) and a validation set (data from the date that was left-out from model construction).We used the balanced data set for the construction and all the data available (unbalanced data set) for the validation set.The process was repeated until all dates were used as validation sets.We conducted four different model validations to evaluate (i) if the GLM model should be applied to all pixels in the area or (ii) applied only to pond pixels (those that could be potentially flooded).And additionally, (iii) if the large tail of pixels with a low probability value of being flooded should be considered as flooded pixels, with a very small fraction of water cover, or alternatively (iv) these pixels should be considered as dry pixels using a threshold probability value.The rationale behind is that zero values (completely dry pixels) are not easily obtained from such GLM models and, thereby, we could expect low values of flooded area to be predicted for field-assessed dry pixels.In that case, the original lowest values of model predictions should be recoded to zero values.To compute the zero-codification, we established a value as threshold that misclassified as partially flooded less than 5% of ground-truth-assessed dry pixels.On the other hand, we considered a pixel as a potentially flooded pixel if it was included as water in the 5 m resolution pond map layer, or in a buffer area of 30-m from the edge of the temporary pond.The radius of the buffer area equals Landsat pixel size to overcome georeferencing biases (note that RMS < 1 pixel).
We also evaluated the temporal differences in model performance.So, we assessed whether the date of image acquisition explained the difference in the magnitude of model errors obtained during the validation phase.We computed an ANOVA analysis with posthoc Tukey test.We considered the log-transformed model absolute residuals (observed fraction-predicted fraction) as response variable and the date of image acquisition as grouping factor.
Table 1.Parameters of the final GLM model used to predict the fraction of water cover in a 30 m pixel from Landsat TM and ETM+ normalized reflectance.F-value, degrees of freedom (df) and p-value are shown.We tested as potential predictors the normalized reflectance of each spectral band (Band, B) and its value correcting for its reflectance in the absence of water (Corrected band, CB).The final GLM model, fitted to data from all ground-truthing dates, was applied to the time-series of Landsat images to obtain a map of the fraction of water in each potentially flooded pixel on each date.We modified model predictions by applying the zero-threshold value yielded by the validation procedure.Pixels in permanent ponds that contained water in August 2006 (136 pixels) were excluded from further pixel-based analyses, and analogously, those permanent ponds were removed from further pond-based analyses (i.e., average hydroperiod, pond annual hydroperiod or pond density).

Spatio-Temporal Variation in the Distribution of Water
We assessed the spatio-temporal differences in flooded area and duration of flooding (temporary pond hydroperiod).We measured flooded area as the sum of the predicted fraction of water cover for all pixels.We measured temporary pond hydroperiod as the ratio between the number of images in which it appeared to be flooded (at least one pixel of water) and the total number of images.We transformed ratio values to hydroperiod values measured in months.For each temporary pond, we computed two different values of temporary pond hydroperiod: average hydroperiod to depict its hydrological behavior over the entire study period, and computed using all remote-sensing images, and annual hydroperiod computed for each particular year.Average hydroperiods were used to assess spatial differences in flooding duration between the ecosections over the entire study period and annual hydroperiods were used to evaluate inter-annual trends in temporary pond flooding duration.We only computed temporary pond annual hydroperiods in years with seven or more images and at least one image in each season (autumn, winter, spring, summer), thus resulting in hydroperiod data for 12 years.

Differences in Hydrologic Behavior among Ecosections
We assessed the flooded area, temporary pond density and temporary pond size in each ecosection at the time of the largest flood event predicted for the study area (February 15th, 1990).We computed an ANOVA with post-hoc Tukey test to test for differences in average hydroperiod among ecosections.We square-root transformed average hydroperiod values in order to achieve normality in model residuals, a critical ANOVA assumption.

Seasonal Hydrologic Behavior of Temporary Ponds
To summarize seasonal variation in flooded surface, we plotted the estimated flooded area for each image in an annual scatter graph.The date of acquisition (number of days from September 1st) was plotted in the X-axis and the total flooded surface in the Y-axis.

Inter-Annual Variation in Hydrologic Behavior
We conducted a trend analysis of annual rainfall, maximum flooded area and pond annual hydroperiod during the entire study period.Analyses were computed both for the entire study area and for each ecosection.In the case of pond annual hydroperiod, we also computed a trend analysis for each temporary pond in order to build a map representing such tendency.

Relationship between Flooded Area and Rainfall
We built a multiple linear regression model to assess the effect of rainfall timing in the predicted flooded area for each Landsat-acquisition date.We evaluated six rainfall-derived predictors: rainfall in the previous (i) 1 to 15 days, (ii) 16 to 30 days; (iii) 31 to 90 days; (iv) 91 to 180 days; (v) 181 to 365 days and (vi) 366 to 760 days.We square-root transformed the dependent variable (flooded area) to achieve normality in model residuals, a critical ANOVA assumption.We conducted an automatic forward stepwise procedure of variable selection.We searched for the standardized regression coefficients in order to compare the relative contribution of each independent variable in the prediction of flooded area.
We built eight different multiple linear regression models, one for each ecosection and one for the entire study area.In the model conducted for the entire study area, we searched for temporal differences in model performance with an ANOVA of the regression residuals with posthoc Tukey test.We conducted two ANOVAs, one testing the year in which the image was taken as grouping factor and the other testing the month.

Model Development and Validation
We monitored 64 different temporary ponds during the ground-truthing campaigns in 2005-2006 and 2006-2007.Since each temporary pond was visited on different sampling dates, we delineated a total of 164 temporary pond perimeters.
The best model consisted in a single predictor (CB5) and significantly explained 47.09% of the deviance (Table 1).The inclusion of the second best predictor (CB4) contributed to a 9.1% increase in explained deviance.However, we did not include CB4 in the final model because its inclusion would yield a reduction of 64% in the available satellite images due to the fact that radiometric correction of band 4 was not available for 112 out of the 174 images due to time constraints.Model validation showed that the best model predictions were obtained when only potentially flooded pixels were considered and the zero-threshold value was applied (Table 2).Threshold value of the predicted fraction of water coverage was set to 0.21 and this implied that temporary ponds smaller than 189 m 2 could not be identified.The mean size of temporary ponds that were flooded but were classified as dry with the threshold value was 0.20 ± 0.21 [S.D] pixels, corresponding to 180 ± 189 m 2 .Applying this threshold both to the observed and predicted values of pixel flooding, we correctly classified 72% of the temporary ponds visited during the ground-truth campaign (Table 3).The analysis of model residuals showed that there were differences in model adjustment among sampling dates (ANOVA Table 2. Model validation with Spearman correlation between the observed fraction of water cover in a given pixel and the ratio predicted by the model.The number of pixels used for the correlation (N), the test statistic (Spearman R) and p-value are shown.The model was validated under four different circumstances: (1) Model predictions of all pixels; (2) Only pixels in potentially flooded areas could be flooded and those outside were assumed dry; (3) All pixels were considered but model predictions below 0.21 were considered as dry pixels and (4) Only pixels in potentially flooded areas with predictions equal and above 0.21 were flooded.

Historical Reconstruction of Water Coverage in the Study Area
Using Landsat imagery, we detected 864 out of 883 temporary ponds flooded at least once over the entire study period (see Supplementary Material Figure 2).The flooded area in any one date ranged between 0.81 ha, equivalent to nine pixels, (October 2005) and 246 ha (February 1990).At the time of a maximum flood (February 1990), we detected 718 flooded temporary ponds.Mean values of annual hydroperiod ranged from 3.7 ± 2.9 [S.D.] months in 1994-1995 to 6.3 ± 3.7 [S.D.] months in 1988-1989.Table 3. Confusion matrix showing the correspondence between ground-truth-assessed temporary ponds and the flooding status predicted with the GLM models after applying the zero-threshold value both for observed and predicted values of pixel flooding (flooded pixel: fraction of water cover ≥ 0.21).The seven ecosections showed different flooded areas and temporary pond densities (Figure 3), although similar distribution of temporary pond sizes (Figure 4) at the time of the largest flood event predicted for the study area (February 1990).Ecosections also differed in the average hydroperiod over the entire study period (ANOVA F 6,916 = 22.918, p < 0.001) (Figure 3).The highest average hydroperiod was found in the ECOTONE MARSHES-SANDS ecosection, where water lasted twice as long as in other ecosections during the entire study period.Among the others, the average hydroperiod in the LOW SANDS and SOUTHERN SANDS was significantly larger than in the rest of the ecosections (Figures 3 and 5).Flooded area presented a large intra-annual variability, mainly in years with annual rainfall higher than the long-term average (Figure 6).As a rule, flooded surface was low during the dry season (June-August) and continued low in September-November until the first autumn-winter rains (Figure 6).Values of flooded surface at the beginning of autumn were similar between rainy and dry years; however, flooded surface in summer months was notably higher in rainy years.Flooded surface tended to increase from September to December and tended to decrease from April-May up to August.We observed a notably increase in flooded surface, especially in rainy years, at the end of autumn months (December).Similarly, we observed a marked decrease in flooded surface between April-May and June-July in rainy years while such decrease was more gradual in dry years.

Inter-Annual Variation in Hydrologic Behavior
Maximum flooded surface and annual pond hydroperiod fluctuated on an inter-annual basis.We observed a significant temporal trend to shorter annual hydroperiod values for the entire study area (y = 6.00 − 0.10t; F 1,10 = 11.632;p = 0.007) and for all ecosections except in the ECOTONE MARSHES-SANDS and in the LOW SANDS (Table 4).We did not observe trends in annual rainfall nor in maximum flooded surface for the entire study area (Figure 7).In contrast, we did observe a negative trend in maximum flooded surface in the drier ecosections: DUNES and DRY SANDS (Table 4).

Relationship between Flooded Area and Rainfall
Focusing on the entire study area, rainfall-derived variables explained the 80% of the variance in the percentage of flooded area (Table 5).The predicted flooded area in a given sampling date was significantly related to all the periods of accumulated rainfall we considered.In particular, rainfall fallen in the previous 31-180 days had the greatest contribution to model, evidencing that ground water has an important influence in the temporary ponds hydrologic regime.Model adjustment varied widely among ecosection.Those of the driest ecosections, DRY SANDS and DUNES, showed the lowest values of explained variance.Such models did not include all the periods of accumulated rainfall.The relevance of groundwater was also evident in all ecosection models since rainfall in the previous 31-90 days showed the largest contribution to the model along with rainfall in the previous 91-180 days in the case of the SOUTHERN SANDS ecosection.
The performance of the rainfall-flood model for the entire study area differed among years (ANOVA F 22,151 = 2.709; p < 0.001) and among months (ANOVA F 11,162 = 2.946; p = 0.001).The rainfall-flood model predicted a smaller extension of flooded area than the one observed with Landsat data in December, which usually corresponds to the month of maximum flooding.The same inconsistency was observed for the most recent years (2001-2007), mainly in 2003-2004, a very rainy year.On the contrary, Landsat data indicated a smaller extension of flooded area than the rainfall-flood model predicted for 2000-2001, a year following a dry period (1998)(1999)(2000) in which the rainfall was scarce in autumn and spring, being concentrated in the months of December and January.

Application of Remote Sensing for the Monitoring of Temporary Ponds
Landsat imagery proved to be useful for reconstructing the retrospective spatio-temporal dynamics of a network of Mediterranean temporary ponds along a 23-year period.With this basis, this study contributes to the understanding of the hydrology of these temporary ponds at a local scale.As a main advantage, the proposed methodology easily enabled mapping the status of small water bodies, many of them at sub-pixel size, with medium spatial resolution imagery (i.e., Landsat) and accurate knowledge on the spatial location of ponds.It constitutes an alternative to spectral unmixing techniques (i.e., linear mixture modelling, fuzzy-c means clustering), which map the fractional coverage of each land cover class in a pixel (i.e., pond cover, tree cover…) based on knowledge of their pure reflectance spectra [50][51][52] or similar techniques which also require spectral records of the cover types of interest [i.e., 22].It should be noted that we could have not applied unmixing techniques since we lacked "pure" land cover classes.The "temporary pond" cover class was expected to be "mixed" itself, as a combination of water, macrophyte and pond bottom reflectance, and to vary from pond-to-pond, according to changes in macrophyte cover and pond depth.The proposed methodology also constitutes an alternative to traditional methods used for pond cartography (i.e., interpretation of aerial photographs), which are usually more time-consuming.Besides, the application of aerial photography to map wetlands has produced both successful [53] and inadequate results [54].Applying traditional photointerpreation techniques to delineate temporary ponds is usually difficult because of the date of acquisition of most orthophotos, in spring or summer, when most temporary ponds are dry.The detection of dry ponds is difficult because there is no difference in the reflectance of the pool bottom from the surrounding [53] and they usually lack an identifiable basin.
The main drawback of our methodology was the decrease in its performance at times of large flooding, when water levels overflowed temporary pond basins and extended over areas that do not frequently flood.Those areas were included in this study through the use of a 30-m buffer around temporary pond basins in order to control for potential positional errors.The decrease in model performance could be explained by the fact that the spectral response of these outside-basin pixels probably had a larger contribution of bottom reflectance than of water reflectance due to the low water depth.We should also be aware that the accuracy of model predictions is expected to get worse towards the beginning of the time-series data, due to changes in Landsat sensors, caused by their degradation over time, i.e. radiometric drift [55], and changes in the vegetation cover in the temporary pond basin and surroundings, as detected in particular ponds [36,57].

The System of Mediterranean Temporary Ponds in Doñana National Park
Mediterranean temporary ponds in Doñana Biological Reserve (6,794 ha) may extend up to 200 ha in very rainy years, thus conforming a large system of temporary ponds of different sizes.Their contribution to the wetland system is remarkable in spite of the fact that they have been frequently overlooked in many previous studies, only focused to the marshes [i.e., 47], permanent ponds [i.e., 35,58] or large temporary ponds with long hydroperiod [i.e., 59,60].This study, which is the first considering such a large period of time, reinforces the complex hydrology of the study area previously reported in more temporally restricted studies [61].For instance, a remarkable fact is that the largest flooded area was recorded in February 1990 despite the period that accumulated more rainfall took place from 1995 to 1997.This study also gives additional support to the ecosection zonation proposed by Montes et al. [41].Ecosections presented differences in flooded area, temporary pond density, average hydroperiod, pond annual hydroperiod and response to rainfall input.Such wide spatial variability evidenced that the system of temporary ponds in Doñana Biological Reserve is a widely heterogeneous system, being this characteristic an evidence of its importance for biodiversity conservation.Ecosystem heterogeneity favors the diversity of their associated species since it provides more niches and diverse ways of exploiting the environmental resources [see 62 for a review].For example, a wide diversity of temporary pond hydroperiods has been reported to benefit macroinvertebrate and amphibian communities at a landscape scale [1,[63][64][65].
The conservation importance of this temporary pond system is enhanced by its temporal variability, since it increases habitat heterogeneity over time.Temporary ponds varied from year-to-year in their flooded area and hydroperiod, both habitat characteristics critical for pond-breeding species, such as amphibians or macroinvertebrates [1,63,[66][67][68].The ecological relevance of such environmental fluctuation is that it provides opportunities for temporal niche partitioning: habitat conditions of a pond will favor different species in different years, depending on their niche requirements [69].As an example, inter-annual differences in the timing of pond filling and desiccation may be critical for pond-breeding species, depending on their phenology.In the case of amphibians, early-breeding species [i.e., Pelobates cultripes or Pelodytes ibericus, following 70] will be favored in years with abundant autumn rainfall whereas years with late pond filling and early desiccation will likely favor the reproduction success of species breeding in ephemeral ponds [i.e.Bufo calamita or Discoglossus galganoi, after 10].
It is of special relevance that we have detected a generalized inter-annual tendency to shorter annual hydroperiods, but not to lower annual flooded area or rainfall input.This fact suggests that annual hydroperiod shortening may have a cause independent of the natural flooding regime of temporary ponds and it is probably related with groundwater dynamics.Notably, the ecotone was not affected, probably because ponds here receive deep aquifer water discharges [71] and hence are less sensitive to changes in water-table depth.A plausible driver of hydroperiod deterioration is groundwater exploitation [72,73], which started in the 1970s, peaking in the 1980s, and currently is more or less stable but with variable emplacements.The exploitation of the aquifer is causing a progressive lowering of the phreatic level [74] and a probably damage to ponds hydrology [36], such as a decrease in the frequency of appearance of temporary ponds [75].Since pumping is concentrated in the more permeable layers of the lower aquifer, it is supposed to affect a large area through a small lowering of the water-table [33].This would explain that a trend that had been already reported for some particular ponds [36] is general for the whole study area, as it was earlier predicted [76] and this study evidences.It should be noted that, along with the potential anthropogenic cause, the progressive lowering of the water-table level also has a natural origin due to slow transient evolution of the aquifer system that cannot be avoided.From an ecological perspective, a progressively wetland desiccation is critical for amphibian communities, as revealed in the amphibian populations decline reported for Yellowstone National Park [77].In particular, in the study area, the reduction in pond annual hydroperiod may severely compromise the medium term population stability of pond-breeding species with an aquatic phase requiring a long period of time to complete metamorphosis, such as P. cultripes.
The flooded area was related both to immediate and previous rainfall inputs, including rainfall input in the previous year.Similar results have been reported in a different temporary pond system, such as the playa-lakes in the Monegros Dessert [15].This result adds to evidence of the relevance of groundwater influxes in the hydrological regime of the ponds in the area [34].Temporary ponds in Doñana National Park are shallow depressional wetlands that flood when the water-table raises above the topographical surface.So, the accumulation of rainfall water is necessary both to initially recharge the aquifer and then to fill the ponds.The amount of water required to completely recharge the aquifer depends on the aquifer level at the end of the previous year and thereby on the annual rainfall during the previous year.Serrano and Zunzunegui [36] provided an analogous result, showing that the rainfall input in a given year was related to pond hydroperiod in the following one.
The conservation value of Doñana National Park for pond-breeding species lies on its heterogeneity across space and time [43], since communities in wetland ecosystems require from spatial variability in order to be able to persist under high environmental variability [78].Areas with such a high density and heterogeneity of natural temporary ponds are not common in Europe, where the number of temporary ponds are probably a mere fraction of what they would naturally had been in the past [79].For all these reasons, we think that the system of temporary ponds in Doñana deserves special attention and conservation measures focused to avoid its destruction or degradation.In particular, we should be concerned for the tendency to lower hydroperiods which may severely compromise their suitability as habitat for pond-breeding species.

Figure 2 .
Figure 2. Scheme of the model development and application procedure.

Figure 3 .
Figure 3. (a) Flooded area and temporary pond density (number of temporary ponds/km 2 ) predicted for each ecosection at the time of the largest flood event in the study area (February 15th, 1990), as predicted from Landsat data.The percentage of flooded area is also shown.(b) Mean and standard error of average hydroperiod for each ecosection.

Figure 4 .
Figure 4. Distribution pattern of temporary ponds size in each ecosection at the largest flood event in the study area (February 15th, 1990).

Figure 5 .
Figure 5. (a) Cartography of average hydroperiod (frequency of water occurrence) over the entire study period.(b) Cartography of the trend in annual hydroperiod for each temporary pond over the entire study period.Hydroperiod trend is estimated as the slope of the regression between annual hydroperiod and year (months/year).

Figure 6 .
Figure 6.Synoptic characterization of intra-annual variation of flooded area as obtained from pooling pond model predictions from the entire time-series data (1984-2007).

Figure 7 .
Figure 7. Inter-annual trends in (a) annual rainfall, (b) maximum flooded area and (c) annual hydroperiod (mean and standard error) for the entire study area.Linear fit and adjusted equation is also shown.Years in grey indicate that no representative or enough images were available to compute annual hydroperiod that year.

Table 4 .
Trend analysis of flooded area and annual hydroperiod in each ecosection.Model parameterization, where t corresponds to "time", F statistics, degrees of freedom (df) and significance (p) are shown.

Table 5 .
Model results for predicting the flooded area from rainfall data.Models were conducted for the entire study area and for each ecosection.Predictor nomenclatures correspond to the number of days the accumulated rainfall was computed (i.e., P 1-15 = Rainfall collected during the days 1 to 15 before the image date).