Late Holocene evolution of playa lakes in the central Ebro depression based on geophysical surveys and morpho-stratigraphic analysis of lacustrine terraces

Available online 26 February 2012


Introduction
The central sector of the Ebro Depression in NE Spain (Fig. 1) constitutes a particularly interesting and challenging area for palaeoclimatic investigations due to several reasons: (1) It is the northernmost semiarid region in Europe, bounded by mountain ranges that were glaciated during the Quaternary. (2) Climate variability during the last glacial cycle may differ substantially from other European areas (Harrison et al., 1996;Prentice et al., 1998;Valero-Garcés et al., 2000a;Davis et al., 2003;González-Sampériz et al., 2008) due to the combined influence of both North Atlantic and sub-tropical climates (Summer et al., 2001). (3) Lacustrine records are largely restricted to saline lakes of poorly understood geomorphic origin and evolution, including the northernmost playa-lakes with active evaporite deposition in Europe.
Studies focused on saline lake sediments provide evidence for a complex palaeohydrological evolution over the last glacial cycle in the semiarid Ebro Depression, with several abrupt and rapid arid/humid transitions (Valero-Garcés et al., 1998, 2000a. However, the chronological control of most lake sequences is rather limited due to the scarcity of organic remains in the sediments and the presence of erosional hiatuses. Holocene palaeoenvironmental interpretations of the Ebro Depression are based mostly on cores drilled in saline lakes located in internally drained depressions: Mediana (Valero-Garcés et al., 2000a,b), Bujaraloz-Sástago (Pérez-Obiol and Roure, 1990;Stevenson et al., 1991;Davis, 1994;Schütt, 1998;Davis et al., 2003;Moreno et al., 2004;González-Sampériz et al., 2008;Mees et al., 2011), Chiprana (Davis, 1994;Valero-Garcés et al., 2000b), Alcañiz (Stevenson et al., 1991;Davis, 1994) and Hijar areas (Davis, 1994). In spite of the large efforts carried out to investigate their stratigraphic record from a palaeoenvironmental perspective, the origin and morphostratigraphic evolution of the enclosed basins remain unclear. The lack of data on the subsurface geometry of the basins precludes evaluating the validity of the different hypotheses proposed to explain the formation of the closed depressions in evaporitic bedrock. Additionally, reconstructing the evolution of these lacustrine basins requires taking into account the scarcely investigated terraces preserved at the lake margins.
As González-Sampériz et al. (2008) point out, stratigraphic records from the bottom of playa-lakes in the Ebro Depression pose significant limitations to palaeoenvironmental studies for several reasons: (1) limited thickness and temporal length of the sequences; (2) significant hiatuses attributable to enhanced wind erosion, difficult to identify in cores; (3) very low pollen content and material datable by the radiocarbon method; (4) significant diagenetic overprinting favoured by the continuous upward discharge of saline groundwater flows and the inherent instability of evaporitic minerals (Warren, 2006); (5) reworking and contamination processes; and (6) limited preservation potential of biological indicators. Recently, detailed geomorphological mapping has revealed the presence of stepped sequences of aggradational lacustrine terraces at the margin of some playa-lakes in the Bujaraloz-Sástago area (Gutiérrez-Elorza et al., 2002). This finding opens new prospects to palaeoenvironmental investigations in playa-lakes, since terrace deposits, in combination with cores from the lake bottoms, may help to improve the completeness and resolution of lacustrine records. Moreover, lacustrine terraces may be used to reconstruct the morpho-stratigraphic evolution of the playa-lakes and identify major aggradation and excavation periods as potential responses to significant palaeoenvironmental changes. Advantages of the terrace deposits, situated at the lake margins and above the water table, with respect to the sediments situated beneath the lake bottoms, include: (1) higher preservation potential of the stratigraphic record than in the lake bottom, where eolian deflation may remove a substantial amount of sediment; (2) less severe overprinting by diagenetic processes; and (3) they may be studied in natural exposures or excavated trenches, allowing the architecture of the deposits to be analysed and increasing the likelihood of finding datable material. Moreover, the identification of the boundary between lake deposits and weathered bedrock can be carried out with a greater level of confidence in exposures; argillaceous karstic residues are frequently misinterpreted as lacustrine deposits. Once the sequence of lake terraces has been established, the main challenge is to determine the stratigraphic and chronological relationships between the terraces and the deposits situated in the lake bottom. Some of the stratigraphic units identified in cores from the bottom of lakes may not have chronostratigraphic equivalents in terraces and vice versa. The main constraint for integrating such data is the availability of reliable numerical dates.
This paper is focused on the Bujaraloz-Sastago endorheic area, and particularly on the largest La Playa-El Pueyo lake system. The main goals of this multidisciplinary work include: (1) evaluating the different hypotheses proposed to explain the origin of the closed depressions developed in evaporitic bedrock in the light of new data on the subsurface geometry of the basin inferred by complementary shallow geophysical techniques (electrical resistivity tomography and seismic refraction); (2) reconstructing the morpho-stratigraphic evolution of the playa-lake on the basis of the spatial distribution of lacustrine terraces and the obtained numerical dates; (3) calculating rates of vertical accretion and eolian deflation using the numerical ages obtained for the most extensive terrace; and (4) making an attempt to integrate the terrace (trenches) and lake bottom (cores) stratigraphic records and relate the palaeohydrological changes recorded by the terraces with other palaeoenvironmental proxies in the Ebro Depression, the Iberian Peninsula and the global context.

Geological setting
The study area is located in the central sector of the Ebro Tertiary Basin, NE Spain (Fig. 1). This sedimentary basin, deeply dissected by the fluvial network, constitutes a topographic depression drained longitudinally by the NW-SE-oriented Ebro River and bounded by the Pyrenees and the Iberian Chain to the north and south, respectively. The investigated La Playa-El Pueyo playa-lake complex form part of a field of closed depressions with a marked WNW-ESE elongation designated as the Bujaraloz-Sástago endorheic area ( Figs. 1 and 2). This internally  Pueyo,3. El Pito,4. La Salineta,5. Guallar,6. Camarón,7. Rebollón,8. La Clota. drained area has been developed on an exhumed structural platform lying at 310-370 m a.s.l., hanging 200 m above the deeply entrenched Ebro River to the south (Fig. 1).
The sedimentary fill in this sector of the Ebro Basin is made up of Oligo-Miocene sediments deposited in evaporite and carbonate shallow lakes and in distal alluvial fan environments (Fig. 1). The strata, showing a very low (b2°) NW to NE dip, form part of the southern limb of a very open WNW-ESE syncline whose pericline is located SE of Bujaraloz (Quirantes, 1978;Arlegui, 1996). The structural platform is underlain by the Lower Miocene Bujaraloz-Sariñena Unit, mostly composed of gypsum, mudstone and limestone (Ramírez, 1997;Solá and Costa, 1997). Salvany et al. (1994Salvany et al. ( , 1996 differentiate two main gypsum-rich units, each underlain by a detrital unit primarily composed of red clays (Fig. 1). From base to top: Middle Detrital Unit (15-20 m), Middle Gypsum Unit (40 m), Upper Detrital Unit (5-6 m), Upper Gypsum Unit (100 m). Most of the closed depressions have developed on the Middle Gypsum Unit, exposed in the southern sector of the platform (Fig. 1). This unit has a much lower proportion of clay than the thicker Upper Gypsum Unit. Moreover, the size and spatial frequency of the depressions decrease towards the east, consistently with the wedging out of the Middle Gypsum Unit and the lateral change to marl and carbonate facies (Salvany et al., 1994(Salvany et al., , 1996. La Playa and El Pueyo playa-lakes occur on the Middle Gypsum Unit, whereas La Salineta has been developed on the Upper Gypsum Unit (Fig. 1). The Miocene bedrock is affected by small-throw normal faults with dominant NW-SE azimuths and subvertical joints. Arlegui (1996, p. 217) and Solá and Costa (1997) report a primary N-S trending joint set and a secondary one with E-W orientation. Gutiérrez-Elorza et al. (2002) identified NE-SW and NNW-SSE prevailing joint directions based on 319 measurements taken in limestone beds exposed in La Playa area.

Climate and vegetation
The Ebro Depression has a continental climate, characterised by very hot summers and cold dry winters. According to the records of Bujaraloz-Petris meteorological station, located 8 km north of the studied playa-lakes, the average annual precipitation and temperature are 360 mm and 14.4°C, respectively. Precipitation shows a bi-modal seasonal pattern, with the highest rainfall in spring and autumn (Rodó et al., 1997). García-Vera (1996) has estimated annual potential evapotranspiration values of 788 and 909 mm applying the Thornthwaite and Blaney-Criddle methods, respectively. The area is characterised by a strong wind with prevailing WNW direction locally designated as Cierzo. This cold and dry wind blows mainly in winter and spring, channelled along the Ebro Depression and controlled by the presence of anticyclonic conditions in the Cantabrian Sea and low pressure in the Mediterranean Sea (Ascaso and Cuadrat, 1981;Puicercús et al., 1997). At Zaragoza meteorological station, located 60 km to the NW, the maximum wind speed recorded over the period 1942-2010 reached 37.5 m/s (135 km/h) in February 1954 with a WNW direction. The data recorded in Bujaraloz anemometric station at 10 m above the ground between March 1991 and January 1993 indicate that the highest mean velocities (6-6.5 m/s) correspond to WNW and NW winds, which are the most frequent ones and represent about 75% of the eolian energy. During this period a maximum WNW wind speed of 27.7 m/s (100 km/h) was measured and the wind velocity exceeded 4 m/s around 45% of the time (Puicercús et al., 1997).
Vegetation in the central Ebro Basin is currently an herbaceousshrubby land, mostly dedicated to winter cereal farming. Broad leaf tree formations are restricted to humid canyons in elevated areas (Molero and Blanché i Vergés, 1986;Longares, 1997). Playa-lake basins typically host vegetation dominated by halophytic plant communities, whose distribution reflects water availability and a salinity tolerance gradient (Aguilella and Riera, 1997).

Hydrogeology
The Upper and Middle Gypsum Units constitute two aquifers separated by the Upper Detrital Unit, which behaves as a leaky aquitard (Fig. 1). These gypsiferous hydrostratigraphic units have a rather low average permeability of around 0.01 m/day, largely controlled by solutionally-enlarged joints. Hydraulic permeability reaches the highest values in the bottom of the depressions, due to enhanced karstification by groundwater flow discharge (García-Vera, 1996;Salvany et al., 1996;Samper-Calvete and García-Vera, 1998). The groundwater flow in this endorheic area is controlled by the topography of the platform; an extensive plateau riddled by enclosed depressions. Water recharge in the aquifers, estimated at 20-45 mm/yr, flows towards and discharges at the bottom of the main basins, forming local and centripetal groundwater flow cells (García-Vera, 1996;Salvany et al., 1996;Samper-Calvete and García-Vera, 1998;Sánchez-Navarro et al., 1998). The groundwater salinity increases progressively along the flow path, changing from a Ca-SO 4 composition in the recharge areas, into a Na-(Mg)-Cl-SO 4 hydrochemical facies in the discharge zones (García-Vera, 1996;Salvany et al., 1996). Brines in the playa-lakes leads to precipitation of salts at the surface and in the subsurface (Pueyo, 1978(Pueyo, /79, 1980Pueyo and Inglés, 1987). Upward flow and capillary rise in the bottom of the playa-lakes, whose extremely flat topography is controlled by the water table zone, are prompted by the high evaporation rate. The ongoing irrigation plan (Monegros II), affecting 200 km 2 of the endorheic area, will have a significant impact on the hydrology and hydrochemistry (water level rise, dilution, pollution) of these highly sensitive systems with unique habitats and endemic species (García-Vera, 1996).

Geomorphology of the playa-lake basins
A total of 99 closed depressions were inventoried by Balsa et al. (1991) in the Bujaraloz-Sástago endorheic area, which covers approximately 250 km 2 . According to these authors, the bottom of these topographic basins has a total area of around 8.8 km 2 . The majority of the depressions are markedly elongated and oriented in the WNW-ESE direction, which is also the prevailing wind trend. La Playa is the largest lake with 3.5 km in axial length and covering 1.72 km 2 , approximately 20% of the cumulative area of the bottom of the depressions (Fig. 2). The structural platform is also carved by poorly hierarchised flat-bottom infilled valleys with dominant WNW-ESE orientation, most of which flow into closed depressions. A significant proportion of the depressions show scarped edges, generally controlled by a laterally extensive limestone bed situated at the top of the Middle Gypsum Unit (Quirantes, 1965). Approximately 20 depressions host ephemeral saline wetlands that get flooded every year (Balsa et al., 1991). These playa-lakes have flat and moist bottoms, indicative of a topography controlled by the water table and the capillary fringe (Rosen, 1994). The largest playas (i.e. La Playa, La Salineta) show a slightly more elevated bottom in the southeastern sector, suggesting higher aggradation in the downwind sector by eolian supply and wave action and/or greater erosion in the opposite sector. Water depth reaches the highest values in winter and rarely exceeds 50 cm (Castañeda, 2002;Castañeda et al., 2005). Field and satellite data obtained from 1985 to 2000 indicate that La Playa is one of the playa-lakes with a longer wet cycle, remaining partially flooded more than six months per year (Castañeda, 2002). These elongated lakes tend to have asymmetric geometry in plan view, with higher width in the downwind half than in the windward one, which may have a pointed margin. The leeward margin of La Playa shows an embayment that seems to represent the southeastward expansion of the lake by eolian erosion of the lacustrine terrace situated between this lake and El Pueyo (Fig. 2). Goudie and Wells (1995) indicate that deflation basins excavated in lake deposits within palaeolacustrine basins typically display two types of morphologies: (1) ovoid, almost subcircular, with the long axis perpendicular to the formative air flow; and (2) cusp-shaped, with the apex pointing upwind and oriented parallel to the formative wind. A significant proportion of the playas in the studied area show similar characteristics to the latter morphological type.
Progressive evaporation of the brines in the playas leads to precipitation of salts (Pueyo, 1978/79;López et al., 1999). Algal mats 3-4 mm thick develop in the playa-lake floor during flooding periods. Saline efflorescences (bloedite, halite and thenardite) precipitate on the playalake floor when the water level recedes. Beneath the surface there is a black sapropelic mud with a high content in sulphides resulting from the reduction of sulphates by bacterial activity. Desiccation cracks develop during dry periods and evolve into saline polygons, eventually framed by teepee-like features as well as bent and overthrusted edges (Fig. 3a). Saline crusts and algal mats locally show blisters formed by the expelling of gases derived from decomposition of organic matter under anoxic conditions and by volume increase related to salt crystallisation (Pueyo, 1978/79;Fig. 3b). A very small scarp-edged pit, approximately 30 cm across and 10 cm deep, was identified in the bottom of La Playa Lake (Fig. 3c). This feature, very rare in the area, may correspond to a cover collapse sinkhole related to the presence of a cavity of limited size in the gypsiferous bedrock.
The margins of the lakes typically have an aureole of halofilous vegetation, and wind-transported particles accumulate in the leeward side of these plants forming elongated nebkha dunes (Fig. 3d). These dunes are frequently made up of sand-sized lenticular gypsum crystals. According to Pueyo (1978/79), these intraclasts result from eolian redeposition of crystals formed in flooded areas and subsequently dried up. Trains of nebkhas may merge resulting in linear dunes several tens of metres long (Fig. 3e). Several processes favour wind deflation in the bottom of the playa-lakes during desiccation periods: (1) Drying cycles involve a significant reduction in cohesion of fine particles, increasing their susceptibility to wind entrainment.
(2) Precipitation of salts at the surface may involve the accumulation of light and loose crystals and the preparation of particles susceptible to deflation by disintegration and salt weathering (i.e. edge of desiccation polygons and collapsed crust blisters; Fig. 3a, b). (3) Trampling by animals (Thomas, 1988;Goudie and Wells, 1995); and (4) accumulation of faecal pellets by worms. Field evidence reveals that this type of biological burrowing activity may predispose a substantial amount of material to wind erosion from the subsurface (Fig. 3f).
Detailed geomorphological mapping was carried out by Gutiérrez-Elorza et al. (2002) in a sector of the endorheic area using a 1:5000 scale topographic map with a contour interval of 1 m. Three lacustrine terrace levels were identified in La Playa and El Pueyo lakes (Fig. 2). This stepped sequence of terraces records alternating periods during which the lakes were dominated by aggradation or deepening. In La Playa-El Pueyo basin, the upper terrace is represented by a few small benches situated 9 m above the playa-lake bottom in the southeastern margin of El Pueyo Lake. The middle terrace, around 6 m above the lake bottom, forms extensive benches in the southeast leeward margin of the playas. The surface of this terrace locally merges with the top of the deposits filling flat-bottomed valleys that drain into the playalakes (Fig. 2). The spatial distribution of this terrace reveals that during its accumulation there was a larger single lake around 2.7 km 2 in the depression where La Playa and El Pueyo are located. Subsequently, differential eolian erosion compartmentalised the lake basin into La Playa and El Pueyo, nowadays separated by a remnant of the middle terrace. The lower terrace, with a limited extent, is situated at about 0.5 m above the lake bottom. This apparently very recent level forms discontinuous benches at the margins of the playa-lakes, mainly in the windward sector. It is not clear whether the current lake bottom is a strath surface cut-across the deposit of the lowest terrace or the aggradation surface of an inset younger stratigraphic unit.
The modern La Salineta Lake, 0.2 km 2 in area, lies within a much larger palaeolake, whose deposits have been partially eroded and form scarps up to 4 m high surrounding the present lake. Two main terraces have been recognised. The upper terrace is best developed at the southeastern end of the basin. Remains of another small terrace at a relative height of 1 m occur at several locations.
A total of 50 yardangs up to 264 m long with a dominant WNW-ESE trend have been mapped in the leeward side on the largest playas, mainly La Playa (1.72 km 2 ), El Pueyo (0.14 km 2 ) and El Pito (0.35 km 2 ) lakes (Gutiérrez-Elorza et al., 2002;Figs. 2 and 4). Yardangs are elongated hills produced by wind erosion in combination with other process such as weathering and gullying. To our knowledge, these are the only yardangs reported in Europe (Goudie, 2007). The spatial distribution of the yardangs is attributed to the higher abrasive capability of the wind in the leeward sector of the larger playas, due to enhanced concentration of particles in the air currents by deflation. Most of these landforms have been developed in gypsiferous bedrock, and a few of them have been carved on unconsolidated terrace deposits (Fig. 4). These landforms provide evidence of the geomorphic effectiveness of wind erosion in the area. Additionally, the yardangs formed in the middle lake terrace deposits strongly suggest that wind erosion has played a significant role in the geomorphic evolution of the playa-yardang systems in recent times.
Several hypotheses have been proposed to explain the origin of the closed depressions in the Bujaraloz-Sástago endorheic area. According to one interpretation, the basins result from differential dissolutional lowering of the ground surface controlled by fractures (Mingarro et al., 1981). These would be solution sinkholes generated by downward vadose flows in the epikarst zone (Williams, 2004;Beck, 2005; Gutiérrez et al., 2008a;Gutiérrez and Cooper, in press). This model is not compatible with the available hydrogeological data, indicating that the playa-lake basins essentially behave as discharge zones. Several researchers explain the depressions as the result of the brittle collapse of cavities generated by structurally controlled subsurface dissolution of the gypsum bedrock (Aramburu, 1904;Quirantes, 1965;IRYDA, 1989). This explanation, which involves the development of cavities beneath the platform, fits with the bedrock collapse sinkholes of several genetic classifications (Williams, 2004;Beck, 2005;Gutiérrez et al., 2008a). However, no significant karstic voids have been found in the boreholes drilled in the area (M.A. García-Vera pers. comm.). The sole small collapse feature identified in La Playa is most likely related to the presence of rare karstic conduits of limited size (Fig. 3c). Another hypothesis involves a mixed karstic and eolian origin (Sánchez-Navarro et al., 1998;Desir et al., 2011). Initially, the depressions may form and evolve as solution sinkholes by percolating water in the vadose zone. Infiltration is favoured by the horizontal topography of the platform (Ibáñez, 1975). Once the bottom of the dolines reaches the water table zone, local groundwater flow cells that discharge in the depressions are established. The underground flows that converge in the basins cause widespread dissolution of the bedrock, leading to gradual subsidence. During certain periods, the strong WNW prevailing winds may cause the erosional lowering of the floor of the basins by deflation. According to this interpretation, the genesis of the depressions involves processes characteristic of solution sinkholes, deflation basins and subsidence sinkholes (Goudie and Wells, 1995;Williams, 2004;Gutiérrez et al., 2008a). From the hydrogeological point of view, this evolutionary model implies the transformation of a recharge basin into a hydrologically closed discharge playa (Rosen, 1994;Yechieli and Wood, 2002).
The origin and hydrogeology of the depressions is particularly important from the stratigraphic and palaeoenvironmental record perspective (Rosen, 1994). The preservation potential of the lake sediments in the bottom of the depressions would be higher if bedrock collapse is the main genetic process. This subsidence mechanism would lead to a lacustrine fill with highly variable thickness and affected by gravitational deformation. Several geophysical surveys have been acquired in La Playa depression aimed at characterising the geometry of the lake deposit and testing the different genetic hypotheses.

The previously studied lacustrine sequences
A number of sequences from playa-lakes in the Bujaraloz-Sástago area have been investigated (Pérez-Obiol and Roure, 1990;Stevenson et al., 1991;Davis, 1994;Schütt, 1998;Davis et al., 2003;Moreno et al., 2004;Valero-Garcés et al., 2004;González-Sampériz et al., 2008;Mees et al., 2011). The last review by González-Sampériz et al. (2008) summarised the vegetation changes and moisture fluctuations since the Last Glacial for the central sector of the Ebro Depression and included a detailed description of two cores retrieved from La Playa and La Salineta playa-lakes. Available chronologies only allow a millennial scale reconstruction. According to that review, the early Holocene was the wettest period, and the common lack of most Holocene sediments in the lake basins would have been caused by an intensification of eolian deflation under arid conditions during the Mid Holocene. Moreover, the records indicate relatively more humid conditions during the late Holocene after 4 ka, with more frequent flooded conditions in the lake basins.

La Salineta
An outcrop of a lacustrine terrace located in the southeastern margin of La Salineta (0-360 cm) and a core (360-465 cm) collected in the lake bottom, at the foot of the terrace riser, were described and sampled by Davis (1994). Chenopodiaceae seeds from the upper 20 cm of the 465 cm long La Salineta section gave a modern AMS age probably caused by soil contamination, and consequently the chronology of this section is unknown.
An 8 m long core (La Salineta Core) was drilled in the upper terrace close to the previous section, providing one of the longest lacustrine records available in the central Ebro Depression (Fig. 5). The core reached the bedrock composed of Miocene limestone . No terrestrial macro organic remains were found in the core, so the chronology is constrained by bulk organic matter AMS radiocarbon ages. The section reaches the Late Glacial Maximum (LGM, around 20 ka cal. BP) and contains the Late Glacial and early Holocene material, although the chronology of the upper half of the sequence is not constrained. AMS samples from the upper 3 m gave modern ages, although pollen indicators suggest that the sequence contains recent Holocene (last millennium) sediments.
The modern Salineta playa-lake bed lies about 4 m below this terrace. An 87 cm long sediment core was recovered by González-Sampériz et al. (2008) in the centre of the lake. These authors described three depositional sequences (S-III, S-II and S-I) bounded by erosional hiatuses. Chronological control is restricted to a radiocarbon age from a pollen concentrate obtained in sequence S-II at 65-67 cm below the surface (2131-1966 cal yr BP). The presence of dolomite, the more positive δ 18 O values pointing to intense evaporative processes, the maximum values of clays and quartz and the dominance of halophytic terrestrial plants, as suggested by δ 13 C org values (−26‰), indicate relatively carbonate-dominated brine and an ephemeral lake system with frequent desiccation periods during deposition of undated S-III (prior to 2000 cal yr BP). Sequence S-II would represent deposition in an ephemeral playa-lake dominated by gypsum-rich sediment formation pointing to a different brine composition after 2000 cal yr BP. Finally, sequence S-I represents deposition during the present-day ephemeral saline lake, characterised by higher content of the more soluble sulphates and the dark-grey, banded to black laminated nature of the sediments. Although the chronological control is poor, sedimentological and palynological data (i.e., arboreal component less than 40% and continuous presence of Cerealia) indicate that these sediments are of recent Holocene age.

La Playa
In the 162 cm long core from the western margin of La Playa, González-Sampériz et al. (2008) differentiated three depositional sequences (P-III, P-II and P-I) bounded by erosional surfaces (Fig. 5). This succession and the pollen record, particularly the ratio between halophytes (Chenopodiaceae) and steppe taxa (Artemisia) versus aquatic plants (Ruppia, Myriophyllum, Potamogeton), record a depositional history characterised by the evolution from a carbonateproducing lake (Units P5, P4 and P3) towards a more sulphateproducing saline lake (Unit P2) some time before 9.7 ka, ending with the present-day ephemeral saline lake system (Unit P1). The upper Unit P1 is characterised by deposition of calcitic, organic-rich mud with abundant gypsum microcrystals and without the halite crust characteristic of other playa-lakes in the area. The main limitation of this record is the poor chronological control; only one age is available from the middle part of depositional sequence PII (80 cm below the surface, 9914-9676 cal yr BP; Fig. 5). The chronology of the depositional sequences and the erosional hiatuses at the boundary between them is unknown. The available data suggest that early Holocene deposits were preserved at the base of the sequence, but most of the Holocene sediments have been eroded and the age interval of the middle sequence is unknown. The top sediments (Unit P1) represent historic and recent deposition taking into account the pollen content (Olea and Cerealia cultivars). The mineralogical study by Mees et al. (2011) based on profile pits and auger holes from 14 basins in this endorheic area shows a similar recent stratigraphy.

Methods
Three trenches were excavated to investigate the deposits of the upper and middle terraces of La Playa-El Pueyo lake system at two sites (Figs. 2,4,7 and 8). The trenching technique is mainly applied to assess the seismogenic potential of active faults in palaeoseismological studies (McCalpin, 2009, in press). However, its practicality has also been proved in the study of sinkholes (Gutiérrez et al., 2009), landslides (Gutiérrez et al., 2008b and various types of surficial deposits; fluvial (Brakenridge, 1985), eolian (Mahan et al., 2009), etc. In this study, after cleaning the vertical trench walls, a reference grid with horizontal and vertical strings spaced 0.5 or 1 m apart was placed on one of the flanks of each trench (Fig. 7b). The stratigraphy of the selected walls was logged on graph paper with the aid of the grid and a tape measure (Fig. 8). Colour pins were nailed at stratigraphic contacts to facilitate the mapping and sampling process.
Sedimentary facies were described in the trench according to sedimentary textures and structures, grain size and sediment composition. Twenty six samples were collected in the 5 m thick sequence of the middle terrace exposed in one of the trenches for clay fraction and elemental composition (Total Organic Carbon, Total Inorganic Carbon, and Total Sulphur). Grain size analyses were performed with a Coulter Grain Analyzer after removing organic matter with oxygen peroxide and a treatment with ultrasounds and clay dispersants. TOC, TIC and TS analyses were performed with a LECO elemental analyzer. Pollen samples were taken in the upper (2 samples) and middle (18 samples) terrace of La Playa-El Pueyo lake system. The chemical treatment of these samples has followed the usual method that includes HF, HCl, KOH and Thoulet density 2.0 heavy liquid (Dupré, 1992). However, the pollen preservation was very low, and most samples were sterile or with very poor taxonomic diversity. Thus, a vegetation or chronological comparison between the lacustrine core  and the trench sequences was not possible.
Seismic refraction and electrical resistivity tomography (ERT) profiles were acquired in La Playa Lake. Seismic refraction utilises seismic compression waves refracted at the interface between units with different propagation velocities and recorded by a seismograph and a linear array of geophones at the surface. The analysis of the first-arrival times of direct and critically refracted waves allows modelling the depth profile of refractor contacts and the propagation velocity (v p ) of each layer. Seismic refraction was selected as an adequate method to investigate the depth and geometry of the boundary between the nonconsolidated lake deposits and the gypsiferous bedrock (i.e. Hetch, 2003;Hoffmann and Schrott, 2003;Schrott and Sass, 2008). Lines of P-wave refraction data were acquired using a 12-channel Geometrics Strataview seismograph, 10 Hz geophones and a sledgehammer energy source. Spacing in the geophone arrays was 5 m, except for the first and last geophones located at 2.5 m from the first and last shot points, respectively. Three shot points were used with each spread, enough to obtain information on the geometry and properties of shallow layers. The span of the traverses was constraint to 55 m both by the available cable length and the ability to generate sufficient energy with the sledgehammer through unconsolidated lake deposits. The recorded seismic data was processed using the SIPT2 iterative inversion refraction analysis software, which generates the first solution by a delay-time method, while subsequent iterations modify the solution through ray tracing (Scott, 1973).
Electrical resistivity tomography (ERT) allows imaging spatial variations in the electrical resistivity of subsurface materials by processing electrical measurements from electrode arrays (i.e. Griffiths and Barker, 1993;Sumanovac, 2006). Bulk electrical resistivity (ρ), the resistance of the ground to the passage of an electric current, is determined by the amount of conducting minerals and specially the chemical composition of interstitial water. It may be used to identify the boundary between different units, like the bedrock-cover interface (Samouëlian et al., 2005;Schrott and Sass, 2008;Van Dam, 2010). It may also image resistivity variations related to changes in the amount and chemistry of impregnating waters within the same geological unit (i.e. Bauer et al., 2006;Cassiani et al., 2006;Leroux and Dahlin, 2006;Poulsen et al., 2010). For example, the electrical resistivity of meteoric water derived from precipitation ranges between 30 and 1000 Ω m, whereas that of sea water is generally about 0.2 Ω m (Bell, 2007). ERT was a priori an adequate method to delineate the rockhead and to infer electrical changes related to permeability and hydrochemical variations. To our knowledge, there are no publications dealing with the application of this method in playa-lakes. The ERT profiles were acquired with a multi-electrode Lund Imaging System (ABEM). The system is composed of a Terrameter SAS 1000, electrode selector ES10-64 and 4 signal cable and reels and 64 steel electrodes. We selected a Wenner electrode array. Although this configuration provides lower horizontal resolution than the dipole-dipole array, it has a greater signal-to-noise ratio, yields tomographs with good vertical accuracy and may also provide a reasonable horizontal resolution, especially if the spacing between electrodes is reduced (i.e. Sasaki, 1992). The interelectrodic spacing (a) used was alternatively 1 or 2 m for the 40 central electrodes and 2a for the 24 electrodes on both ends of the line. This ensured enough resolution, especially in the central sector of the profile, without compromising the required investigation depth. The measured apparent resistivity (ρ a ) data were processed using the 2D finite-difference inversion commercial software RES2DINV of Loke and Barker (1996). The programme employs a 2D model of rectangular cells, producing through inversion of the raw data a true resistivity value for each node and an electrical resistivity tomograph.

Geophysical investigation
The seismic refraction and ERT profiles have been acquired along an E-W trending trace across La Playa, including the middle terrace at the NE margin (Fig. 6a). In order to use subsurface direct data for the interpretation of the profiles, the selected trace met the location of a borehole drilled in the lake bottom  and the trenches excavated in the middle terrace (Fig. 6a, b). A total of 4 ERT profiles (T1a, T1b, T2 and T3) and 2 seismic refraction traverses (SR1, SR3) have been acquired in three sections of the selected trace. ERT profiles T1a (80 m) and T1b (160 m), centred at the scarp of the middle terrace next to the trenches, investigate both the terrace and the lake bottom margin. Profile T1b was acquired with an electrode spacing of 2 m, whereas the spacing in profile T1a was 1 m in order to get higher resolution. The location of the 55 m long traverse SR1, performed on the terrace surface, essentially coincides with the eastern half of the electrical resistivity profiles T1a and T1b. ERT profiles T2 and T3, both 80 m long and with an electrode spacing of 1 m, were carried out in the central and western sector of the lake bottom, respectively. The core investigated by González-Sampériz et al. (2008) was obtained between profiles T2 and T3. The 55 m long SR2 seismic traverse roughly coincides with the eastern half of the ERT profile T3.
The seismic refraction profile SR1 carried out in the middle terrace reveals two main seismic units with markedly different velocities (see dashed lines in Fig. 6c). The upper unit, with homogeneous behaviour and attributable to the terrace deposit and probably also to the underlying clay-rich karstic residue, has an average velocity of 400 m/s. At the top of this unit there is a thin layer with a lower average velocity of 291 m/s that corresponds to the uppermost deposits of the terrace disturbed by ploughing. The higher thickness and channel-like base of this subunit in the eastern sector coincides with a very gentle erosional trough, suggesting the presence of higher porosity alluvial deposits and/or man-made ground. The lower unit, which corresponds to the bedrock, has a dip-corrected velocity of 3664 m/s. According to the seismic refraction profile SR1, the rockhead has an overall planar geometry and is situated at ca. 5 m in depth in the vicinity of the terrace scarp, coinciding with the position of the top of the fresh bedrock exposed in the trench dug at the terrace scarp. The ERT profiles T1a and T1b show a sharp resistivity change at approximately the same depth beneath the terrace. This geophysical contact shows a slight dip towards lake and profile T1a displays an apparent west-facing step at the reference point 8 that may be related to local changes in the karstic residue. The resistivity of the terrace deposits varies significantly with the distance to the surface, both the terrace tread and scarp. The material situated close to the ground has resistivity values of 6-25 Ω m and decreases away from the surface to 1-6 Ω m. This resistivity distribution pattern may be explained by differences in soil moisture between the surface material and the subsoil. The sections of the ERT profiles T1a and T1b situated at the lake bottom next to the terrace scarp do not allow inferring with confidence the position and geometry of the coverbedrock contact. The complex resistivity pattern recorded in this sector may be related to different factors, including: (1) significant lateral changes from lacustrine to colluvial facies associated to the 6 m high terrace riser. The deposits at the foot of the scarp may incorporate fallen blocks of terrace deposit as suggests the high resistivity anomaly shown between the reference lines 20 and 14 in profile T1a.
(2) Variations in electrical resistivity related to changes in water content and salinity, controlled by changes in porosity and the position of the water table and capillary fringe; and (3) the probable presence of a fresh water-brine water interface at the lake margin dipping away from the basin (Samper et al., 1993;García-Vera, 1996;Yechieli and Wood, 2002). The main findings of the seismic refraction traverse (SR2) and the ERT profiles (T2 and T3) carried out in the bottom of the lake include ( Fig. 6c): (1) The unconsolidated lake deposits plus the underlying karstic residue have average v p =301 m/s and ρ b 1 Ω m. This unit has an overall planar base situated at 1-2 m below the surface, in agreement with the position of the top of the bedrock identified in three boreholes drilled in this lake. González-Sampériz et al. (2008) found the top of the bedrock at 1.62 m in a core located between profiles T2 and T3. Pérez-Obiol and Roure (1990) and Stevenson et al. (1991a,b) reached the rockhead in boreholes perforated in La Playa bottom at 1.1 and 2.3 m below the surface, respectively. The low-relief irregularities that display the base of this low electrical resistivity unit may be related to lateral changes in thickness and composition of the karstic residue. That is, a slightly irregular rockhead typical of mantled karst settings.
(2) The seismic refraction traverse SR2 reveals two seismic units within the bedrock at approximately 6 m beneath the surface, that may correspond to two zones with different degree of karstification or to two lithological units of the horizontally lying Miocene sequence. The upper one has an average v p = 1549 m/s and the lower one v p = 4290 m/s. Bulk resistivity in the bedrock ranges from ρ =2 to 12 Ω m. The higher electrical conductivity of the bedrock beneath the lake than at the margin of the depression may account for the presence of higher salinity phreatic water in the former. Profile T2 performed in the central sector of the lake, shows vertical fringes with lower resistivity attributable to zones of preferent brine discharge related to fracture and dissolution enhanced permeability.

Lake terraces. Stratigraphy
The middle terrace was investigated by means of two adjacent trenches dug in the NE margin of La Playa (E735927 N4589264; Figs. 2, 4, 7, and 8). Initially, we dug a trench by hand in the terrace scarp, extending from the top of the terrace to the base of the karstic residue situated between the lacustrine deposit and the bedrock (Figs. 7a and 8). This N70E trending trench with stepped base was 10 m long, 1 m deep and 0.6 m wide. It provided an exposure with limited lateral extension (0.5-2 m) of the 5 m thick terrace deposit. Subsequently, a backhoe trench with N160E orientation was excavated next and parallel to the terrace scarp (Figs. 7b, c, and 8). The aim of this trench was to obtain a wider picture of the upper deposits of the terrace, in order to better characterise their geometry and architecture and to increase the likelihood of finding datable material. The depth of this 11 m long trench was limited to 3.1 m due to safety reasons. This backhoe trench was connected to the hand-dug trench by a subsidiary trench in order to establish a physical correlation between the upper stratigraphic units identified on both trenches (Figs. 7c and 8). A synthetic stratigraphic log has been produced for the terrace; the lower 188 cm were described in the hand-dug trench and the upper 312 cm in the higher quality exposure of the backhoe trench (Fig. 9).
Three main units have been identified in the trenches (Fig. 8): (1) bedrock composed of slightly weathered and horizontally lying alabastrine nodular gypsum with marl intercalations; (2) a karstic residue 45 cm thick made up of highly cohesive dark grey marl with residual gypsum particles up to 1 cm long. The top of both the bedrock and the karstic residue are slightly inclined towards the lake, as it was interpreted from the ERT profiles (Fig. 6). (3) A 5 m thick lake terrace deposit, in which 25 layers have been differentiated. The uppermost unit, 39 cm thick, corresponds to lacustrine sediments disturbed by ploughing. The hand-dug trench exposed in the upper part a sediment-filled conduit 50 cm across cross-cutting several layers (Fig. 8). This feature corresponds to a biogenic burrow, probably related to a fox or rabbit den. Pollen content in the 18 samples collected in the trench is very low (less than 30-40 grains) and they contain only pine, some Chenopodiace, Poaceae and isolated grains of Artemisia and Fabaceae, preventing any palynological analysis.
The terrace deposit is primarily composed of tabular and horizontally-bedded massive or vaguely laminated gypsiferous sands and silts (Fig. 8). A significant proportion of lenticular gypsum crystals occurs in some beds, most likely precipitated in the playa-lake bottom or within the sediments from interstitial water (Pueyo, 1978/79;Pueyo and Inglés, 1987). The gypsum crystals may be in place or may have been reworked as intraclasts by eolian or wave currents. Some beds (4, 7, 17, 21, and 22) that may include granule-sized particles show an undulated erosional base recording higher competence water flows at the lake margin (Fig. 9). Layer 17 corresponds to a 50 cm wide and 8 cm deep channel filled with massive dark grey-brown sandy silt with charcoal fragments up to 5 mm long and abundant bioturbation (Fig. 8). The top of layer 3 includes centimetre-size silt rip-up clasts probably reworked from mud cracks. Layer 18 displays faint undulated lamination attributable to ripples. Several depositional sequences 11-35 cm thick can be recognised in the middle part of the terrace deposit. These sequences typically consist of light-coloured gypsiferous sand or silt with an erosional or planar base and a darker and finer-grained gypsiferous facies with abundant bioturbation at the top. The base of the darker upper part, with higher content in organic matter, may be net or gradational. The bioturbation-related porosity is largely filled by white secondary carbonate. These depositional sequences may record flooding events in the lake, starting with an erosional water flow event, a lake level increase with progressive dilution favouring the proliferation of biota in the lake bottom and the accumulation-preservation of organic matter (higher productivity and/or preservation of organic matter).
Grain size and compositional analyses also identify four main units above the rockhead, including the basal 45 cm karstic residue, the top 39 cm sediments altered by farming, and two main units within the lake sediments in the terrace (Fig. 9). The lower terrace unit (210-391 cm depth) is characterised by finer clayish silt sediment (up to 50% clay size particles) with the lowest organic (less than 1%) and carbonate (less than 1%) content and the highest sulphur content (up to 1.5%) of the whole sequence. These features suggest a sulphate-rich brine lake more conducive to gypsum deposition. The transition to the upper unit at about 210 cm depth is rather sharp and the upper sediments in the trench are dominated by coarser particles (40% clay size) with relatively higher organic (up to 1.5%) and carbonate (up to 1.5%) content and lower sulphur content (b0.5%). The upper 1 m thick unit is slightly coarser. Lower gypsum content, the relatively coarser nature of the sediments and the presence of numerous fining upward sequences are indicative of a less saline environment and the occurrence of frequent flooding episodes in the basin.
A 4 m long and 1 m deep trench with N40E orientation was dug by hand in a bench ascribed to the upper terrace located southeast of El Pueyo Lake (E736799 N4588450; Figs. 2 and 4). This terrace remnant, at about 9 m above the lake bottom, is situated at the windward edge of a bedrock yardang and is surrounded by an extensive inset surface of the middle terrace. The trench exposed massive yellowish white gypsiferous silts 75 cm in thickness overlying an erosional base carved in Miocene limestone bedrock. The deposit includes a discontinuous bed of angular cobble-size limestone clasts. This deposit is clearly more indurated (cemented) than those of the younger middle terrace. The geomorphic setting and the sedimentological characteristics suggest deposition at the lake margin affected by both lacustrine sedimentation processes and slope wash responsible for the accumulation of clasts derived from the nearby bedrock outcrop. The limited thickness of the deposit may be related to its marginal position and postdepositional erosional lowering of the terrace surface.

Chronology and morpho-sedimentary evolution
Charcoal samples from the middle terrace with enough organic carbon content for AMS 14 C dating were collected from layers 1, 11, 14, 17, 21 and 23 at 20, 221, 256, 294, 330, 385 and 405 cm above the base of the deposit, respectively (Figs. 8 and 9; Table 1). These samples yielded the following calibrated ages consistent with their stratigraphic order: 3887-3699, 3512-3383, 3638-3465, 3085-2925, 2742-2461, 2158-1995, and 487-311 cal yr BP (2 sigma calibrated age ranges with a relative probability greater than 0.80). The anomalous Table 1 Radiocarbon AMS dates for the middle La Playa terrace including label of the charcoal samples, location, laboratory and code, conventional radiocarbon ages, calibrated dates (using CALIB 6.0.1 and the data set intcal 09.14c; Reimer et al., 2009) and vertical distance from the base of the terrace deposit to the sampling point. In bold, 2 sigma calibrated age ranges used in the chronological model, with relative probabilities higher than 0.80.  10. Chronological model for the middle terrace based on AMS radiocarbon ages derived from samples collected in the trenches. Age ranges represent the 2 sigma calibrated 14 C ranges with relative probabilities higher than 0.80. The ages that do not fit the model corresponds to a sample collected from the terrace scarp at some distance from the trenches and a sample obtained 95 cm below the terrace tread. The deviation of the lower sample is most likely related to erroneous correlation or the sampling of reworked material. The anomalous modern age of the higher sample is attributed to soil contamination. modern age of the uppermost sample, collected at 95 cm below the terrace surface, is attributed to soil contamination. Relatively large pieces of charcoal were collected from an approximately 50 cm thick bed exposed at the terrace scarp at some distance from the trenches, providing a calibrated age of 3269-3060 cal yr BP. The sampled sediments were visually correlated with the upper part of layer 1. However, the large deviation between the obtained numerical age and the chronological model constructed with the better stratigraphically constrained trench-derived ages suggests that such correlation was erroneous (Fig. 10). Most likely the analysed charcoal corresponds to a younger lacustrine bed or to a colluvial deposit derived from more recent lake deposits.

Label and location
Unfortunately, we did not find any suitable material for 14 C dating in the trench dug in the upper terrace. Two samples (UTL1 and UTL2) were collected for optically stimulated luminescence (OSL) dating from the massive gypsiferous silts situated above and below the layer of limestone clasts, at about 45 and 15 cm above the base of the deposit, respectively (Fig. 11a). Details of the OSL dating procedures and equipment employed in this study are provided in the online Supplementary information.
Sample UTL2 did not yield sufficient quartz grains for equivalent dose (D e ) estimation. However, we were able to obtain finite D e estimates for 133 of the 900 grains analysed from sample UTL1 (Fig. 11b). Sample UTL1 yielded a single-grain OSL age of 128 ± 9 ka and an Fig. 11. a: OSL sampling position in the hand-dug trench excavated in the upper terrace. The two PVC tubes represent replicates of sample UTL1. Sample UTL2 was collected 30 cm above sample UTL1, but did not yield sufficient quartz grains for dating. b: Single-grain equivalent dose (D e ) distribution for sample UTL1, shown as a frequency histogram with ranked D e estimates. Field water content, expressed as a % of dry mass of mineral fraction and assigned a relative uncertainty of ±10% to accommodate the effect of any potential variations during the burial periods. c Measurements made on dried and powdered samples by low-level beta-particle counting using a Risø GM-25-5 beta counter. External beta dose rates were determined by comparing the total beta counts obtained for each sample against those obtained for a geological standard (Nussloch loess GEOPT13) having an independently determined beta dose rate. Allowance was made for the effect of grain size (Mejdahl, 1979) and HF acid etching (Brennan, 2003) on beta-dose attenuation. Beta dose rates were adjusted for moisture attenuation using present-day water contents. d Mean ± total uncertainty (68% confidence interval), calculated as the quadratic sum of the random and systematic uncertainties. e Calculated from in situ Na:Tl gamma spectrometry using the energy windows approach of Aitken (1985) (see further details in Arnold et al., 2012). Gamma dose rates were derived from radionuclide concentrations of 40 K, 238 U and 232 Th using the published conversion factors of Adamiec and Aitken (1998) and Stokes et al. (2003). f Determined using the approach of Prescott and Hutton (1994), adjusted for site altitude, geomagnetic latitude, and thickness and water content of sediment overburden, and assigned a relative uncertainty of ± 10%. g Total dose rates includes an assumed internal dose rate of 0.03 Gy/ka (Bowler et al., 2003), with a relative uncertainties of ± 30%. h Number of single grains (SG) or multi-grain aliquots (MG) that passed the SAR rejection criteria  and were used for D e determination/total number of grains or aliquots analysed. i Relative standard deviation of the D e distribution after allowing for measurement uncertainties, calculated using the central age model (CAM) of Galbraith et al. (1999). j Final burial dose estimate, calculated using the 4-parameter minimum age model (MAM-4) of Galbraith et al. (1999) based on the goodness-of-fit criterion outlined by Arnold et al. (2009). k Total uncertainty includes a systematic component of ±2% associated with laboratory beta-source calibration. associated multi-grain aliquot OSL age of 127 ± 18 ka (both calculated using the 4-parameter minimum age model (MAM-4) of Galbraith et al., 1999; Table 2). The obtained age of the analysed sample is much older than might be expected considering its morphostratigraphic position. The upper terrace is situated just 3 m above the middle terrace, so it is highly improbable that the two deposits are separated by >100 ka and there are no evidence of other Late Pleistocene or Early Holocene high-stand terraces. Moreover, the preservation of unconsolidated lake deposits at the windward edge of a yardang over a 128 ka period seems unlikely. The terrace might correspond to a strath surface developed on much older deposits. Additional discussions about the reliability of the OSL age are provided in the online Supplementary information. The mapped sequence of terraces reveals that the La Playa-El Pueyo lake basin has been dominated by an overall deepening trend throughout its evolution (Fig. 12). The terraces record interruptions in this trend at times during which the lake was dominated by aggradation. Unfortunately, we only have one OSL numerical age for the upper terrace (ca. 128 ka) whose reliability is difficult to assess independently and the chronology of the lower terrace level is poorly constrained; though it is likely to be recent Holocene in age and probably historical. Regarding the middle terrace, at the trench site, the accumulation of lacustrine deposits 5 m thick occurred from ca. 3.9 ka to soon after 2 ka, yielding a maximum average aggradation rate of 2.6 mm/yr. Using the ages of the samples collected in units 1 and 23 we can estimate an average vertical accretion rate of 2.1 mm/yr. According to the chronological model (age vs. distance to base), aggradation was much more rapid during deposition of the lower half of the sequence; 8.3 mm/yr between layers 1 and 14, and 0.9 mm/yr between layers 14 and 23 (Fig. 10). Considering the approximate area of La Playa-El Pueyo lakes during the deposition of the middle terrace (2.73 km 2 ) and assuming that 5 m of sediments were accumulated on average during this aggradation phase, we can estimate a total sedimentary input of 13.65 million m 3 . This volume, incorporated in a minimum time span of 1.9 ka, yields a maximum sedimentary input rate of 7184 m 3 /yr and a specific value of 26.31 m 3 /ha/yr. Subsequent to the accumulation of the middle terrace, the lake bottom underwent entrenchment attributable to wind deflation during periods of low water tables. This erosional lowering phase was interrupted by a brief aggradation period in recent times, probably historical, recorded by the lower terrace. Differential erosion resulted in the segmentation of La Playa and El Pueyo lakes (Fig. 12). Considering that the top of the deposits of the middle terrace is situated 5.95 m above the lake bottom and that the erosional phase started sometime after 2 ka, we can estimate a minimum mean vertical erosion rate due to wind deflation of 3 mm/yr. This minimum rate also applies to the yardangs carved in the deposits of the middle terrace. Additionally, considering the current area of La Playa and El Pueyo lakes (1.72 + 0.14 km 2 ) and the relative height of the terrace (5.95 m), we can estimate that eolian deflation has evacuated a volume of around 11.07 million m 3 in La Playa-El Pueyo lake system in a time period shorter than 2 ka. This means a maximum long-term erosion rate of 5533 m 3 /yr and a specific rate of 30 m 3 /ha/yr.
The long-term deflation rates estimated for La Playa-El Pueyo lake system compares well with those calculated in several arid regions of the world, generally using yardangs carved in Holocene lake deposits. Fig. 12. Geomorphic evolution of La Playa-El Pueyo lake system and graph illustrating the successive aggradation and entrenchment phases recorded by the lacustrine terraces. It is not clear whether the current lake bottom represents a strath cut-across the lower terrace or a fill inset in the latter. McCauley et al. (1977) indicate rates of 20 mm/yr for the Lop Nor yardangs in Central Asia. The data presented by Anderson et al. (2002) on the Laguna del Perro, a deflation basin within the pluvial Estancia Lake in New Mexico, indicate a surface lowering rate of ca. 3 mm/yr. Excavation rates ranging between 0.4 and 4 mm/yr have been reported in lakes for post-Neolithic pluvial periods (Cooke et al., 1993 and references therein). Goudie et al. (1999) estimate wind erosion rates between 0.5 and 1.7 mm/yr over the last 4.5-5 ka at Kharga Oasis, Egypt. In Sebkha Mellala, Algeria, a deflation rate of 0.4 mm/yr has been estimated for the last 6-7 ka (Boyé et al., 1978). In the last 2 ka, differential wind erosion has operated at rates of 0.5-2 mm/yr in relict mantled pediments in Algeria (Williams, 1970). Present-day deflation rates in the Bodélé Depression, Chad, a pluvial lake filled with diatomites which constitutes the planet's largest single source of dust, are of the order of 3 mm/ yr (Washington et al., 2006). In Ebinur Lake, NW China, measured contemporaneous deflation rates range from 2.4 to 9.6 mm/yr (Liu et al., 2011).

Discussion and conclusions
Two main interpretations have been proposed to explain the origin of the lake basins in Bujaraloz-Sástago area: (1) collapse of large karstic cavities developed within the evaporite bedrock; and (2) widespread subsurface dissolution and subsidence in combination with eolian deflation. The overall tabular geometry of the lake fill, with no significant thickness changes inferred from the ERT and seismic profiles acquired in La Playa, strongly support the second alternative. Once the bottom of solution sinkholes reaches the water table, local groundwater flows that discharge in the lakes cause widespread dissolution of the bedrock and gradual settling of the lake sediments. The gypsum removed from the bedrock as solutes precipitates in the lake contributing to vertical accretion. Under favourable conditions, when the lake dries out, the lacustrine sediments are evacuated by wind deflation causing the erosional lowering of the basin floor. This genetic model is also supported by: (1) the limited thickness of the deposits underlying the lake, including significant erosional hiatuses (Moreno et al., 2004;González-Sampériz et al., 2008), (2) the occurrence of terrace deposits at the margins of the lakes that do not seem to have correlative units in the lake bottom, (3) the presence of yardangs carved in bedrock and terrace deposits in the leeward side of the largest playas, (4) the rapid development of nebkha dunes and (5) the absence of significant cavities in the boreholes drilled in the area.
The studied lake basins have been dominated by an overall deepening trend during their evolution, episodically interrupted by net aggradation periods recorded by lacustrine terraces (Fig. 12). The upper terrace, situated at 9 m above the lake bottom, records the oldest accumulation phase with preserved geomorphic expression. Most likely the limited thickness of the deposit (75 cm) found in the trench excavated in this terrace is related to the marginal position of the site and to post-depositional erosion. We have obtained a single OSL age of 128 ka from this level. Further chronological control is needed to assess the validity of this age estimate.
The middle terrace is situated at 6 m above the lake bottom and at the trenching site the deposit of this morpho-stratigraphic unit is 5 m thick. This indicates that the lake bottom underwent a minimum deepening by deflation of 8 m between the accumulation of the upper and middle terraces. The middle terrace records a net aggradation phase of at least 5 m between 3.9 ka and soon after 2 ka. The lowest and upper dated layers allow the estimation of a maximum average vertical accretion rate of 2.6 mm/yr. The chronological model reveals that aggradation occurred at a much higher rate during the lower half of the sequence than in the upper one (8.3 vs. 0.9 mm/yr). The spatial distribution of the middle terrace indicates that during the accumulation of this unit La Playa and El Pueyo lakes were integrated into a larger single lake.
The subsequent excavation period was interrupted by a brief aggradation phase recorded by the lower terrace, situated a 0.5 m above the lake bottom. During this major erosional phase, differential eolian deflation excavated the depressions of La Playa and El Pueyo lakes, separated by remnants of the middle terrace. This involved a reduction of the lake area from 2.73 km 2 to 1.86 km 2 . Considering that the erosional phase started soon after 2 ka, we can estimate a minimum mean vertical erosion rate due to wind deflation of 3 mm/yr. This minimum rate also applies to the yardangs carved in the deposits of the middle terrace. The apparently high long-term deflation rate estimated for La Playa-El Pueyo lake system compares well with those calculated in several arid regions of the world using yardangs carved in Holocene lake deposits, as we have mentioned above (Washington et al., 2006 and references therein).
The changes in the morpho-sedimentary behaviour of these internally drained playa-lakes depend on the relative contribution of the processes that lead to accumulation of sediments in the lake bottom and those that cause its deepening. The processes that determine the complex mass balance in these basins are largely controlled by the presence of evaporite bedrock, groundwater flows that discharge in the lakes and strong winds with great geomorphic effectiveness, as reveals the presence of yardangs.
Vertical accretion in the lakes may be related to two main types of sedimentary processes: (1) autochthonous deposition related to evaporitic precipitation of salts from solutes largely supplied by groundwater flow; and (2) allochthonous deposition of extraclasts transported by surface runoff and the wind. La Playa and El Pueyo terminal lakes, covering 1.86 km 2 , act as the ultimate sediment sinks for their catchment area of about 35 km 2 . Additionally, air currents may undergo a reduction in velocity below the transport threshold due to flow separation in these topographic depressions, leading to accumulation of wind-transported particles.
Deepening of the playa-lakes may be caused by two processes: (1) Groundwater flows that discharge in the lakes may cause dissolution of the bedrock with the consequent subsidence of the lake bottom. This process, albeit continuous, acts at a low rate, since the highly mineralised underground water that discharges in the lakes has a limited karstification capability. (2) Eolian deflation by winds with velocities higher than the entrainment threshold. This process can only operate when the surface of the playa is vegetation free and dry and there is material susceptible to be entrained by air currents. Several processes such as drying cycles, salt weathering, trampling and biogenic burrowing predispose particles to wind erosion (i.e. Goudie and Wells, 1995). Lowering of the lake bottom by eolian deflation is limited by the position of the water table and the capillary fringe (elevation to which the groundwater ascends by capillary force), controlling the development of extremely flat surfaces (Stokes, 1968;Rosen, 1994;Yechieli and Wood, 2002). Once the surface meets the water table, the basins cannot get any deeper, but expand laterally through backward retreat of their margins. According to Rosen (1994), in order to have net aggradation in a playa, the surface must be flooded more often than it is dry.
The relative contribution of the processes listed above may be significantly altered by climate and land use changes. Under more arid conditions the playas remain dry during longer periods, favouring the lowering of the surface by deflation. Conversely, lower water availability in semiarid areas may involve a significant reduction in the vegetation cover, with the consequent increase in sediment flux by runoff and deflation; that is higher detrital input to the terminal lakes (Knox, 1984;Bull, 1991). A similar impact may be expected from anthropogenic reductions in the vegetation cover (deforestation, fires, and overgrazing). In contrast, a denser vegetation cover and more frequent flooding under a relatively more humid climate would hinder both external detrital input and eolian deflation. A wetter climate would also favour dissolution-induced subsidence. According to this conceptual model, a change to drier conditions could lead to net erosion or aggradation depending on the relative role of the different processes on the mass balance, mainly deflation vs. external detrital input. However, it is widely accepted that major deepening phases in playa-lakes and pans occur during more arid periods, when the bottom of the depressions is exposed more frequently to wind deflation and the conditions are more favourable to the preparation of wind-exportable particles (i.e. Goudie and Wells, 1995;Gutiérrez-Elorza et al., 2005;Washington et al., 2006). This is also our preferred interpretation for the studied lakes based on the following circumstances: (1) The middle terrace records the existence of a single lake 47% larger than the present-day La Playa and El Pueyo playas, suggesting a more positive hydrological budget. (2) Although La Playa-El Pueyo lake basin has a large catchment area, the connectivity (coupling) between the slopes, drainage network and terminal lakes is rather limited due to the dominantly flat topography of the platform, dissected by poorly ranged valleys with low topographic gradient (Fig. 4). Higher sediment yield from the slopes could involve increased aggradation in the valleys but limited sediment transfer to the lakes. (3) Runoff has a very limited efficiency within the lakes due to the extremely flat topography, hindering sediment redistribution from the margins towards the centre.
A number of publications present chronological information on Holocene aggradation and incision phases in ephemeral valleys within the Ebro Depression (Burillo et al., 1985;Cuchí and Soriano, 1991;Arauzo and Gutiérrez, 1994;Macklin et al., 1994;Gutiérrez and Arauzo, 1995;Peña, 1996;Peña et al., 2004;Harvey and Gutiérrez-Elorza, 2005;Sancho et al., 2008) (Fig. 13). The variable chronologies suggest considerable spatial variations in the timing of the aggradation/incision sequences, most probably related to the superposition of climatic and local factors (human activity, extreme hydrological events, local base-level changes and hillslope-channel coupling conditions; i.e. Constante et al., 2010) on the response of these highly sensitive systems. The available data indicate that during the period of accumulation of the middle lacustrine terrace, some ephemeral fluvial systems were dominated by both aggradation and incision processes. There seems to be a temporal cluster of aggradation phases (2.4-1.8 ka BP) correlative to the uppermost part of the middle terrace.
A recent generation of talus flatirons has been AMS 14 C dated at three sites in the Ebro Depression, recording an accumulation phase on the slopes between 3 and 2.5 ka BP. This aggradation period, contemporaneous with the lower half of the middle lake terrace, is attributed to warm/ wet conditions that favoured an increase in the vegetation cover preventing incision on the slopes (Gutiérrez-Elorza et al., 2010). This interpretation is in agreement with a more positive hydrological budget in the lakes, with more frequent flooding preventing deflation. The aggradation period recorded by the middle terrace of La Playa partially overlaps with a phase of increased frequency of large magnitude floods inferred from slackwater flood deposits (2880-2430 cal yr BP) and a cluster of dates from alluvial floodplain deposits (2750-2150 cal yr BP) identified by Thorndycraft and Benito (2006) in Spain.
The onset of the Mid Holocene period of higher aridity in the Iberian Peninsula at around 6 ka BP correlates with decreasing moisture in the Fig. 13. Dated aggradation phases in ephemeral valleys within the Ebro Depression. Roman numbers correspond to relative ages based on geoarchaeological criteria. Arrows and question marks indicate poorly constrained onset and/or resumption of accumulation period.
Alps (Magny et al., 2006) and the end of the African Humid Period in northern Africa (Kropelin et al., 2008), and it has been interpreted as a reflection of decreasing summer solar radiation in the northern hemisphere. The 4.2 ka arid event has been characterised in eastern (Weiss et al., 1993;Arz et al., 2006), central (Drysdale et al., 2006;Magny et al., 2007) and western regions (Martín-Puertas et al., 2009) as a period of intense droughts in the Mediterranean that ended the long Mid Holocene arid period and inaugurated the relatively more humid late Holocene.
The review of available lacustrine, pollen and archaeological records in the Central Ebro valley by González-Sampériz et al. (2008) identified the early Holocene as the more humid period during the Holocene, contrasting with the aridity of the Middle Holocene interval, and a moisture recovery during the late Holocene. Mid Holocene deposits are frequently absent in the Central Ebro valley lakes as a result of intense eolian erosive processes. The lack of archaeological remains associated with the Middle Holocene (Neolithic) also supports increased aridity that would have impeded human settlement. Palynological sequences from archaeological and lacustrine sites in northeastern Spain have documented significant changes in vegetation during the last millennia, although in most cases they have been ascribed to human impact more than to climate or hydrological changes. A list of the most significant pollen sites in NE Spain includes: Cova Punta Farisa (Burjachs, 1993), Majaladares Cave (López-García, 1992), Salada Mediana (Valero-Garcés et al., 2000a), Laguna Guallar (Davis, 1994), La Clota, El Rebollón and El Camerón (Pérez-Obiol and Roure, 1990), La Playa and La Salineta , Salada Pequeña and La Estanca (Davis, 1994), Gallocanta (Davis, 1994;Burjachs et al., 1996;Roc et al., 2002), Tozal de Macarullo, Tozal de Andrés and El Prao (González-Sampériz and Sopena, 2002), Taravilla Lake , Cabezo de la Cruz (Iriarte, 2009), Estanya Lake  or El Cabo and La Cabañeta (González-Sampériz, 2004) sites, among others. Although the chronological resolution is different for these records, all of them show a recent period characterised by: i) less arboreal pollen content (AP) dominated by pines, some junipers and evergreen oaks, interpreted as a result of deforestation processes; ii) the beginning of continuous curves of Olea, Cerealia and/or Vitis and the presence of isolated grains of Castanea or Juglans, showing the development of cultivated areas; iii) the increase of nitrophilous and herbaceous taxa, frequently related to pastoralism and agricultural activities (Brassicaceae, Cichorioideae, Asteroideae, Plantago, Rumex, Chenopodiaceae, etc.). Thus, although the main patterns of the recent Holocene vegetation landscape evolution at a regional scale are known, the investigated sequences lack sufficient chronological and geographical resolution to reconstruct a more precise palaeoclimatic scenario for the last millennia.
Large palaeohydrological fluctuations during the last 4 ka have been described in numerous sequences in the Iberian Peninsula. Lacustrine (e.g., Zoñar Lake: Martín-Puertas et al., 2008, Taravilla Lake: Moreno et al., 2008, Estanya Lake: Morellón et al., 2009, Redon Lake: Plá and Catalan, 2005, marine (Frigola et al., 2007;Martín-Puertas et al., 2010) and fluvial records ) document a large regional climate variability and millennial scale moisture fluctuations, particularly the mid Holocene arid period with the 4 ka BP crisis, the cold Iron Ages followed by the Iberian-Roman Humid Period, and the dry and warm Medieval Climate Anomaly (Moreno et al., 2011, in press;Morellón et al., in press). At a millennial scale, the playa-lake records in the Central Ebro Basin have shown the occurrence of an arid period prior to 4 ka and a recovery of average saline lake levels during the last millennia (Valero-Garcés et al., 2000a,b,c, 2004González-Sampériz et al., 2008). But, as we have described above, the large palaeohydrological variability during the Late Holocene detected in other Mediterranean and Iberian regions has not been identified in the Central Ebro Basin lake sequences because they do not have the adequate chronological resolution.
Both, La Playa middle terrace and La Salineta terrace document large palaeohydrological fluctuations in the lake basins with intense and very significant periods of aggradation and excavation during the last 4000 years. However, the correlation between both sequences is not possible due to the absence of a more constrained chronological framework for La Salineta terrace. The occurrence of an aggradation phase in La Playa Lake after 3.9 ka is coherent with a number of Iberian and Mediterranean records that show more humid conditions after the 4.2 ka event. In Montcortès Lake (Lleida, NE Spain), a humidity increase occurred between 3800 and 2350 cal yr BP , similar to that at Zoñar Lake (Córdoba, S Spain) after 3.5 ka (Martín-Puertas et al., 2008. This period corresponds with the cold and wet phase of the Iron Age in northern Europe (Leira, 2005) between 3000 and 2000 cal yr BP (1050-50 BC), marked humid conditions in some alpine lakes (Magny et al., 2007) between 3950 and 3850 cal yr BP, the onset of the humidity increase in the Mediterranean area (Sadori et al., 2004) and evidence of higher flood activity (2865-2350 cal yr BP) in the Spanish fluvial record . This period includes the Bronze Age (4-2.7 ka BP) and the first part of the Iron Age, when several areas of the Ebro Depression show an intense occupation pattern (some caves but more frequently open sites in lowlands and middle altitudes) and therefore, a high population density (Utrilla and Rodanés, 1997), probably in response to a more humid period.
The erosional phase that occurred after ca. 2 ka is coherent with relatively more arid conditions in the Iberian Peninsula after the Roman Period and also during the Medieval Ages (Moreno et al., 2011, in press): in Lake Montcortès, a lake level drop between 2350 and 1910 cal yr BP (400 BC-40 AD) has been inferred by  coinciding with a reduction in river flood frequency on the Iberian Peninsula between 2350 and 2000 cal yr BP (Macklin et al., 2006) and lower lake levels in Zoñar Lake between 2100 and 1900 cal yr BP (Martín-Puertas et al., 2008).