This work is close beyond the common applicability of the sequence stratigraphy principles. Nevertheless one aim is to demonstrate its validity for a small ice-dammed palaeo-lake. Since base-level variability is smoothly related with an oscillatory motion, the chronostratigraphy chart emerges as a crucial tool to produce enough quantitative data for a Fast Fourier Transform (FFT) analysis, even if this chronostratigraphy is based on local stratigraphic constraints. The exposed case is related to sixth order stratigraphic cycles, even higher, related to high-frequency global events occurring at the end of the Pleistocene. The glacio-lacustrine record was controlled by an intermittent ice-damming from the main glacier located in Andorra. Tributary glacier tongues advance and retreat accordingly to the main glacier, however asynchronously. The system Tract evolution and the related unconformities reflect the glacier motion imprint. The glacier advances produce an unconformable subglacial till, but also a repeated TST-HST sedimentary evolution over it and their related conformable surfaces. The long-term base-level evolution was leaded by axial tilt cycle, however precession did not appear at the FFT analysis, indicating a strong high latitude climatic influence on the palaeo-lake’s evolution. The sub-orbital periodicities have been observed in the range of the Heinrich events, but also in the range of the Mediterranean salinity anomalies (periodicity of 3–5 ky).
- Andorra palaeo-glaciers
- ice-dammed palaeo-lake
- orbital cycles
- Heinrich events
- Mediterranean moisture
The used approach allows to use the principles of sequence stratigraphy in order to study the Valira del Nord sedimentary infill, which has been controlled by the activity of the valley glaciers.
The unstable glacier motion is identified from its footprint in the sedimentary record. This non-stability was forced by global events . Thus, the derived sedimentary changes can be then considered as a climate proxy. In such a case, the sedimentary infill is ruled by sixth order stratigraphic cycles  or higher frequency global events, like Heinrich events. In such natural phenomena the icebergs break off from the glaciers crossing the Northern Atlantic Ocean influencing the circulation patterns of thermohaline currents.
The sedimentary infill from the La Massana has been assimilated to eustasy cycles produced in a passive margin basin. In such a case the deltaic build-up depends exclusively on base level changes. In the case of La Massana, an ice-dammed palaeo-lake was produced at each significant glacier advance.
The mapping resolution derives from available profiles (about 232 profiles) par basin surface (about 11 km2), here close to 50 profiles par km2. Stratigraphic data is accessible in an online database . The chronological resolution derives from previous published dating , from about 17 absolute dating in a span of time of 46 cal ka (one par 3 ka).
Previous woks  supposed a regular timing of the depositional sequence build-up, however this assumption was supported with too few datings. Although the impact of global events in continental areas, like Heinrich events, was still not well known . However increasing data, especially geochemistry and absolute datings, a first correlation to the global events could be performed by [6, 7]. Nowadays, since the chronostratigraphic framework has been settled at La Massana  it is possible to assess numerically the impact of the global events.
This work aims to explore about the main forcing causes that ruled the system tracts deposition in this ice-dammed palaeo-lake.
The Principality of Andorra (42°30′N, 1°30′E) is located in the southern slope of the eastern Pyrenees (Figure 1). The height of about one half of its surface (233 km2) exceeds 2000 m a.s.l., and the landscape has been mainly controlled by the glacial geomorphology and by glacial geology. The tributary valley system has a “Y” shape, with two main basins joining in the city of Andorra (1000 m a.s.l.). The highest peaks occur in the Valira del Nord valley, the NW tributary (145 km2), reaching 2942 m a.s.l. (Coma pedrosa Peak), not far from the first easternmost 3000 peak from the Pyrenees, the Pica d’Estats peak (3143 m a.s.l.) on the Franco-Spanish border. In former times different glaciers coming from the Valira del Nord valley converged with those coming from the NE valleys (211 km2), and form the main glacier of Andorra (Figure 2). The NE tributary embraces the Valira d’Orient (166 km2) and the Madriu valley (45 km2), the last one listed in the UNESCO’s world human heritage. This main glacier overflowed the Spain-Andorra border (838 m a.s.l.) at the Upper Pleistocene .
The studied area is located in the lowest part of the Valira del Nord tributary (Figure 3), involving the village of Erts (42° 33′ 41′′ N–1° 29′ 48′′ E, 1360 m a.s.l.), and the towns of La Massana (42° 32′ 44′′ N–1° 30′ 53′′ E, 1230 m a.s.l.) and Ordino (42° 33′ 19′′ N–1° 31′ 60′′ E, 1300 m a.s.l.). The Valira del Nord population (15,000 people) lives in 1% of its basin surface and urban zones lack in space. Understanding the stratigraphy and the extent of the palaeo-lake sediments is quite relevant for Engineering-Geologists who work in La Massana. Glaciolacustrine deposits produced during the Last Glacial Period influence mainly on building bearing capacity and hazards of excavations. For this reason when the spatial distribution knowledge is improved, the variability and predictability nature of sediments is reduced. Geological knowledge improves the safety factors to apply on such projects and reduces costs. No-negligible financial efforts are actually employed for geotechnical surveying, what benefit the sequence stratigraphy approach by the available profiles. Fortunately or unfortunately, it depends on the lecturer’s opinion, there is no any other region in the southern slope of the Pyrenees with a comparable building pressure than in Andorra, making La Massana palaeo-lake a case unique in this region.
However, unique are also other ice-damming cases in the Pyrenees, like Llestui, Llauset, Cerler, Linás de Broto [9, 10], or in their respective mountain massif, like Pias in NW Spain , but in none of such cases had been correlated to global climatic events previously as in La Massana .
La Massana palaeo-lake was dammed by the largest glacier from Andorra at MIS 3 , long time after the maximum ice extent (MIE) in Andorra (Figure 4). The MIE is previous to 59 ka, just before a widespread glacier recession that starts before 40 ka cal BP . The maximum ice extent (MIE) in the Pyrenees is diachronous to the global last glacial maximum (LGM) of the northern hemisphere ice sheets [14, 15, 16, 17, 18].
This large period of glacier recession in the Pyrenees at the MIS 3 start is a classical boundary from the local würmian stratigraphy of prehistoric excavations, the würm II - würm III boundary [19, 20]. A general glacier advance before 32 ka occur in Andorra , also in the SW Pyrenees , and in the inner Iberia could be the beginning of glacier tongues presence . The starting point of the glaciolacustrine sedimentation in La Massana palaeo-lake is related to that general glacier advance at MIS 3. However all glaciers recede again short after . Only the SE Pyrenees glaciers recover their original extent at MIS 2 , but those from the SW Pyrenees not .
Two different but complementary disciplines of the sequential stratigraphy can be clearly differentiated: 1) the analytical way and 2) the synthetic way. The last of both aims to develop a time scale of global changes that can be inserted into the global chronostratigraphic scale. The first one, that is the analytical concept, is understood as stratigraphic modeling of facies associations.
In any cases the main goal is the recognition of units (set of strata) within the basin infill materials, limited by surfaces that mark a genetic change conditions affecting the entire basin. Its recognition is a mandatory object in all basin analysis.
The applicability of the sequence stratigraphy to large ice-dammed palaeo-lakes has been settled since the end of the last century . Shortly after sequence stratigraphy principles were used in Andorra . In such works the sequence stratigraphy terminology is widely applied such as depositional sequences, system tracts and parasequences. However, since theses pioneer works the use of the sequence stratigraphy terminology to describe the sedimentary architecture from glaciated valleys in Iberia starts to be common [1, 8, 17, 26].
The term of depositional sequence were introduced by [27, 28], and corresponds to a relatively concordant stratigraphic succession of genetically related strata, where the base and the top are discontinuities, as well as their correlative surfaces though the basin. However, the essential criteria for the recognition and subdivision of sequences were established by . Two kind of sequence limits are established by discontinuities (SB 1 and SB 2, Surface Boundary), what define also the depositional sequence type.
The System Tracts were introduced by [27, 28] to name the set of contemporaneous depositional systems formed under the same base level conditions. The System Tracts sedimentary models were established according to their transgressive retrograde character (TST, transgressive system tract), regressive prograde (FSST, falling stage system tract & LST, lowstand system tract), transgressive or prograde character (HST, highstand system tract).
A set of relatively concordant strata limited by flooding surfaces are parasequences . Three types of parasequences can be distinguished: (1) Those that had a progradational character, (2) those that are retrogradational, and (3) those that had an aggradational character. Such sedimentary build-up is according to the material supply contribution, but also in function to the base level evolution.
Since this approach is based on the magnitude of base-level changes and related boundary formation (independently to the cycle span), sequence ordering are somewhat relative. Ordering is based on the principle that a sequence cannot involve a boundary that has an equal or greater magnitude than the one at the start of the sequence , but also the magnitude of base-level changes leads the hierarchical rank of the bounding surface (SB, sequence boundary) and unconformities (US, unconformity surface), so temporal and spatial scales are not relevant , and thus applicable to this case of study.
Since glacial valleys form by external processes, the sedimentary infill at their troughs, can be assimilated to the one produced in a passive margin . Thus, only type 1 discontinuity can able to be produced in constrained glacial valleys like La Massana. Sequences starting on such kind of discontinuity had three system tracts, a FSST-LST, the TST and HST. However for large ice-dammed palaeo-lakes type 2 discontinuity has been announced by .
In the case of constrained ice-dammed palaeo-lakes from valley glacier, the glacial tongues were small comparing to the mountain range rot (>25 km) , what inhibit quick glacio-isostatic rebound and vertical accommodation, reducing the possible surface boundaries to an unique case, the type 1 case. For large ice-dammed lakes the principles are the same, however several sedimentary assemblages types can coexist , in parallel to sedimentary accommodation  or uplift (glacio-isostatic rebound).
A particular unconformity has been defined  related to the glacier motion. When the glacier tongue overrode its own proglacial deposits, deformation-compaction and eventually erosion would be produced . In such a case unconformity and diamicton are synchronously produced (Figure 5), and thus sequences that include diamictons had two supplementary system tracts, the TST and the HST .
3.3 Particular assemblages
Glaciolacustrine deltas in ice-dammed valley glaciers [1, 26] are commonly of Gilbert type . Their bottomset is much larger than the foreset and topset, what mean that most of the available outcrops are related with dense turbidites and rhythmites (Figure 6). For that reason the study of fine grain-size facies and their relative position are important. Ice-dammed lakes sedimentation is dominated by turbidity currents when the ice dam seals the end of the valley [36, 37].
If the damming increases an induced transgression occurred by ice thickening, then lacustrine facies migrated upward into the valley. Glacier induced unconformity drive base level to shift, similarly as by an eustacy cause . Once the ice ends to grow, the base level remains stable and delta progression is due by large outwash from meltwaters [17, 26]. Mass flows and grain flows on the foreset slope are then produced . Inversely graded gravel and sand facies were found by  filling scours from previous mass flow channels, similar to the high stand turbidites , but also lenticular bodies that had Bouma’s truncated series  were described at the bottomset .
If the damming ceased by glacier recession, the base lake level drops producing very marked erosion and incisions . Debris flows and denser turbidites are then deposited in the lowstand fans .
4. Sedimentary architecture from La Massana palaeo-lake
Two deltas have been distinguished (Figure 6), the Hortals and Erts deltas, forming an staircase system that works sequentially. The sedimentary build-up of the Hortals delta trunk the bottom set of its tributary’s delta (Erts), however truncation was unstable by the glaciers front oscillations. When the local glacier stabilized its front at La Massana palaeo-lake, the Hortals delta gets feed with outwash sediments coming from the glacier’s meltwaters. When the Ordino glacier advanced over the Hortals delta its tributary (Arinsal valley) was then flooded and a highstand proglacial delta was built at the locality of Erts, feed by meltwaters of Arinsal glacier . Glacier recession is produced in the reverse order, the last to come was the first to quit the palaeo-lake’s margin.
Almost seven depositional sequences (SD) have been described in the valley of Arinsal , each divided from the others by type 1 discontinuities (Figure 7). In all cases, the nature of the discontinuity leads by subaerial exposure and colluvium deposits (FSST). In some of the recessional phases, the LST local base level could lead by the main glacier of Andorra, the Valira d’Orient glacier. Then coarse turbidites outward from the prodelta feeding with sediments the bottom set (SD2 and SD4) .
In all depositional sequences the TST starts when the ice volume in the Valira d’Orient trunk glacier increase and flood its tributary valley. As damming persisted, parasequences became retrograde and low-density turbidites where produced. Ice-damming steady-state and maximum flooding surfaces (msf) are both linked, regarding the start of the HST. Extensive rhythmites are then deposited in the palaeo-lake. However if damming progresses a higher base-level reach in the palaeo-lake, subglacial till sedimentation also progresses by the advance of the local glaciers. This subglacial till is a sedimentary unconformity (Unconformity Surface, US) correlative to a transgressive surface (TS), or maximum regressive surface (mrs), promoting the start of a new TST .
The two starting sequences (Figure 7) had relative large sedimentary hiatus starting abruptly in a Forced Stand System Tract (FSST). That happen also at the start of the third and fifth sequences (Figure 7), however the time span is lesser at the start of the second and the forth sequence (Figure 7).
The relationship between water-depth changes is a complex matter, and not always elucidated . Being aware to this, only absolute height of the base-level (meters above the sea level) is here considered. When applying this to the general trend of the palaeo-lake base-level for all HST and all LST (Figure 8) a normal regression evolution  is identified at both. From this is possible to infer that deltas topsets grow faster than their bottomsets , or differential flooding is forded by the glaciers (local ice-tongues melt faster). Either cases promotes the formation of long clinoforms through time [40, 41], resulting in thin bottomsets and thick delta fronts (foresets and topsets).
Several TST and HST had been distinguished  within the sequences. Three intra-sequence unconformities (US) are identified regarding the glaciers motion. The trunk glacier (Andorra) produces the first unconformity (USa), and its correlative surface is a maximum regressive surface (mrs). The biggest local glacier (Ordino) produces a second unconformity (USb) when it advance, and its correlative surface is a transgressive surface (TS). The last unconformity (USc) is produced when the Arinsal (the smallest glacier) advance into the palaeo-lake. However at the younger sequences (SD7 and SD6) the Arinsal glacier did not reach the palaeo-lake, only USa and USb were formed .
5. Depicting base-level rates
Low-order polynomial curves had been used to fit long-term trends in time series data, like first-order functions (linear) and second-order polynomial quadratic type function:
Wich is a parabola involving a single maximum value. Progressive increase in complexity, nth-order polynomial equation had n-1 maxima and minima values, that can be used to extract trends from a log evolution (Figure base level), like the palaeo-lake base-level. The rate of the base-level change for a given time-step is given by the first derivative of the polynomial Equation , while the second derivative the variation of the rate for a given time-step (acceleration).
For La Massana glaciolacustrine record two maxima are almost needed to ensure the general trend for the palaeo-lake base-level, one at the MIS 2 and the other one at the MIS 4 related to the MIE glacial phase . The minimal base-level happen when the general glacier recession occur at the beginning of the MIS 3 . Then a third-order polynomial equation (n-1 = 2) must be invoked first (Figure 9). The first derivative (base-level rate) is a parabola (Figure 10), so a single maximum (a minimum in this case) is expected. The second derivative (variability of the rate) is here a linear function (Figure 10) and the acceleration is proportional to its inclination.
A Massana the sense of the acceleration change at 32 ka to a negative decrease, in fact acceleration increases. Thus the base-level rate (velocity) starts to increase, close to the MIS 3 – MIS 2 boundary.
6. Depicting cycles from the sedimentary record
The chronological differences between sequences are related to a high frequency cycles  in the palaeo-lake. Time-based hierarchy emphasizes here general orbital mechanisms driving the existing glacier to dam the valley, however the damming intermittency is ranging suborbital frequencies .
Fast Fourier Transform spectral analysis from the base-level rate (Figure 9) exhibit significant correlation to the orbital 41 ka axial tilting period (Figure 11), which affect the total annual solar radiation at high latitude (Imbrie & Imbrie, 1986). However, Andorra’s latitude is not so high, but the main trend of the ice-damming acts as it was. Equinoxes precession are not present in the FFT spectra, a 23 ka time period orbital cycle , affecting the total annual solar radiation at low latitude . This would imply a strong climate forcing from higher latitude climate over Andorra.
Faourier Transform spectral analysis from the base-level itself (Figure 7) exhibit significant correlation to sub orbital periods (Figure 12). The 7.4 kyrs period would be related to Heinrich events , which affect the North Atlantics polar front position. Higher frequency cycles with a period of 4.1 kyrs and 5.3 kyrs are also depicted from the FFT base-level data, linked to strong dry westerly winds blowing over the Gulf of Lion and the strength variability of the western Mediterranean stormy tracks .
Theiterranean can produce regional moisture increases by itself, most likely focused on winter. WMDW (Western Mediterranean Deep Water) formation is switched-on by outbreaks of the cold and relatively dry Mistral over the Gulf of Lions. Then Bernoulli aspiration of deep waters outflows through the Strait of Gibraltar . WMDW is produced in winter, however it may not occur every year . Nevertheless WMDW formation under extreme climate forcing outside the Gulf of Lion  could explain the paleoclimatology of the Eastern Pyrenees.
Since a strong climate forcing from high latitudes rules the Andorra (SE Pyrenees) palaeo-climate, a feedback response may exist within the Mediterranean sea. This high latitude influence must be related with the production of strong westerly winds . Westeries blowing over the Gulf of Lion enhance cyclonic circulation, inducing the WMDW to form in the Balearic basin. Within the gyre centre, mixing convective plumes reach great depths (−2000 m) where the net vertical transport is zero, acting more likely mixing elements . Divergent Ekman upwelling promotes pycnocline shallowing between the surface and intermediate waters . A violent mixing phase starts with an intense evaporation and cooling diminishing the density gradient between intermediate waters in the western Mediterranean (WM), and these from the surface . Quickly the mixed waters sink and spread horizontally at great depth into the Balearic basin what form the WMDW . Once the WMDW reach the Alboran Sea, the WMDW is a westward current along the Moroccan coast promoting evaporation and cooling . Mediterranean Outflow through the Strait of Gibraltar involves dense deep outflow from the Strait of Sicily coming from the eastern Mediterranean (TDW, Tyrrhenian Deep Water) .
Progressive surface buoyancy increase reduces new deep-water formation, what might cause reduction of Bernoulli aspiration and potential Mediterranean deep-water stagnation . WMDW properties change only due to diffusion (low Bernoulli aspiration) producing salinity anomalies in WM  when WMDW salinity is higher than in the MIW (Mediterranean Intermediate Water). Between the middle Holocene and the LGM the modeled salinity anomalies appear in a periodicity of 3–5 kyrs . Such periodicities are in great accordance to the high frequency cycles observed in the La Massana’s palaeo-lake FFT spectra.
However, variability to such general WM pattern increases by massive external water incomings, such as deglaciation phases or levantine monsoon influences. In the first case cold Alpine melt-water flow to the Gulf of Lion enhancing salinity differences causing remarkable changes in the Bernoulli aspiration . Levantine monsoon flooding reduces significant differences in salinity and WMDW strength, organic layers deposition in the WM area are then favored when deep-water injections are inhibited . Both cause Mediterranean sea-level increases and in turn a great Gibraltar’s strait opening. This effect favor salinity contrasts reduction and low exchanging velocities by continuous freshening of Atlantic inflow .
Larger valley would produce larger glacier as its de case for the Andorra glacier , but this cannot explain a faster ice growing at the smaller valleys. A sequential westward advance of the glaciers is observed  as follows: the Valira d’Orient glacier firstly advance and trunk La Massana, then the Ordino’s glacier advance and trunk the Arinsal valley, and at the end the Arinsal glacier advance. Reverse glacier motion also occurred, in an eastward glacier recession. Glacier ongoing-outgoing motion strongly influenced the glaciolacustrine record. Last In-First Out glaciers motion produce a base-level pile-up sequence (a LIFO sequence), in where TST-HST repeat several times in a single depositional sequence.
An eastward valley trend of snow owing to easterly moisture incomes of the Mediterranean sea may explain such westward motion of glaciers. Persistence of snow deflation experienced on the windward slopes  would explain a faster ice-grow of the eastern glacier of Andorra. This could be the reason why the Andorra glacier block out always its tributary valley first, even several millennia before than the local glaciers advance until the vicinity of the palaeo-lake. Moisture incomings could have a strong Mediterranean influence  at the beginning of each sequence, drifting to the west and producing the observed intra-sequence unconformities evolution (LIFO sequence).
The starting time-period of such dynamics is around 32 ka, when negative acceleration and base-level rate increase (Figure 10). This coincides when D-O events  become less frequent [8, 17]. An increase of the base-level rate entails also better snow accumulation, just in a steady decline of insolation . Long-term evolution of the palaeo-lake base-level is leaded by the 41 kyr orbital cycle.
Precession did not appear at the FFT spectra, somewhat surprising since Andorra’s latitude is not so high. This might indicate a strong high latitude climate influence in the palaeo-lake’s evolution.
Suborbital periodicities are also observed in the palaeo-lake base-level FFT spectra, in the range of Heinrich events , but also in the range of Mediterranean salinity anomalies (3–5 kyrs). Such anomalies are linked to the Mediterranean Levantine monsoon/flooding that reduces the Western Mediterranean Deep Water strength (WMDW). This inhibits the Tyrrhenian Deep Water (TDW) to advance toward the western Mediterranean over the Strait of Sicily, and the WMDW are not able to reach the Alboran sea. However only during persistent strong westerly winds situation blowing on the Gulf of Lion enhance the WMDW current, inverting the stagnation situation of the western Mediterranean.
Dry mistral winds coming from the Atlantics crossing over Iberia to the Mediterranean force the general retreat of the glaciers from the Pyrenees, almost those from the southern slope. Type 1 surface boundaries are then produced in the La Massana palaeo-lake. Ice-damming and system tracts progression could happen only when persistent moisture incomings feed Andorra, apparently in opposition to a low WMDW. However any Mediterranean moisture influence over the SE Pyrenees does not progress with the same intensity toward the Western Pyrenees .