<figure><img src="image_1.png" /></figure>

HAL
open science

# Origins of thousands of small enclosed depressions in Southwest France: an integrated geomorphological, geostatistical and geophysical study

Léa Bussière, Myriam Schmutz, Alain Dupuy

► To cite this version:

Léa Bussière, Myriam Schmutz, Alain Dupuy. Origins of thousands of small enclosed depressions in Southwest France: an integrated geomorphological, geostatistical and geophysical study. 2024. [hal-04615553](https://hal.science/hal-04615553)

HAL Id: hal-04615553

https://hal.science/hal-04615553

Preprint submitted on 18 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# A karst rather than periglacial origin for small enclosed depressions of the Landes Triangle, southwest France

Léa Bussière, Myriam Schmutz, Alain Dupuy

► To cite this version:

Léa Bussière, Myriam Schmutz, Alain Dupuy. A karst rather than periglacial origin for small enclosed depressions of the Landes Triangle, southwest France. Earth Surface Processes and Landforms, 2025, 50 (6), 10.1002/esp.70091. hal-04615553v2

HAL Id: hal-04615553

https://hal.science/hal-04615553v2

Submitted on 3 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Graphical Abstract

# Origins of thousands of small enclosed depressions in Southwest France: an integrated geomorphological, geostatistical and geophysical study

Léa Bussière, Myriam Schmutz, Alain Dupuy

<figure><img src="image_2.png" /></figure>

Small enclosed depressions
in the Landes Triangle

Findings

Origin of the depressions:

● karst processes (n=1366)

● aeolian processes (n=976)

● unknown (n=90)

Region with high karstification potential (in the groundwater recharge area)

Legend

Villagrains-Landiras anticline axis

Coversand (Landes Triangle)
# A karst rather than periglacial origin for small enclosed depressions of the Landes Triangle, southwest France

Léa Bussière<fnref n="1" />,<fnref n="2" /> ID | Myriam Schmutz<fnref n="1" /> | Alain Dupuy<fnref n="1" />

## Correspondence

Léa Bussière, EPOC, CNRS-UMR 5805, Bordeaux-INP, 1 Allée F. Daguin, 33607 Pessac, France.

Email: lbussier@uni-heidelberg.de

## Funding information

IdEx Bordeaux; Sentinelle Nord, Université Laval; Université de Bordeaux; Université Laval

## Abstract

Small topographic depressions offer insight into past and present groundwater recharge processes. The Landes Triangle (SW France) is dotted with over 2,400 such depressions, whose origins and hydrogeological implications remain unclear. This paper aims to disentangle their origin, with emphasis on the Villagrains-Landiras anticline sector: a key recharge zone, which also features the highest density of depressions. We first evaluate three historical hypotheses with open GIS data at the regional scale, then detail the near-surface structure of a depression near the anticline with electrical resistivity tomographies and ground-penetrating radar. We identify two groups of depressions based on their morphology, distribution and geological context. One aligns with an aeolian origin. The other (encompassing 97% of the depressions near the anticline) suggests a karst origin, with more circular shapes and proximity to karst, streams, faults and neighbouring depressions. However, 18% of these lie outside the crypto-karst area derived from GIS references and may rather be of periglacial origin. Yet, our geophysical survey revealed subsidence but no evidence of cryogenic sediment or solifluction. We conclude that a karst origin remains the best hypothesis in the anticline area, and highlight that sole reliance on GIS references may underestimate crypto-karst extent. Further geological investigation is therefore needed to map its full extent and related water pathways near the anticline.

## KEYWORDS

Aeolian processes (deflation), electrical resistivity tomography (ERT), ground-penetrating radar (GPR), karst processes (doline), periglacial processes (thermokarst)

## 1 | INTRODUCTION

In a context of increasingly frequent and long-lasting droughts, groundwater resources are of growing strategic importance. The analysis of landforms provides insights into flow paths and paleoclimate, both essential for inferring past and present groundwater recharge and guiding sustainable management strategies (Jasechko, 2016; Letz et al., 2021; Priestley et al., 2020).

In particular, small enclosed depressions from a few metres to a kilometre in diameter are ubiquitous features that can originate from a large variety of genetic processes, strongly determined by the lithological context and climate (Pardo-Igúzquiza & Dowd, 2021). Typically found in karst landscapes (Sauro, 2019), they are also widespread in glaciated environments, periglacial areas, or aeolian drylands (Grosse, Jones, & Arp, 2013; Hugenholtz & Wolfe, 2006; Prowse et al., 2006), spanning diverse hydrogeological settings. They may indicate complex

**Abbreviations:** BDLISA, Database of aquifer system boundaries; BSS, Subsoil databank of the French geological survey; BP, Before present (i.e., before 1950); BRGM, French geological survey; CBNSA, National botanical conservatory of South Aquitaine; CD40, Landes departmental council; DOI, Depth of investigation; CEN, Natural spaces conservatory; ERT, Electrical resistivity tomography; GPR, Ground-penetrating radar; HI, Hypsometric index; IGN, French national geographic institute; K<sub>G</sub>, Gravelius compactness coefficient; MPQ, Mio-Plio-Quaternary; PNRLG, Landes de Gascogne natural regional park; RMS(E), root mean square (error); SANDRE, French national service for water data and common repositories management; VL, Villagrains-Landiras.

This is an open access article under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2025 The Author(s). *Earth Surface Processes and Landforms* published by John Wiley & Sons Ltd.

<footnote n="1">¹EPOC, CNRS-UMR 5805, Bordeaux-INP, Pessac, France</footnote>

<footnote n="2">²Département de Géologie et de Génie Géologique, Université Laval, Québec City, Québec, Canada</footnote>
Highlights

Origins of thousands of small enclosed depressions in Southwest France: an integrated geomor-
phological, geostatistical and geophysical study

Léa Bussière, Myriam Schmutz, Alain Dupuy

* The Landes Triangle is dotted with 2,400+ shallow depressions of unknown origin.

* They fall into two distinct groups in terms of morphology, distribution and geology.

* Proportions of depressions vulnerable to karst processes are provided for each group.

* In central and western parts of the region, they are likely of aeolian origin.

* Elsewhere, crypto-karst origin is favored.
underground systems (Ford & Williams, 2007), function as focal points for groundwater recharge - including in permafrost areas (Rains, 2011) - and contribute to paleoclimatic reconstructions (Brosius et al., 2021).

However, due to morphologic similarities, identifying the geomor-
phic nature of shallow or relict enclosed depressions may prove chal-
lenging in the absence of clear contextual data (Goudie &
Wells, 1995; Majecka et al., 2019; Pissart, 2003). In some cases, even
the natural origin of such depressions is uncertain (Etienne
et al., 2011; Kołodyńska-Gawrysiak & Poesen, 2017). Misinterpreta-
tion may lead to the overlooking of possible hazards or the formula-
tion of inaccurate assumptions about paleoclimate. For instance, the
'viviers' on the Hautes Fagnes Plateau (Belgium), were long consid-
ered pingo scars from the Younger Dryas, implying past continuous
permafrost. However, Pissart (2003) later demonstrated that they are
instead thermokarst ponds associated with discontinuous permafrost.
Similarly, the 'spungs' of the Pine Barrens (New Jersey, USA) have
been successively interpreted as relic traces of discontinuous perma-
frost (thermokarst basins), continuous permafrost (pingo scars) and
seasonally frozen ground (deflation hollows) (French &
Demitroff, 2001). Since permafrost acts as an aquitard (Walvoord &
Kurylyk, 2016), such a distinction leads to contradictory reconstruc-
tions of groundwater recharge, affecting hydrogeological interpreta-
tions (Jiráková et al., 2011).

In this article, we focus on the origin of the ‘lagunes’ of the Landes Triangle, in southwest France. These shallow enclosed depressions are scattered across the Quaternary coversand of the Aquitaine sedimentary basin. Their highest density is found along the axis of the Villagrains-Landiras (VL) anticline (Legigan, 1979), a key recharge area for the regional, multilayered aquifer (Labat, 2021; Saltel, 2008), making them highly relevant to regional hydrogeology. To date, a total of 2,432 depressions have been inventoried (Figure 1), based on two morphological criteria: their shallow depth (approximately 2 m or less) and their circular or elliptical shape in aerial view. Many of them are filled with freshwater for at least part of the year and form ecologically valuable ponds. However, they exhibit notable diversity in sedimentary context, morphology and hydrological functioning.

Since the 1950s, three possible origins for these depressions have
been proposed. Initially, they were interpreted as periglacial features
resulting from permafrost thawing at the end of the Last Glacial
period (Legigan, 1979). Later, researchers suggested they were the
surface expression of karst processes affecting the underlying calcare-
ous formations (Becheler, 2014; Texier, 2011). Others have attributed
their formation to aeolian deflation (Corbier, 2015; Enjalbert, 1960;
Sitzia, 2014).

A periglacial origin would imply ground subsidence following either the melting of buried ice or the thawing of ice-rich permafrost. Indeed, in continuous permafrost, the melting of ice-cored hills known as pingos results in ramparted, enclosed depressions measuring 100 to 500 m across and a few metres deep (Ballantyne, 2018; Harris, Brouchkov, & Guodong, 2017). These features are generally scattered across the landscape, with densities of less than one per square kilometre (Pissart, 2002). Thermokarst lakes can also form through the melting of ice-wedge networks or through the broad, inhomogeneous subsidence of ice-rich ground (Grosse, Jones, & Arp, 2013). In discontinuous permafrost, the thawing of ice-rich permafrost mounds, sometimes referred to as lithalsas, produces similar ramparted depressions, known as thermokarst ponds. These range from 2 to 150 m in

diameter and are a few metres deep (Arlen-Pouliot & Bhiry, 2005; Seppälä, 2011). Unlike pingo scars, thermokarst ponds often occur in clusters ranging from a few to dozens of depressions.

The covered karst hypothesis would involve the dissolution or
collapse of underlying carbonate rocks or evaporites (Gökkaya
et al., 2021; Keskin & Yilmaz, 2016; Sauro, 2019). In such a context,
small enclosed depressions can correspond to buried sinkholes, suffu-
sion sinkhole or dropout sinkholes (Kranjc, 2013; Waltham &
Fookes, 2003). These landforms commonly have circular to subcircular
plan geometry, with a diameter from one to hundreds of metres and a
bowl- or funnel-shaped concave profile.

Lastly, aeolian depressions result from the deflation of lake sedi-
ments, sand ridge dunefields or fluvial drainage lines (Goudie &
Wells, 1995). In particular, blowouts are saucer-, cup- or trough-
shaped depressions formed by wind erosion on a preexisting sand
deposit (Hesp, 2002).

Previous studies on the 'lagunes' depressions were conducted by
Legigan & Marambat (1993), Texier (2011) and Faure & Galop (2011).
These researchers performed palynologic analysis and ¹⁴C dating on
shallow peat layers from the infilling of five depressions: three located
away from the VL anticline and two in its vicinity. The results revealed
a contrast between the two groups, with peat deposits in the former
dating to approximately 6,000 to 8,000 years before present (BP),
while those in the latter were only 2,000 to 3,000 years BP. Using
cross-sectional stratigraphic observation of four depressions (both
close and far from the anticline), Texier (2011) further identified sig-
nificant sedimentary deformations at their centres and margins. These
included mass flows, shears, ball-and-pillows and centimetric collapses
framed by short vertical strike-slip steps. They concluded that such
deformations could be attributed to both karst and periglacial pro-
cesses. However, Sitzia (2014) and Corbier (2015) later observed that
the location and elongated shape of depressions in the southern half
of the Landes Triangle are often influenced by dunes, suggesting an
aeolian origin. This raises the possibility that, despite being grouped
under a common term, the 'lagunes' may actually encompass distinct
geomorphic features.

In the study presented here, we combine GIS analysis with near-
surface geophysical investigations to specifically disentangle the pos-
sible karstic, aeolian and periglacial origins of these depressions and
assess their role in vertical water flows near the Villagrains-Landiras
anticline. After a short description of the study site, we compile
geospatial information to examine the morphology of all depressions
and their spatial distribution as a function of topography, lithology of
the surficial deposit, proximity to aeolian and karst features and prox-
imity to faults and folds. We partition the depressions into clusters
that share similar properties regarding these variables. Then, we ana-
lyse the vertical distribution of the karst and impervious layers, and
classify each depression with respect to its susceptibility to
crypto-karst processes (i.e., the development of karst beneath the
superficial layer of permeable sediment). Special attention is paid to
the 368 depressions located in the municipality of Saint-Magne, which
lies just above the axis of the VL anticline and hosts the highest den-
sity of depressions. In light of these results, we eventually select one
depression in a non-karst context near the VL anticline, and perform a
sedimentological and geophysical investigation to get insights into its
morphogenesis. This includes seven 2D electrical resistivity tomogra-
phies (ERT) and ground-penetrating radar (GPR) profiling. The findings
# Origins of thousands of small enclosed depressions in Southwest France: an integrated geomorphological, geostatistical and geophysical study

Léa Bussière<sup>a,b,*†</sup>, Myriam Schmutz<sup>a</sup> and Alain Dupuy<sup>a,‡</sup>

## ARTICLE INFO

Keywords:

GIS

Ground penetrating radar (GPR)

Electrical Resistivity Tomography (ERT)

Deflation

Doline

Thermokarst

## ABSTRACT

Small topographic depressions can provide information on past and present groundwater recharge processes. The Landes Triangle (SW France) is dotted with more than 2,400 of them, valued for their ecological significance. Yet, their origin and hydrogeological implications remain poorly understood. This paper aims to elucidate their formation process, with special emphasis on the Villagrains-Landiras (VL) anticline sector: a key recharge zone for regional aquifers, which also features the highest density of depressions. We first assess the likelihood of three historical hypotheses with open GIS data, then detail the near-surface structure of a depression in the VL sector with electrical resistivity tomographies and ground penetrating radar. We identify two distinct groups of depressions in terms of morphology, distribution and geological context (SE1 and SE2). SE2 shows characteristics consistent with aeolian origin. SE1 includes 97% of the VL anticline depressions, and exhibits circular morphology and proximity to karst features, natural streams, faults, and neighboring depressions, suggesting a karst origin. However, 18% of SE1 may rather be of periglacial origin, as they do not lie in the crypto-karst area derived from GIS references. Yet, the geophysical survey specifically designed to target periglacial features in this context showed only evidence of subsidence, and none of cryogenic sediment or solifluction. We conclude that a karst origin remains the best hypothesis in the VL anticline area, and highlight that sole reliance on GIS references may underestimate crypto-karst extent. Further geological investigation is therefore needed to fully capture its extent near the VL anticline, and map related water pathways.

## 1. Introduction

In a context of increasingly frequent and long-lasting droughts, groundwater resources are of growing strategic importance. The analysis of superficial landforms contributes to their sustainable management, by offering a means to identify preferential flow path and to gain information on current or past landscape dynamics that influence recharge processes. In particular, the observation of numerous shallow enclosed depressions in the quaternary coversand of the Landes Triangle (Aquitaine sedimentary basin, Southwest France) has triggered lively debates since the 1950s, concerning the nature of the sedimentary processes by which they appeared, and the information it provides on current and past hydrogeology. Indeed, as the coversand is the upper layer of the regional multilayered aquifer, any particular process affecting this formation could influence the recharge of the underlying aquifers and the quality of groundwater. The main concerns arise from the fact that the highest density of depressions is found along the axis of the

Villagrains-Landiras (VL) anticline (Legigan, 1979), one of the recharge areas of the multilayered aquifer (Labat, 2021; Saltel, 2008). Three hypotheses have emerged to explain the onset of these depressions: periglacial processes related to the thawing of permafrost during the Last Glacial (Legigan, 1979), karst processes affecting the underlying calcareous formations (Becheler, 2014, 2015; Texier, 2011), or eolian deflation (Enjalbert, 1960; Sitzia, 2014; Corbier, 2015). As they imply different assumptions about the subsurface processes and the paleoclimatic context of the region, it is essential to examine these hypotheses.

To date, 2,432 depressions have been inventoried (Figure 1), based on two morphological criteria: their shallowness (about 2 m depth or less) and their regular shape, either circular or elliptic in aerial view. Palynologic analysis and ¹⁴C dating were performed on shallow peat layers from the infilling of five of them: three located away from the VL anticline and two in its vicinity (Legigan and Marambat, 1993; Texier, 2011; Faure and Galop, 2011). The results showed a contrast between the two groups: the peats in the former were dated around 6,000 to 8,000 years before present (BP) while the others only 2,000 to 3,000 years BP. Using cross-sectional stratigraphic observation of four depressions (both close and far from the anticline), Texier (2011) further highlighted the presence of strong sedimentary deformations in the center and the border of the depressions. These include mass flows, shears, ball-and-pillows, and centimetric collapses framed by short vertical strike-slip

<figure><img src="image_3.png" /></figure>

<footnote><sup>a</sup>EPOC, CNRS-UMR 5805, Bordeaux-INP, 1 Allée F. Daguin, Pessac, 33607, France</footnote>

<footnote><sup>b</sup>Département de Géologie et de Génie Géologique, Université Laval, Québec City, Québec, Canada</footnote>

<footnote>*Corresponding author</footnote>

<footnote>lbussier@uni-heidelberg.de (L. Bussière)</footnote>

<footnote>ORCID(s): 0000-0002-2321-9449 (L. Bussière); 0000-0001-5363-7599 (M. Schmutz)</footnote>

<footnote n="1">¹Present address: Professorship in Hydrology and Climatology, Department of Geography, Heidelberg University, INF 348, Heidelberg 69120, Germany.</footnote>

<footnote n="2">²Present address: A. Dupuy, Program Director Groundwater and Global Change, BRGM-French Geological Survey, F-45060 Orléans, France.</footnote>
<figure><img src="image_2.png" /></figure>

Sources: Andrieux et al. (2016a) (periglacial features); BD Carthage (hydrology); BDCAVITE (cavities); BDLISA v.2 (coversand); BDLISA Zones karstiques (karst layers); Bourbon et al. (2019) (dunes); CBNSA, CD40, CEN Aquitaine and PNRLG (depressions); IGN (municipalities). Photo credits: Spadrone for Dt40.

<caption>**FIGURE 1** Distribution of the depressions in relation to the main karst, aeolian and periglacial features near and around the Landes Triangle. [Color figure can be viewed at wileyonlinelibrary.com]</caption>

provide new insights into the origin of the depressions in the Landes Triangle and especially near the VL anticline. Discussion of the results is followed by a comparison of our findings with previous studies, and by a reflection on the limitations of the available geospatial data and possible perspectives to overcome them.

## 2 | REGIONAL SETTING

The Landes Triangle is defined as the spatial extent of the Quaternary coversand in the Aquitaine basin. This aeolian deposit occupies a triangular surface of about 13,000 km², delimited by the Atlantic Ocean to the west, the Adour River to the south and the Garonne River to the north. Its thickness is greatest near the Atlantic coastline (50 to 60 m) and decreases down to a few metres towards the east. The main karst, aeolian and periglacial features found in the Landes Triangle are synthesized in Figure 1. No map of the karst rockhead exists for the region; however, Figure 2 synthesizes the stratigraphic position of the karst based on the “Karst zones” subset of the BDLISA (v2) - the reference database for aquifer system boundaries. Depth measurements to the corresponding karst layer are provided locally, according to the subsoil databank (BSS) of the French geological survey.

In terms of general geology, the deep substratum is composed of thick sedimentary layers accumulated over successive cycles of marine transgression and regression from the Triassic to the Miocene (Platel & Dubreuilh, 1986). The carbonate layers deposited during the marine phases were subsequently eroded and covered by clay-marl or even sand deposits during the regression phases, progressively building up a multi-layered aquifer. The geological structure exhibits pronounced asymmetry on either side of a NW-SE flexure that divides the Landes Triangle into two halves. To the north, sedimentary layers from the Mesozoic to the Neogene are between 500 and 2000 m thick (Platel & Dubreuilh, 1986) and are predominantly horizontally stratified, with occasional folds (Bourbon, 2019). In contrast, to the south, the sedimentary complex is significantly thicker, ranging from 4,000 to 6,000 m, and has been extensively folded as a result of Pyrenean orogenesis.

After the late Miocene, the basin evolved toward a continental context, with fluvial and aeolian sedimentation within a large delta (the future Landes Triangle), where sandy-gravelly and clayey layers were successively deposited. The Plio-Quaternary complex is structured into five sequences, each at least 10 m in thickness. However, their spatial distribution is patchy due to either erosion or a lack of deposition. As a consequence, few areas exhibit the entire lithological column. Schematically, each sequence starts with coarse sand or
<figure><img src="image_4.png" /><figcaption>Shallow depressions
in the Landes Triangle</figcaption></figure>

Sources: BD Carthage (hydrology); BDLISA v.2 (coversand); CBNSA, CD40, CEN Aquitaine and PNRLG (depressions); IGN (municipalities).

<caption>**Figure 1:** Synthetic map of the depressions generated from the aggregation of various recent inventories.</caption>

steps. They concluded that such deformations were both in agreement with karst and periglacial processes.

In the study presented here, we combine GIS analysis and near-surface geophysical investigation to help discriminating the process by which these depressions appeared, and to understand how they currently contribute to the vertical water flows near the Villagrains-Landiras anticline.

After a short description of the study site (section 2), we compile geospatial information to examine the morphology of all the depressions, and their spatial distribution as a function of: the topography, the lithology of the surficial deposit, the proximity to aeolian and karst features, and the proximity to faults and folds (sections 3.1 and 4.1). We partition the depressions into clusters that share similar properties regarding these variables (sections 3.2 and 4.2). Then, we analyze the vertical distribution of the karst and impervious layers, and classify each depression with respect to its vulnerability to crypto-karst processes (i.e., the development of karst beneath the superficial layer of permeable sediment) (sections 3.3 and 4.3). Special attention is paid to the 368 depressions located in the municipality of Saint-Magne, which lies just above the axis of the VL anticline and hosts the highest density of depressions. In light of these results, we eventually select one depression in a non-karst context near the VL anticline, and perform a sedimentological and geophysical investigation to get insights into its

morphogenesis (sections 3.4 and 4.4). This includes seven 2D electrical resistivity tomographies (ERT) and ground penetrating radar (GPR) lines. The findings provide new insights into the origin of the depressions in the Landes Triangle and especially near the VL anticline. They are discussed in section 5.1, followed by a confrontation of our results to previous studies (section 5.2), and by a critical reflection on the limitations of the available geospatial data and on possible perspectives to overcome them (section 5.3).

## 2. Regional setting

The spatial extent of the Land Triangle corresponds to that of the quaternary coversand in the Aquitaine basin. This aeolian deposit occupies a triangular surface of about 13,000 km², delimited by the Atlantic Ocean to the west, the Adour River to the south and the Garonne River to the north. The main karst, aeolian and periglacial features found there are synthesized in Figure 2.

Beneath the coversand, the substratum is composed of thick sedimentary layers accumulated during successive cycles of marine transgression and regression (Platel and Dubreuilh, 1986). The carbonate layers deposited during the marine phases were subsequently eroded and covered by clay-marl or even sand deposits during the regression phases, progressively building up a multi-layered aquifer. Numerous
<figure><img src="image_3.png" /></figure>

Sources: BD Carthage (hydrology); BDCAVITE (cavities); BDLISA v.2 (stratigraphy); BDLISA Zones karstiques (karst layers); BSS BRGM (well logs); IGN (municipalities).

<caption>**FIGURE 2** Stratigraphic position of the karst rockhead across the Landes Triangle, with verified depth and observed karst cavities. [Color figure can be viewed at wileyonlinelibrary.com]</caption>

gravel in a clayey matrix at its base, then gradually evolves to clayey sands and terminates with clay. The sequence ends with the aeolian sedimentation, which started from 400 ka BP. The parabolic and bar-chanoid dunes that currently overlain the coversand were formed after 13 ka BP (Bertran et al., 2020).

Together with this overlying coversand, the Mio-Plio-Quaternary (MPQ) complex constitutes the shallow aquifer, which reaches up to 100 m in thickness.

The overall stratigraphy is affected by a network of faults and folds oriented N110–120°E (Bourbon, 2019). These features were formed during two key tectonic events: extensional tectonics during the opening of the Bay of Biscay, and compressional tectonics during the Pyrenean orogeny (Platel & Dubreuilh, 1986).

Numerous caves and sinkholes are found where carbonate layers outcrop outside the Landes Triangle, at topographic heights (Figure 2). Inside the Landes Triangle, these karst units are found at depth, except along anticline folds, where they nearly outcrop (Cabaret et al., 2020; Labat et al., 2021; Perrin, Noury, & Cartannaz, 2014). Some natural caves and suffusion sinkholes have therefore been iden-tified in the eastern part of the Landes Triangle (Figure 2), where they are accompanied by sinking streams and springs. While most of them are a geological heritage, two new collapses appeared in 2011 and 2012 in the south-east.

According to global paleoclimatic models (Lindgren et al., 2016), the Landes Triangle lies outside the permafrost extent of the Last Glacial Maximum (26.5–19 ka BP). On the other hand, local proxies reported in Figure 1 (Andrieux et al., 2016; Andrieux, Bertran, & Saito, 2016) suggest that discontinuous permafrost may have extended down to this area during the Last Permafrost Maximum (25–17 ka BP), under favourable conditions (Bertran et al., 2014; Vandenberghe et al., 2014). However, the most recent publications reinterpret these proxies as the expression of an older permafrost episode or a deep seasonal freezing (Bertran et al., 2017). According to dune orientations and to grain-size gradients, dominant winds blew from the west and northwest during the Late Pleniglacial, then from the west and southwest during the Younger Dryas (Sitzia et al., 2017).

Today, the region is characterized by a temperate oceanic climate. The mean annual temperature is 14.2°C (over the period 1991–2020), and precipitation is relatively consistent throughout the year, averaging 925 mm annually. The wind regime changes seasonally, shifting from southeasterly winds in winter to westerly and northwesterly winds in summer (Sitzia et al., 2017). The mean wind speed decreases eastward and southward, ranging from 5 m.s⁻¹ at the northwestern vertex of the Landes Triangle to 2 m.s⁻¹ at the eastern vertex.
<figure><img src="image_5.png" /></figure>

Sources: Andrieux et al. (2016a) (periglacial features); BDLISA v.2 (coversand); BDLISA Zones karstiques (karst layers); BDCAVITE (caves); BD Carthage (hydrology); Bourbon et al. (2019) (dunes); IGN (municipalities).

<caption>**Figure 2:** Geological context: main karst, aeolian and periglacial features near and outside the Landes Triangle.</caption>

natural caves are found where the carbonate layers outcrop outside the Landes Triangle, at topographic heights. Inside the Landes Triangle, these karst units are found at depth, except along anticline folds, where they nearly outcrop (Cabaret et al., 2020; Labat et al., 2021; Perrin et al.).

After the late Miocene, the basin evolved toward a continental context, with fluvial and aeolian sedimentation within a large delta (the future Land Triangle), where sandy-gravelly and clayey layers were successively deposited. Together with the overlying coversand, this Mioplio-Quaternary (MPQ) complex constitutes the shallow aquifer, which reaches up to 100 m in thickness. The aeolian sedimentation started from 400 ka BP, but the parabolic and barchanoid dunes that currently overlain the coversand were formed after 13 ka BP (Bertran et al., 2020).

The overall stratigraphy is affected by a network of faults and folds inherited from extensional tectonics during the opening of the Bay of Biscay and compressional tectonics during the Pyrenean orogeny (Platel and Dubreuilh, 1986).

According to global paleoclimatic models (Lindgren et al., 2016), the Landes Triangle lies outside the permafrost extent of the Last Glacial Maximum (26.5-19 ka BP). On the other hand, local proxies reported in Figure 2 (Andrieux et al., 2016a,b) suggest that permafrost can have occurred down to this area during the Last Permafrost Maximum (25-17 ka BP), under favorable conditions (Vandenberghe et al., 2014; Bertran et al., 2014). However, the most recent

publications reinterpret these proxies as the expression of an older permafrost episode or a deep seasonal freezing (Bertran et al., 2017).

Today, the region is characterized by a temperate oceanic climate with 14.2°C mean annual temperature (over the period 1991-2020) and relatively continuous precipitations throughout the year, averaging 924.9 mm per year.

# 3. Materials and Methods

## 3.1. Selection of descriptive variables with principal component analysis

The 2,432 depressions were mapped by cross-checking four independent inventories carried out between 2005 and 2014, and updated on a case-by-case basis up to 2018 (Table 1). In case of inconsistency regarding the shape of the depressions between two inventories, the most recent mapping was retained.

For each depression, seven variables were measured:

* the lateral distances to the closest karst feature, dune, natural stream, fault or fold, and depression;

* the position of the depression on the hypsometric curve of its watershed (arbitrarily called “hypometric index” (HI) in this study);

* and the Gravelius compactness coefficient ($K_G$) of the depression.
# 3 | MATERIALS AND METHODS

## 3.1 | Selection of descriptive variables

The 2,432 depressions were mapped by cross-checking four independent inventories carried out between 2005 and 2014, and updated on a case-by-case basis up to 2018. In case of inconsistency regarding the delineation of the depressions between two inventories, the most recent mapping was retained.

For each depression, seven variables were measured. The first five are lateral distances, respectively from the depression to the nearest: (1) karst feature, (2) dune, (3) natural stream, (4) fault or fold and (5) neighbouring depression. Indeed, karst features and dunes are surface expressions of past karst and aeolian processes, respectively. Therefore, the proximity of a depression to either feature provides insights into the geomorphological context in which it formed. The closer a depression is to a karst or aeolian feature, the more likely its formation was influenced by the corresponding process. The proximity to streams and faults is known to favour karst processes (Keskin & Yılmaz, 2016; Perrin et al., 2015). Lastly, measuring the lateral distance to the nearest depression helps determine whether a given depression is part of a cluster or isolated.

These variables are complemented by two geomorphological criteria: (1) the position of the depression on the hypsometric curve of its watershed (referred to as the “hypometric index” or HI in this study), and (2) the Gravelius compactness coefficient ($K_G$) of the depression.

HI represents the percentage of the watershed's surface that lies at a lower elevation than the depression. It reflects the hydrological context of the depressions – a parameter that influences karstic, aeolian and karstic processes alike. The lower the index, the wetter the surface conditions.

$K_G$ is commonly used to describe the shape of watersheds (e.g., Sassolas-Serrayet, Cattin, & Ferry, 2018). It represents the ratio between the depression's perimeter and that of the circle of equal area:

$$K_G = P/(2 \sqrt{\text{PI}^* S}) \quad (1)$$

with $K_G$: Gravelius compactness coefficient (without dimension); P: depression's perimeter (m); and S: depression's area (m²). $K_G$ value is 1 for a perfect circle and increases as the shape becomes more elongated or irregular.

The geospatial data needed to compute these variables were obtained from governmental sources listed in Table 1.

The covariance table (Table 2) of these seven parameters was checked to evaluate the redundancy of information. The maximum correlation strength between two variables is -0.51, which remains quite low. This occurs between the distances to the nearest karst feature and dune, reflecting the geographical distribution of calcareous units around the border of the sandy unit (Figure 1). The other variables show no significant correlation with each other (all correlation coefficients are < 0.25 in absolute value). This means that each variable brings new and relevant information about the diversity of depressions.

<table><caption>TABLE 1 Data source for the calculation of the seven variables.</caption><thead><tr><th>Geospatial data</th><th>Database (producer, date)</th><th>Data source</th></tr></thead><tbody><tr><td rowspan="2">Karst features</td><td>BDLISA v.2 (BRGM, 2018)</td><td>http://sandre.eaufrance.fr</td></tr><tr><td>BDCavités (BRGM, 2021)</td><td>https://www.georisques.gouv.fr</td></tr><tr><td>Dunes</td><td>1/25,000th numerical geological map of Aquitaine (Bourbon, 2019)</td><td>https://sigesaqi.brgm.fr</td></tr><tr><td>Streams and watershed</td><td>BD Carthage v.3.0 (SANDRE, 2019)</td><td>http://services.sandre.eaufrance.fr</td></tr><tr><td>Faults and folds</td><td>1/25,000th numerical geological map of Aquitaine (BRGM, 2019)</td><td>https://sigesaqi.brgm.fr</td></tr><tr><td>DEM</td><td>RGE ALTI v.2.0 (IGN, 2018)</td><td>https://portail.pigma.org</td></tr></tbody></table>

## 3.2 | Clustering strategy with K-means

The clustering strategy was then based on the non-supervised K-means method (James et al., 2021; Mariani, Tweneboah, & Beccar-Varela, 2021). The depressions were first randomly assigned to a pre-specified number (K) of clusters. The centroids of these clusters were then computed based on the seven variables. Each depression was reassigned to the cluster with the nearest centroid, and this process was repeated iteratively until the clustering stabilized. It is important to note that the K-means approach has four limitations: (1) the algorithm always converges, regardless of the existence of a clear structure in the dataset, (2) the resulting clusters are sensitive to the presence of outliers, (3) the resulting clusters depend strongly on the first random step and (4) even when relevant, it says nothing about the underlying rationale for the clustering (James et al., 2021). To overcome these limitations, we implemented three additional processing steps.

First, to test the existence of a real K-structure, we cross-checked the results of 1,000 independent K-means clusterings and assigned a depression to a given subset (cluster) only if it was assigned to it more than 90% of the times. Below 90% affiliation, the depression's assignment was considered unknown. This first step has been tested on dummy datasets to evaluate its performance. Attempts to cluster a random dataset into K = 2 or K = 3 subsets resulted in more than 90% of the points being categorized as “unknown”. This strategy is thus robust and seems unable to cluster the data in absence of real structuration within the dataset. On the contrary, tests on dummy structured datasets also revealed that it may fail to detect the underlying clustering when the distinction among subsets is weak. In such cases, it results in a high proportion of unassigned data. This clustering strategy is consequently restricted to the detection of strong heterogeneities within the dataset.

Then, we filtered out the depressions with “unknown” assignment, ran again 1,000 independent K-means clusterings and selected the one with the best clustering score (i.e. with the lowest total within-cluster distances). This second step is necessary to eliminate the distortion effect produced by the filtered-out data (outliers) on the clusters' shapes. The purpose of running a large number of
**Table 1**

List of the four independent inventories used to map the depressions. The * mark indicates that the inventories have been regularly updated until 2018. CEN: Conservatoire d'Espaces Naturels (natural spaces conservatory); PNRLG: Parc Naturel Régional des Landes de Gascogne (Landes de Gascogne natural regional park); CBNSA: Conservatoire Botanique National Sud Aquitaine (national botanical conservatory of South Aquitaine); CD40: Conseil départemental des Landes (Landes departmental council).

<table><thead><tr><th>Author</th><th>Spatial extent</th><th>Date</th><th>Depressions count</th></tr></thead><tbody><tr><td>CEN</td><td>Gironde department</td><td>2008*</td><td>1051</td></tr><tr><td>PNRLG and CBNSA</td><td>Gironde department and PNRLG</td><td>2008-2014*</td><td>1273</td></tr><tr><td>CD40</td><td>Landes department</td><td>2005 - 2009*</td><td>462</td></tr><tr><td>Saltel</td><td>Saint-Magne</td><td>2008</td><td>393</td></tr></tbody></table>

**Table 2**

Data source for the calculation of the seven variables.

<table><thead><tr><th></th><th>Database (producer, date)</th><th>Data source</th></tr></thead><tbody><tr><td rowspan="2">Karst features</td><td>BDLISA v.2 (BRGM, 2018)</td><td>http://sandre.eaufrance.fr</td></tr><tr><td>BDCavités (BRGM, 2021)</td><td>https://www.georisques.gouv.fr</td></tr><tr><td>Dunes & Faults and folds</td><td>1/25,000th numerical geological map of Aquitaine (Bourbon, 2019)</td><td>https://sigesaqi.brgm.fr</td></tr><tr><td>Streams and watershed</td><td>BD Carthage v.3.0 (SANDRE, 2019)</td><td>http://services.sandre.eaufrance.fr</td></tr><tr><td>DEM</td><td>RGE ALTI v.2.0 (IGN, 2018)</td><td>https://portail.pigma.org</td></tr></tbody></table>

HI represents the percentage of the watershed's surface that lies at a lower elevation than the depression. $K_G$ represents the ratio between the depression's perimeter and that of the circle of equal area:

$$K_G = \frac{P}{2\sqrt{\pi S}} \qquad (1)$$

with $K_G$: Gravelius compactness coefficient (without dimension); $P$: depression's perimeter (m); and $S$: depression's area (m²).

The references used to calculate the variables are listed in Table 2.

We then tried to reduce the number of variables by carrying out a principal component analysis (PCA - see e.g. Mariani et al. (2021b); James et al. (2021)). This method consists in building new variables (or “descriptors”) by merging the seven original variables into linear combinations, in a way to limit repetition. The descriptors, although less explicit, are consequently more synthetic.

## 3.2. Clustering strategy with K-means

The clustering strategy was then based on the non-supervised K-means method (see e.g. Mariani et al., 2021a; James et al., 2021). In this method, the depressions are first randomly split into K unique and non-overlapping subsets in the p-dimension space of the (new) variables. Then, the K subsets are iteratively adjusted so that each one gather most similar depressions and is as far as possible from the other subsets. It is important to notice that the K-means approach has four limitations: (1) the algorithm always converges, regardless of the existence of a clear structure in the dataset,

(2) the resulting clusters are sensitive to the presence of outliers, (3) the resulting clusters depend strongly on the first random step, and (4) even when relevant, it says nothing about the underlying rationale for the clustering (James et al., 2021). The following three processing steps were therefore added to overcome these limitations.

First, to test the existence of a real K-structure, we cross-checked the results of 1,000 independent K-means clusterings, and assigned a depression to a given subset (cluster) only if it was assigned to it more than 90% of the times. Below 90% affiliation, the depression's assignment was considered unknown. This first step has been tested on dummy datasets to evaluate its performance. Attempts to cluster a random dataset into $K = 2$ or $K = 3$ subsets resulted in more than 90% of the points being categorized as “unknown”. This strategy is thus robust and seems unable to cluster the data in absence of real structuration within the dataset. On the contrary, tests on dummy structured datasets also revealed that it may fail to detect the underlying clustering when the distinction among subsets is weak. In such cases, it results in a high proportion of unassigned data. This clustering strategy is consequently restricted to the detection of strong heterogeneities within the dataset.

Then, we filtered out the depressions with “unknown” assignment, ran again 1,000 independent K-means clusterings and selected the one with the best clustering score (i.e. with the lowest total within-cluster distances). This second step is necessary to eliminate the distortion effect produced by the filtered-out data on the clusters shape. The purpose of running a large number of independent clusterings is to reduce the effect of the K-means random initialization on the result.
<table><caption>TABLE 2 Correlation table of the seven quantitative variables measured on the depressions. (FoF: faults or folds; Depres.: depressions; K<sub>G</sub>: Gravelius compactness coefficient; HI: hypsometric index).</caption><thead><tr><th></th><th>Karst</th><th>Dune</th><th>Stream</th><th>FoF</th><th>Depres.</th><th>K<sub>G</sub></th><th>HI</th></tr></thead><tbody><tr><th>Karst</th><td>1</td><td>-0.51</td><td>0.25</td><td>0.04</td><td>0.34</td><td>0.13</td><td>0.10</td></tr><tr><th>Dune</th><td></td><td>1</td><td>0.03</td><td>-0.03</td><td>-0.24</td><td>-0.19</td><td>0.14</td></tr><tr><th>Stream</th><td></td><td></td><td>1</td><td>0.21</td><td>0.17</td><td>-0.01</td><td>0.19</td></tr><tr><th>FoF</th><td></td><td></td><td></td><td>1</td><td>0.09</td><td>-0.01</td><td>0.08</td></tr><tr><th>Depres.</th><td></td><td></td><td></td><td></td><td>1</td><td>0.05</td><td>-0.00</td></tr><tr><th>K<sub>G</sub></th><td></td><td></td><td></td><td></td><td></td><td>1</td><td>-0.05</td></tr><tr><th>HI</th><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td></tr></tbody></table>

independent clusterings is to reduce the effect of the K-means random initialization on the result.

Eventually, appropriate statistical tests were performed to compare the clusters in pairs and get insights into the rationale for the clustering. For each of the seven variables, we tested the null hypothesis (HO): “there is no significant difference between the clusters regarding the variable i” versus the alternative one (H1): “there is a significant difference”. The HO hypothesis was rejected only when the p-value of the statistical test was lower than 0.001. The choice of the appropriate statistical test was based on the distribution of the variable of interest and the difference in its variance between clusters. These were respectively evaluated using a Shapiro-Wilk normality test and a Fischer-Snedecor test, with a significance level of 5%.

## 3.3 | Susceptibility assessment to karst processes through stratigraphic analysis

The susceptibility of each depression to karst processes was assessed following the methodological framework developed by Perrin, Noury, & Cartannaz (2014); Perrin et al. (2015), on the basis of a literature review (Doctor & Doctor, 2012; Farrant & Cooper, 2008; Galve et al., 2009; Lamelas et al., 2008; Orndorff, Weary, & Lagueux, 2000). This methodology postulates that undercover karstification is more efficient when water easily percolates down to the karst layers (promoting dissolution), and when dissolution products are efficiently evacuated. Specifically, five criteria are evaluated: (1) permeability of the sedimentary cover; (2) thickness of the sedimentary cover; (3) strength of water table fluctuation; (4) proximity to loosing streams; (5) proximity to geological discontinuities. The first two criteria assess water percolation efficiency, while criteria (3) to (5) evaluate the removal of dissolution products.

To simplify the approach, we built a decision tree solely on the first two criteria (Figure 3), while the last three (less decisive) were taken into account during the interpretation phase. The assessment was made on all the Landes Triangle depressions, but in a more detailed manner near the VL anticline (in Saint-Magne municipality), where a map of the thickness of the MPQ complex was provided by the French geological survey (Corbier and Karnay, 2010).

In both cases, impervious and potential karst layers (listed in Supporting information, Table S1) were identified according to the corresponding lists provided by the national aquifer reference system (BDLISA, version 2). An additional unit of shelly sand, sandstone and limestone (labelled “CodeEH: 320AA01” in BDLISA) was considered a

karst layer, as it shows a high susceptibility to karst near the Roquefort anticline (Perrin, Noury, & Cartannaz, 2014), and evidence of intense ghost-rock karstification along the VL anticline (Labat et al., 2021).

In France, Perrin, Noury, & Cartannaz (2014) have shown that covered karst generally occurs where karst layers are overlain by less than 40 m of permeable sediments. This conclusion is based on the examination of approximately 600 covered karst features across the territory. Since the BDLISA lacks information on layer thickness, we compiled 194 verified stratigraphic logs from the BSS database (French geological survey) across the Landes Triangle to estimate how the stratigraphic position of karst provided by the BDLISA relates to the actual depth of the karst layers. On average, the 40 m threshold is reached when the karst layer occupies the 5th position in the stratigraphy, with a mean depth of 48 ± 10 m (95% confidence interval). However, depth variations at this stratigraphic position are substantial, ranging from 9 to 150 m. Notably, the depth is significantly lower when the karst layer lies directly beneath the MPQ complex, averaging 30.5 m compared to 74.5 m in other cases. These observations were therefore included in the decision tree (Figure 3).

## 3.4 | Sedimentological and geophysical investigation of a depression in non-karst context

Shallow sediment sampling, ground-penetrating radar (GPR) profiling and electrical resistivity tomography (ERT) were carried out along 7 survey lines (L1 to L7) each 31.5 m in length, arranged as shown in Figure 4. The objective was to examine the near-surface sedimentary structure of the depression. Lines L1 and L7 crossed the depression's centre at two different angles, to reveal potential symmetries within its 3D structure. This information was essential for interpreting 2D ERT cross-sections. Lines L2 to L6, running parallel to L1, were designed to focus on the 3D structure of one half of the depression. To balance acquisition time with accuracy, the spacing between these lines increased progressively from the centre of the depression - where greater precision was required - towards its edges. Lastly, line L6 lay outside the depression (delimited by the amphibious grass vegetation belt shown in Figure 4) to serve as a reference for the general stratigraphy in areas without depressions. Absolute elevation of the ground surface was measured with a differential GNSS antenna every 50 cm along each line, with a maximum cumulative uncertainty of 3 cm in all three spatial dimensions. The antenna was connected to a mobile GNSS base installed a few metres away from the depression (44°32′46.37268" N, 0°41′31.86930" W,
Eventually, appropriate statistical tests were performed to compare the clusters in pairs. For each variable, we tested the null hypothesis ($H_0$): “there is no significant difference between the subsets regarding the variable $i$” versus the alternative one ($H_1$): “there is a significant difference”. The $H_0$ hypothesis was rejected only when the p-value of the statistical test was lower than 0.001. The choice of the appropriate statistical test was based on the statistical distribution of the variable of interest, and on the difference of its variance between the two subsets. These were respectively evaluated by a Shapiro-Wilk normality test and a Fischer-Snedecor test, with 5% significance.

### 3.3. Vulnerability assessment to karst processes through stratigraphic analysis

The vulnerability of each depression to karst processes was assessed following the national methodological framework developed by Perrin et al., on the basis of a literature review (Farrant and Cooper, 2008; Galve et al., 2009; Lame-las et al., 2008; Orndorff et al., 2000; Doctor and Doctor, 2012). This methodology postulates that undercover karstification is more efficient when water easily percolates down to the carbonate layers (promoting dissolution), and when dissolution products are efficiently evacuated. Specifically, five criteria are evaluated: water percolation efficiency is gauged by the permeability (1) and the thickness (2) of the sedimentary cover, whereas removal of dissolution products is gauged by the strength of water table fluctuation (3) and by the distances to loosing streams (4) and geological discontinuities (5).

To simplify the approach, we built a decision tree solely on the first two criteria (Figure 3), while the last three (less decisive) were kept in mind for the interpretation phase. The assessment was made on all the Landes Triangle depressions, but in a more detailed manner near the VL anticline (in Saint-Magne municipality), where more detailed information was available.

In both cases, impervious and potential karst layers (listed in supplementary material S1) were identified according to the corresponding lists provided by the national aquifer reference system (BDLISA, version 2). An additional unit of faluns, sandstone, and limestone (labelled “CodeEH: 320AA01” in BDLISA) was considered a karst layer, as it shows a high susceptibility to karst near the Roquefort anticline (Perrin et al.), and evidence of intense ghost-rock karstification along the VL anticline (Labat et al., 2021).

In the absence of impervious unit, Sinclair and Stewart (1985) noted that undercover karst efficiency drops when the thickness of the overlying layer exceeds 60 m. Since the BDLISA provides no information on layer thickness, we considered that this 60 m threshold was necessarily reached when the karst layer was buried under 5 layers. This approximation is deliberately large in order to overestimate rather than underestimate the vulnerability to karst. In Saint-Magne, the cover layer thickness was more precisely evaluated with the geological model of the MPQ complex

<figure><img src="image_6.png" /><figcaption>Figure 3: Decision tree used to assess the karst vulnerability of the Landes Triangle depressions, on the basis of the methodological framework proposed by Perrin et al., MPQ: Mio-Plio-Quaternary.</figcaption></figure>

produced by the French geological survey (Corbier and Karnay, 2010).

### 3.4. Sedimentological and geophysical investigation of a depression in non-karst context

Shallow sediment sampling, ground penetrating radar (GPR) profiling and Electrical Resistivity Tomography (ERT) were carried out along 7 survey lines (L1 to L7) of 31.5 m length, arranged as shown in Figure 4, to examine the near-surface sedimentary structure of the depression. Absolute elevation of the ground surface was measured with a differential GPS every 50 cm along each line, with a maximum cumulative uncertainty of 3 cm in all three spatial dimensions.

#### 3.4.1. Sampling and analysis of shallow sediments

Shallow sediment samples were collected in the depression at four different distances from its center (Figure 4). The samples were taken at 20 cm depth intervals, up to a maximum of 2.30 m. They have been oven-dried at 30°C during one week, then sieved through 14 meshes ranging from 4 mm to 45 µm. Grain size distribution was interpreted according to the Blott and Pye (2001) classification, and the Trask's sorting coefficient $So$ (Trask, 1930):

$$So = \sqrt{\frac{d_{25}}{d_{75}}} \qquad (2)$$

where $d_i$ is the $i^{th}$ percentile of the grain size distribution.
<figure><img src="image_4.png" /><figcaption>FIGURE 3 Decision tree used to assess the karst susceptibility of the Landes triangle depressions, on the basis of the methodological framework proposed by Perrin et al. (2015). MPQ: Mio-Plio-quaternary; *data only available in Saint-Magne municipality.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4 Location of survey lines L1 to L7 in the depression #135. Labels of lines L2 and L3 have been omitted for readability. Red dots represent the shallow sediment samples collected in distinct vegetation belts. [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

$\sigma = 0.0066 \text{ m})$. The real-time kinematic data were then post-processed using a combination of three reference GNSS bases, located in La Brède, Bordeaux and Lacanau.

### 3.4.1 | Sampling and analysis of shallow sediments

Shallow sediment samples were collected in the depression at four different distances from its centre (Figure 4). The samples were taken at

20 cm depth intervals, up to a maximum of 2.30 m. They were oven-dried at 30°C during one week to preserve the structure of the clayey fraction, then sieved through 14 meshes ranging from 4 mm to 45 µm. Grain size distribution was interpreted according to the Blott & Pye (2001) classification and Trask's sorting coefficient $So$ (Trask, 1930):

$$So = (d_{25}/d_{75})^{(1/2)} \quad (2)$$

where "$d_i$" is the $i^{th}$ percentile of the grain size distribution.

### 3.4.2 | Ground-penetrating radar

GPR profiling was achieved with a pulseEKKO PRO GPR (Sensors & Software) in October 2019, at two nominal frequencies: 100 and 250 MHz. The transmitter and receiver antennas were maintained a fixed distance apart and oriented perpendicularly to the survey lines. Only the 250 MHz antennas were shielded. Distinct survey parameters were used for each frequency. At 100 MHz, the antenna spacing was 1 m, and data were collected at a 0.25 m step size, with a 200 ns time window and 64 repeated measurements (stacks). Based on the average radar velocity in common soil materials (0.1 m.ns⁻¹), these settings allow for an estimated penetration depth of 10 m. At 250 MHz, the antennas were mounted on a SmartCart (Sensors & Software), with a spacing of 0.38 m. Data were collected at a 0.05 m step size with a 100 ns time-window, enabling a penetration depth of approximately 5 m under the same velocity assumption. The number of stacks at this frequency was automatically adjusted by the GPR system as a function of the travel speed.

Data were processed with the software Reflexw (Sandmeier geo-physical research), starting with time-zero correction, followed by
<figure><img src="image_7.png" /><figcaption>**Figure 4:** Location of survey lines L1 to L7 in the depression #135. Red dots represent the shallow sediment samples, collected in distinct vegetation belts.</figcaption></figure>

### 3.4.2. Ground penetrating radar

GPR profiling was achieved with a pulseEKKO PRO GPR (Sensors & Software) in October 2019, at 2 nominal frequencies: 100 and 250 MHz. The transmitter and receiver antennas were maintained a fixed distance apart and oriented perpendicularly to the survey lines. Depending on the frequency, Table 3 details the values for (1) the antenna separation, (2) the distance interval between two consecutive radar traces along one line, (3) the time window, and (4) the number of repeated measurements (stacks) to produce one trace. At 250 MHz, the antennas were mounted on a SmartCart (Sensors & Software), so the traces were automatically acquired and the stacks adjusted by DynaQ method as a function of the travel speed.

<caption>**Table 3**</caption>

GPR survey parameters as a function of nominal frequency ($d_{Rx-Tx}$: antenna separation; TW: time window).

<table><thead><tr><td>Nominal frequency (MHz)</td><td>d<sub>Rx-Tx</sub> (m)</td><td>Step size (m)</td><td>TW (ns)</td><td>Number of stacks</td></tr></thead><tbody><tr><td>100</td><td>1</td><td>0.25</td><td>200</td><td>64</td></tr><tr><td>250</td><td>0.38</td><td>0.05</td><td>100</td><td>DynaQ</td></tr></tbody></table>

Data were processed with the software Reflexw (Sandmeier geophysical research) to perform successively time-zero correction, topographic correction and signal amplification by manual gain. Radar velocity was estimated by diffraction adjustment fitting in sandy areas, and by collecting GPR data immediately above a buried metallic electrode (at known depth) in the peaty central part of the depression.

### 3.4.3. GPR interpretation strategy

The output of GPR processing is a lateral juxtaposition of all traces measured along a given profile, called

radargram. The consistent measurement of a peak in radar intensity from one trace to another draws a reflector, which represents a contrast in dielectric permittivity. Following Annan (2005) and Everett (2013), we interpreted such a contrast as a change in water content or lithology. Whenever the change is accompanied by an attenuation of the radar signal, we considered that the electrical conductivity of the medium increases (e.g., due to an increase in water content or a decrease in grain size). Finally, in the absence of attenuation, the absence of reflector is assumed to reflect a homogeneous medium.

### 3.4.4. Electrical resistivity tomographies

ERTs were carried out separately along each line in September 2019, in dipole-dipole array, by means of sixty-four metal electrodes grounded at regular 50 cm intervals. We used a DAS-1 impedance meter (Multi-Phase Technologies, LLC) with 1 Hz alternating current. The dipole-dipole sequence involved 1,653 quadrupoles with a minimum electrode spacing of 1 m. Two reciprocal measurements were systematically performed on each quadrupole by swapping the injection electrodes. We also ensured that the contact resistance of the electrodes was lower than 3,000 ohm, and that a minimum of 2.5 mA was received at the electrodes. Data with more than 1% standard error between reciprocal measurements were filtered out, as were those with effective electrode currents below 5 mA.

ERT inversions were performed with the software RES-2DINV (Loke and Barker, 1995, 1996) with robust data constrain and robust model constrain, as well as topographic correction. This software solves the inverse problem (see e.g. Binley, 2015; Binley and Kemna, 2005) by means of an optimization function based on the incomplete Gauss-Newton least-squares equation. The subsurface was discretized in finite elements with a width of 25 cm. The convergence criterion was a decrease of less than 5% in root mean square error (RMSE) between two consecutive iterations.

Reliability of the resulting models of electrical resistivity was evaluated by the depth of investigation (DOI) method proposed by Oldenburg and Li (1999). Parts of the models with DOI index values above 0.1 were blanked to avoid misinterpretation.

### 3.4.5. ERT interpretation strategy

We assumed that the salinity of the porous medium does not vary on the scale of the survey lines, and that therefore any change in resistivity is attributable to a change in water content, lithology, or both. Especially, a decrease in water content (other things being equal) would result in an increase in resistivity, as would an increase in grain size (and vice versa).

## 4. Results

### 4.1. Principal component analysis

According to the covariance table (Table 4), the maximum correlation strength between two variables is -0.51, which remains quite low. It is obtained between the distance
signal amplification with manual gain (applied exclusively to the data collected at 250 MHz), and concluding with topographic correction. Radar wave velocity was estimated by diffraction adjustment fitting in sandy areas, and by collecting GPR data immediately above a buried metallic electrode (at 60 cm depth) in the peaty, central part of the depression.

### 3.4.3 | Electrical resistivity tomographies

ERTs were carried out individually along each line in September 2019, in Wenner array, by means of sixty-four metal electrodes grounded at regular 50 cm intervals. The aim was to delineate the geometry of sedimentary deposits in the depression and its surrounding area, with a focus on identifying potential traces of solifluction or cryogenic sediments at depth. These targets typically produce lateral variations in the subsurface stratigraphy, which in turn affect electrical resistivity. The Wenner array, characterized by good vertical sensitivity, is well suited for detecting such variation (Binley, 2015). Moreover, the proximity of the injection and potential electrodes ensures a strong signal.

We used a DAS-1 impedance metre (Multi-Phase Technologies, LLC) with 1 Hz alternating current. The Wenner sequence involved 144 quadrupoles with a minimum electrode spacing of 1 m. Two reciprocal measurements were systematically performed on each quadrupole by swapping the injection electrodes. We also ensured that the contact resistance of the electrodes was lower than 3,000 Ω, and that a minimum of 2.5 mA was received at the electrodes. Data with more than 1% standard error between reciprocal measurements were filtered out, as were those with effective electrode currents below 5 mA.

ERT inversions were performed with the software RES2DINV (Loke & Barker, 1995, 1996) with robust data constrain and robust model constrain, as well as topographic correction. This software solves the inverse problem (Binley, 2015; Binley & Kemna, 2005) by means of an optimization function based on the incomplete Gauss-Newton least-squares equation. The subsurface was discretized in finite elements with a width of 50 cm and a length that increases with depth, varying from 0.25 m to 4.9 m. The convergence criterion was a decrease of less than 5% in root mean square error (RMSE) between two consecutive iterations.

The reliability of the resulting electrical resistivity models was assessed using the depth of investigation (DOI) method proposed by Oldenburg & Li (1999), which helps identify regions of the models that are weakly constrained by the data. This method involved performing two separate inversions on the same dataset, each starting from a distinct homogeneous initial model: one with a resistivity value 10 times lower than the average observed apparent resistivity, and the other 10 times higher. The DOI index R quantifies the discrepancy in resistivity between the two inversions. It was calculated for each cell using Equation 3, where $q_{m1}$ and $q_{m2}$ are the electrical resistivities of the two initial models, and $q_1$ and $q_2$ are the cell's electrical resistivities obtained from the corresponding inversions.

$$R = (q1 - q2)/(qm1 - qm2) \quad (3)$$

A high DOI index indicates that the computed value is strongly influenced by the initial model rather than being reliably determined

by the data. To prevent misinterpretation, model regions with DOI index above 0.1 were blanked.

### 3.4.4 | ERT interpretation strategy

We assumed that salinity in the porous medium remained constant at the scale of the survey lines. Therefore, any resistivity changes (ρ, in ohm.m) were attributed to variations in water content, lithology, or both. Specifically, a decrease in water content is known to increase resistivity, as does an increase in grain size (and vice versa).

## 4 | RESULTS

### 4.1 | Clustering

The results of the clustering strategy (based on cross-checking 1,000 independent K-means runs) on the depressions dataset is shown in Table 3 for cases K = 3 and K = 2. For K = 3, 76% of the depressions were not successfully assigned to a cluster, while the remaining 24% were all assigned to the same one. The attempt to group the depressions into three strongly different homogeneous subsets was thus unsuccessful. In contrast, with K = 2, the clustering consistently assigned 96% of the depressions to a specific subset. The first subset (SE1) included 1,403 depressions, while the second (SE2) included 939. The 2-subsets clustering was therefore validated.

After filtering out the 90 unassigned depressions, 1,366 depressions belonged to SE1 and 976 to SE2. Their spatial distribution is shown in Figure 5, and their statistical comparison is presented below.

Finally, according to the Shapiro-Wilk test, none of the seven variables followed a normal distribution. Therefore, the appropriate statistical test for comparisons was the Wilcoxon test. This analysis revealed significant contrasts between the SE1 and SE2 subsets for each variable. In particular, depressions in SE1 are significantly closer to karst features, natural streams, faults or folds and other depressions, while they are farther from dunes. They are also situated at lower elevations within their watersheds and exhibit a more circular shape compared to those in SE2.

### 4.2 | Susceptibility to karst processes

There is a strong dissymmetry between the SE1 and SE2 subsets regarding susceptibility to karst processes. Although 55% of the

<table><caption>TABLE 3 Output of the clustering strategy on the Landes Triangle depressions, after cross-checking 1,000 independent runs of K-means algorithm with K = 3 clusters (right) or K = 2 clusters (left). (prop.: proportion; Unkn.: unknown).</caption><thead><tr><th colspan="3">(a) Clustering with K = 3</th><th colspan="3">(b) Clustering with K = 2</th></tr><tr><th>Cluster</th><th>Count</th><th>Prop.</th><th>Cluster</th><th>Count</th><th>Prop.</th></tr></thead><tbody><tr><td>SE1</td><td>572</td><td>24%</td><td>SE1</td><td>1,403</td><td>58%</td></tr><tr><td>SE2</td><td>0</td><td>0%</td><td>SE2</td><td>939</td><td>39%</td></tr><tr><td>SE3</td><td>2</td><td>0%</td><td>Unkn.</td><td>90</td><td>4%</td></tr><tr><td>Unkn.</td><td>1858</td><td>76%</td><td></td><td></td><td></td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Figure 5:** Scree plot of the principal component analysis, showing the proportion of variance explained by the new descriptors in the depressions dataset (left: original, right: cumulative).</figcaption></figure>

to closest karst feature and the distance to the closest dune, and reflects the geographical distribution of the calcareous units around the border of the sandy unit (Figure 2). The other variables are not significantly correlated with each other (all the correlation coefficients are less than 0.25 in absolute value). This means that each variable brings new and relevant information about the diversity of depressions. For this very reason, the PCA was unable to come up with a limited number of descriptors that could reflect their diversity: to recover 80% of the observed variance, one must keep at least the first five new descriptors issued from the PCA (Figure 5). This reduction in the number of descriptors (from 7 to 5) is not sufficient enough to compensate the related loss of information. Therefore, the following work (including the clustering of depressions) will continue to rely on the seven original variables.

<caption>**Table 4**</caption>

Correlation table of the seven quantitative variables measured on the depressions (FoF: Faults or folds; Depres.: Depressions).

<table><thead><tr><th></th><th>Karst</th><th>Dune</th><th>Stream</th><th>FoF</th><th>Depres.</th><th>K<sub>G</sub></th><th>HI</th></tr></thead><tbody><tr><th>Karst</th><td>1</td><td>-0.51</td><td>0.25</td><td>0.04</td><td>0.34</td><td>0.13</td><td>0.10</td></tr><tr><th>Dune</th><td></td><td>1</td><td>0.03</td><td>-0.03</td><td>-0.24</td><td>-0.19</td><td>0.14</td></tr><tr><th>Stream</th><td></td><td></td><td>1</td><td>0.21</td><td>0.17</td><td>-0.01</td><td>0.19</td></tr><tr><th>FoF</th><td></td><td></td><td></td><td>1</td><td>0.09</td><td>-0.01</td><td>0.08</td></tr><tr><th>Depres.</th><td></td><td></td><td></td><td></td><td>1</td><td>0.05</td><td>-0.00</td></tr><tr><th>K<sub>G</sub></th><td></td><td></td><td></td><td></td><td></td><td>1</td><td>-0.05</td></tr><tr><th>HI</th><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td></tr></tbody></table>

## 4.2. Clustering

The result of the clustering strategy (by cross-checking 1,000 independent K-means run) on the depressions dataset is shown in Table 5 for cases $K = 3$ and $K = 2$. For $K = 3$, 76% of the depressions are not successfully assigned to a cluster, while the remaining 24% are all assigned to the same one. The attempt to group the depressions into three strongly different homogeneous subsets is thus a failure. On

the contrary, with $K = 2$, the clustering consistently assigns 94% of the depressions to a given subset: 1,403 to the SE1 one and 939 to the SE2 one. The 2-subsets clustering is therefore validated.

After filtering out the 90 unassigned depressions, 1,366 depressions belong to SE1 and 976 to SE2. Their spatial distribution is shown in Figure 6, and their statistical comparison is presented below.

<caption>**Table 5**</caption>

Output of the clustering strategy on the Landes Triangle depressions, after cross-checking 1,000 independent runs of K-means algorithm with $K = 3$ clusters (right) or $K = 2$ clusters (left) (Prop.: Proportion; Unkn.: Unknown).

<table><thead><tr><th>Cluster</th><th>Count</th><th>Prop.</th></tr></thead><tbody><tr><td>SE1</td><td>572</td><td>24%</td></tr><tr><td>SE2</td><td>0</td><td>0%</td></tr><tr><td>SE3</td><td>2</td><td>0%</td></tr><tr><td>Unkn.</td><td>1858</td><td>76%</td></tr></tbody></table>

(a) Clustering with $K = 3$.

<table><thead><tr><th>Cluster</th><th>Count</th><th>Prop.</th></tr></thead><tbody><tr><td>SE1</td><td>1403</td><td>58%</td></tr><tr><td>SE2</td><td>939</td><td>39%</td></tr><tr><td>Unkn.</td><td>90</td><td>4%</td></tr></tbody></table>

(b) Clustering with $K = 2$.

Finally, according to the Shapiro-Wilk test, none of the 7 variables follows a normal distribution. The appropriate comparison test is thus always the Wilcoxon's test. This revealed that SE1 and SE2 subsets show significant contrasts for each variable. Especially, the SE1's depressions are significantly:

* closer to karst features, natural streams, faults or folds, and other depressions;

* farther from dunes;

* located at a lower elevation within their watershed;

* and more circular than the SE2's ones.

## 4.3. Vulnerability to karst processes

There is a strong dissymmetry between the SE1 and SE2 subsets regarding vulnerability to karst processes. Although 45% of the Landes Triangle depressions lie in a context favorable to karstification, this proportion rises to 82% in SE1 subset, and drops to only 18% in the SE2 one (Figure 7).

In Saint-Magne, most of the depressions belong to SE1 subset (358), compared with just 5 in SE2 and 5 unassigned. According to stratigraphic data (layers permeability, thickness, and susceptibility to karst), 66% of them lie in a context suitable for crypto-karst processes, with a “high” to “very high” degree of karst vulnerability (Figure 8). However, the remaining 34% are located in a context that is not favorable to crypto-karst, due to the absence of potential karst layer within the five first geological layers (6 depressions), or to the presence of an impervious layer between the MPQ complex and the deep potential karst layer (119).

## 4.4. Investigation of depression #135

The depression #135 - selected for further geophysical investigation - is located in a non-karst context along the VL anticline, and belongs to the SE1 subset. This is a circular depression 19.6 m in diameter and 1.85 m deep. Geological
Clustering of
the shallow depressions

<figure><img src="image_6.png" /></figure>

Sources: BD Carthage (hydrology); CBNSA, CD40, CEN
Aquitaine and PNRLG (depressions); IGN (municipalities).

<caption>**FIGURE 5** Partition of the depressions into two homogeneous clusters (n: depressions count). [Color figure can be viewed at wileyonlinelibrary.com]</caption>

Landes Triangle depressions were found to lie in a context favourable to karstification, this proportion rises to 82% in SE1 subset and drops to only 18% in the SE2 one (Figure 6).

In Saint-Magne, most of the depressions were found to belong to SE1 subset (358), compared with just 5 in SE2 and 5 unassigned. According to stratigraphic data (layers permeability, thickness and susceptibility to karst), 66% of them lied in a context suitable for crypto-karst processes with a “high” to “very high” degree of karst susceptibility (Figure 7). However, the remaining 34% were located in a context that was not favourable to crypto-karst, due to the absence of a potential karst layer within the five first geological layers (6 depressions), or to the presence of an impervious layer between the MPQ complex and the deep potential karst layer (119).

**4.3** | **Investigation of depression #135**

Depression #135 (Figure 7) is a circular depression 19.6 m in diameter and 1.85 m deep. It was selected for further geophysical investigation because it was found to lie in a non-karst context along the VL anti-cline, while belonging to the SE1 subset. In addition, geological soundings and piezometric monitoring of the shallow MPQ aquifer were

available within 300 to 500 m, providing valuable benchmarks for the hydrogeological context. Finally, the surrounding area was free of electromagnetic noise or metallic structures that could alter the measured signal.

4.3.1 | Shallow sedimentary structure

The lithological information issued from the analysis of shallow samples collected near the line L1 is synthetized in Figure 8.

In its centre, the depression was partially filled with highly mineralized peat, whereas the host sediment in its vicinity was a medium (250–500 µm) to coarse sand (500 µm - 1 mm). The sand grains were moderately sorted (So in 1.35–1.87), except just beneath the peat, where they were poorly sorted (So in 1.87–2.75).

4.3.2 | Ground-penetrating radar

All radargrams are provided in Supporting information (Figure S1, and a selection is shown in Figure 9 (top and middle). Radar speed evaluation was 0.06 m.ns⁻¹ within the peat and 0.11 m.ns⁻¹ within the
<figure><img src="image_9.png" /></figure>

Sources: BD Carthage (hydrology); CBNSA, CD40, CEN
Aquitaine and PNRLG (depressions); IGN (municipalities).

<figure><img src="image_10.png" /><figcaption>**Figure 7:** Number (and proportion) of depressions vulnerable to karst processes in each subset.</figcaption></figure>

sounding and piezometric monitoring of the shallow MPQ aquifer are available within 300 to 500 m and provide

benchmarks for the hydrogeological context. In addition, the surrounding of the depression is free of electromagnetic noise or metallic structures that could alter the measured signal.

### 4.4.1. Shallow sedimentary structure

The lithological information issued from the analysis of shallow samples collected near the line L1 is synthetized in Figure 9.

In its center, the depression is partially filled with highly mineralized peat, whereas the host sediment in its vicinity is a medium (250-500 µm) to coarse sand (500 µm - 1 mm). The sand grains are moderately sorted (So in 1.35-1.87), except just beneath the peat, where they are poorly sorted (So in 1.87-2.75).

### 4.4.2. Ground penetrating radar

All radargrams are provided in Supplementary material S2 and a selection is shown in Figure 10 (top and middle). Radar speed evaluation is 0.06 m.ns⁻¹ within the peat and 0.11 m.ns⁻¹ within the sand. In Figure 10 and S2, the speed 0.06 m.ns⁻¹ was taken as a reference to convert the vertical axis of the radargrams from two-ways travel time to depth. The depth of investigation drops from 3 m at 100 MHz, to 2 m at 250 MHz. No signal amplification
sand. In Figure 9 and S1, the speed 0.06 m.ns⁻¹ was taken as a reference to convert the vertical axis of the radargrams from two-ways travel time to depth. The depth of investigation dropped from 3 m at

100 MHz, to 2 m at 250 MHz. Regardless of frequency, the position and shape of the reflectors on the radargrams from lines L1 to L4 and L7 were very similar.

<figure><img src="image_7.png" /></figure>

**FIGURE 6** Number (and proportion) of depressions vulnerable to karst processes in each subset. [Color figure can be viewed at wileyonlinelibrary.com]

### 4.3.3 | Tomographies of electrical resistivity

All inverse models of electrical resistivity are provided in Supporting information (Figure S1), and a selection is shown in Figure 9 (bottom). Data filtering removed three points from line L1, one from L4, eleven from L5, six from L6 and none from the other lines. On every model, the RMSE was lower than 4.6%, and the DOI was about 3 m. The interpretation is given together with that of the GPR results in the following paragraph.

### 4.3.4 | Interpretation of the GPR and ERT images

The strongest GPR reflector on the 100 MHz radargrams of line L1 to L4 (and L7) highlights a gently curved contact between the depression's infill and the substrate. This contact fits precisely with both the base of the conductive unit (ρ < 600 Ω.m) on the electrical resistivity models and the boundary between the peat and sand layers observed during sediment sampling (see the dashed line on Figure 9). Therefore, GPR and ERT methods efficiently allow to map the 3D structure of the peat layer.

<figure><img src="image_8.png" /></figure>

Sources: Andrieux et al. (2016) (periglacial features); BDLISA Zones karstiques (karst layers); Bourbon et al. (2019) (dunes); CBNSA, CEN Aquitaine and PNRLG (depressions); IGN (municipalities); Saltel (2008) (anticline axis).

**FIGURE 7** Degree of susceptibility of Saint-Magne depressions to karst processes (n: depressions count). [Color figure can be viewed at wileyonlinelibrary.com]
<figure><img src="image_11.png" /></figure>

Sources: Andrieux et al. (2016) (periglacial features); BDLISA Zones karstiques (karst layers); Bourbon et al. (2019) (dunes); CBNSA, CEN Aquitaine and PNRLG (depressions); IGN (municipalities); Saltel (2008) (anticline axis).

<figure><img src="image_12.png" /><figcaption>**Figure 9:** Stratigraphic interpretation of the shallow samples collected across depression #135, along with the stratigraphic column extracted from the geological model of the MPQ complex (with layers names).</figcaption></figure>

was needed at 100 MHz. Regardless of the frequency, the radargrams of lines L1 to L4 and L7 are very similar.

### 4.4.3. Tomographies of electrical resistivity

All inverse models of electrical resistivity are provided in Supplementary material S2 and a selection is shown in Figure 10 (bottom). One hundred seventy-nine data points were filtered out along the survey line L1, but none on the other lines. On every model, the RMSE is lower than 4.6%

and the DOI is about 3 m. The interpretation is given together with that of the GPR results in the following paragraph.

### 4.4.4. Interpretation of the GPR and ERT images

The main (thickest) GPR reflector on the 100 MHz radargrams of line L1 to L4 (and L7) highlights a gently curved contact between the depression's infill and the host formation. This contact fits precisely with both the base of the conductive unit ($<$ 600 ohm.m) on the electrical resistivity models, and the boundary between the peat and sand layers observed during sediment sampling. Therefore, GPR and ERT methods efficiently allow to map the 3D structure of the peat layer.

The top of peat has a blurry aspect on the 100 MHz radargrams, and displays multiple identical reflectors at 250 MHz. These are GPR artifacts, caused by the bouncing of radar waves between the antenna and the peat, due to the high contrast in dielectric permittivity $ε$ between air ($ε$ = 1) and wet peat (60 ≤ $ε$ ≤ 80). However, the bottom part of the peat is clearly visible at 250 MHz and depicts a homogeneous unit, as shown by the absence of minor reflectors.

On the other hand, the host formation is highly resistive ($>$ 5,000 ohm.m) near the surface and shows many minor reflectors on the 250 MHz radargrams. These correspond to boundary reflections in the layered sand. The sharp decrease in electrical resistivity with depth (down to 500-1,000 ohm.m) matches the position of the water table
<figure><img src="image_9.png" /><figcaption>**FIGURE 8** Stratigraphic interpretation of the shallow samples collected across depression #135, along with the stratigraphic column extracted from the geological model of the MPQ complex (with layers names). [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

On the 100 MHz radargram of survey line L1, a strong reflector divides the peat unit into two layers. The upper layer exhibits very low resistivity (ρ < 200 Ω.m), while the lower one is slightly more resistive (200 < ρ < 600 Ω.m). This contrast marks the transition from saturated to dry peat at depth, caused by the downward propagation of a wetting front following a rain event. At 250 MHz, the bottom part of the peat is clearly visible and depicts a homogeneous unit, as shown by the absence of weak reflectors.

On the other hand, the substrate is highly resistive (ρ > 5,000 Ω. m) near the surface and shows many weak reflectors on the 250 MHz radargrams. These correspond to boundary reflections in the layered sand. The sharp decrease in electrical resistivity with depth (down to 500-1,000 Ω.m) matches the position of the water table measured in the well PQ-08515X0053 (250 m away) at the time of the ERT survey. It represents the transition to the saturated zone. More interestingly, while the substrate weak reflectors are subparallel and

<figure><img src="image_10.png" /><figcaption>**FIGURE 9** Results of the geophysical investigation along the survey lines L1 (left) and L4 (right). From top to bottom: (1) GPR at 250 MHz, (2) GPR at 100 MHz and (3) ERT models. All the graphs are at the same scale (assuming a wave speed of 0.06 m.ns⁻¹ for the radargrams), with vertical exaggeration by a factor of two. Black arrows locate the intersection with the survey line L7. The displayed water tables are those measured specifically on the day of the corresponding survey in well PQ-08515X0053 (250 m away). [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 10:** Results of the geophysical investigation along the survey lines L1 (left) and L4 (right). From top to bottom: (1) GPR at 250 MHz, (2) GPR at 100 MHz, and (3) ERT models. All the graphs are at the same scale (assuming a wave speed of 0.06 m.ns⁻¹ for the radargrams), with vertical exaggeration by a factor of two. The displayed water tables are those measured specifically on the day of the corresponding survey in well PQ-08515X0053 (250 m away).</figcaption></figure>

measured in the well PQ-08515X0053 (250 m away) at the time of ERT survey. It reflects the transition to saturated zone. Likewise, the lateral increases in electrical resistivity to very high values (> 10,000 ohm.m) near the surface on L4, L5 and L6 relate to moisture variations rather than to lithological ones, because they are not correlated with any change on the radargrams. More interestingly, while the host minor reflectors are subparallel and horizontal away from the depression (as visible at 250 MHz on lines L5 and L6), they tend to converge below the pond, before being truncated by the major reflector. Figure 11 illustrates this phenomenon along L1. This reveals that ground subsidence occurred below the pond, either before or during the peat deposition. Despite most of the minor reflectors depicted in the sand unit just below the main reflector at 100 MHz are artifacts, some of them are genuine and indicate that this subsidence has spread a little deeper (Supplementary material S3). In addition, the step-like shape of the main reflector on either side of the peat infill suggests that successive small collapses occurred during the peat deposition.

# 5. Discussion

## 5.1. Updated hypothesis on the origin of the depressions

Despite three morphogenetic processes could explain their onset, only two types of depressions have been observed. Those in the SE2 group have irregular shape and proximity to dunes, which are consistent with the aeolian origin by deflation proposed by Sitzia (2014). Blowout hollows indeed develop by wind erosion in dunes (Hesp, 2002) and are typically described as semi-circular or elongated shallow depressions, respectively called “saucer-shaped blowouts” or “trough blowouts” (Cooper, 1958, 1967). In the Landes Triangle, they may result from the erosion of low amplitude dunes (0.3-0.7 m high) of Pleniglacial age (62-42 ka BP) in the northern half, or of more recent parabolic dunes (12.9-11.7 ka BP) in the southern half (Sitzia et al., 2017). Interestingly, SE2 depressions are evenly distributed between these north and south halves (49% vs 51%, respectively).
horizontal away from the depression (as visible at 250 MHz on lines L5 and L6), they tend to converge below the pond, before being truncated by the strong reflector. Figure 10 illustrates this phenomenon along L1. This reveals that ground subsidence occurred below the pond, either before or during the peat deposition. Most of the weak reflectors depicted in the sand unit at 100 MHz, just below the strongest reflector, are ‘multiples’ caused by repeated reflections of the GPR signal on shallower reflectors. However, some of them correspond to deep reflection boundaries and indicate that this subsidence has spread a little deeper (Supporting information, Figure S2). In addition, the step-like shape of the strongest reflector on either side of the peat infill suggests that successive small collapses occurred during the peat deposition.

## 5 | DISCUSSION

### 5.1 | Updated hypothesis on the origin of the depressions

Despite three morphogenetic processes that could explain their onset, only two types of depressions have been observed. Those in the SE2 group have irregular shape and proximity to dunes, which are consistent with the aeolian origin by deflation proposed by Sitzia (2014). Blowout hollows indeed develop by wind erosion in dunes (Hesp, 2002) and are typically described as semi-circular or elongated shallow depressions, respectively called “saucer-shaped blowouts” or “trough blowouts” (Cooper, 1958, 1967). In the Landes Triangle, they may result from the erosion of low amplitude dunes (0.3–0.7 m high) of Pleniglacial age (62–42 ka BP) in the northern half, or of more recent parabolic dunes (12.9–11.7 ka BP) in the southern half (Sitzia et al., 2017). Interestingly, SE2 depressions are evenly distributed between these north and south halves (49% vs 51%, respectively).

Yet, surface morphology alone is not sufficient to confirm blowout origin, as blowout hollows can display a high variability of shape and surface, ranging from pits to broad basins (Smith, 1960). The only characteristic criterion is the presence of a deposition lobe and erosional walls, resulting from deflation (Hesp, 2002). González-Villanueva et al. (2011) and Neal & Roberts (2001) have shown that

proper identification of these features can be achieved by GPR. Similarly, sedimentary structures typical of deflation could be investigated in SE2 depressions with GPR at 100 MHz and 250 MHz, which have proven effective for delineating sand beddings near the depression #135.

Depressions in the SE1 group, on the other hand, exhibit characteristics consistent with karst origin, such as circular morphology, distribution in clusters and proximity to karst features, natural streams and faults or folds. Notably, 82% of them lie over shallow potential karst layers buried under permeable sediments- a context favourable for crypto-karst. According to Waltham & Fookes' (2003) classification of closed depressions in carbonate karst, they could correspond to buried dolines or suffosion dolines.

However, their circular morphology, distribution in clusters and location in permeable sediments could also be consistent with a periglacial origin. They specifically show similarities to thermokarst ponds observed in Northern Canada after the thawing of ice-rich permafrost mounds (Bussière et al., 2022; Calmels, Allard, & Delisle, 2008; Fortier & Aubé-Maurice, 2008). The mixture of sand and Brach's clay at the top of the Belin formation of the MPQ complex could act as a frost-susceptible sediment: a necessary component for frost-heave (Konrad, 1999, 2005; Peppin & Style, 2013). Therefore, given the paleoclimatic context of the region, our alternative hypothesis was that the portion of the SE1 group located outside of the karst context may have originated from the thawing of frost mounds during the Last Permafrost Maximum. In such a case, we should be able to detect typical solifluction features and frost-susceptible sediments close to the surface (Bussière et al., 2022; Pissart, 2000).

Near the VL anticline, however, our geophysical survey, aimed at detecting periglacial features, found only evidence of subsidence and no trace of frost-susceptible sediment or solifluction feature typical of permafrost. This new finding contradicts the periglacial hypothesis and therefore, supports the most recent publications in demonstrating the absence of permafrost in southwest France during the Last Glacial Maximum (Saltel et al., 2019; Stadelmaier et al., 2021). The detection of subsidence features definitely eliminates the aeolian hypothesis as well. Alternatively, this close to the anticline axis, the MPQ layers become thinner, the karst layers shallower and the faults more abundant: three elements likely to strengthen water percolation and crypto-karst. The step-like shape of the strongest reflector across depression #135 bears strong similarities to the sedimentary structures obtained under laboratory conditions when simulating ground collapse below a sand deposit (Al-Halbouni et al., 2018; Geyer, Folch, & Martí, 2006). A crypto-karst origin is therefore prone to explain the origin of the depressions near the VL anticline. Deeper geophysical investigation designed to target karst features (like gravimetry) would allow verification.

### 5.2 | Validity of our morphogenetic classification: comparison with previous work

In this section, we verify our rough morphogenetic classification - opposing the SE1 group of karst depressions to the SE2 group of aeolian ones - against local investigations available in the literature. As far as we know, only five academic publications provide observations at the depression scale (Figure 11). Additionally, some dating carried out

<figure><img src="image_11.png" /><figcaption>**FIGURE 10** Zoom on the 250 MHz radargram at the edge of the pond along L1. To the right, the interpretation depicts the termination of weak reflectors (blue arrows) on the strong reflector (dashed red line). As in Figure 9, there is a vertical exaggeration by a factor of two. [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 11:** Zoom of the 250 MHz radargram along L1, showing the termination of minor reflectors (thin blue lines) on the major reflector (thick red line) at the edge of the pond. As in Figure 10, there is a vertical exaggeration by a factor of two.</figcaption></figure>

Yet, surface morphology alone is not sufficient to confirm blowout origin, as blowout hollows can display a high variability of shape and surface, ranging from pits to broad basins (Smith, 1960). The only characteristic criterion is the presence of a deposition lobe and erosional walls, resulting from deflation (Hesp, 2002). González-Villanueva et al. (2011) and Neal and Roberts (2001) have shown that proper identification of these features can be achieved by GPR. Similarly, sedimentary structures typical of deflation could be investigated in SE2 depressions with GPR at 100 MHz and 250 MHz, which have proven effective for delineating sand beddings near the depression #135.

Depressions in the SE1 group, on the other hand, exhibit characteristics consistent with karst origin, such as circular morphology, distribution in cluster, and proximity to karst features, natural streams and faults or folds. Notably, 82% of them lie over shallow potential karst layers buried under permeable sediments- a context favorable for crypto-karst. According to Waltham and Fookes (2003)'s classification of closed depressions in carbonate karst, they could correspond to buried dolines or suffosion dolines.

However, their circular morphology, distribution in cluster, and location in permeable sediments could also be consistent with a periglacial origin. They specifically show similarities to thermokarst ponds observed in Northern Canada after the thawing of ice-rich permafrost mounds (Bussière et al., 2022; Calmels et al., 2008; Fortier and Aubé-Maurice, 2008). The mixture of sand and Brach's clay (at the top of the Belin formation) could act as a frost-susceptible sediment in the MPQ complex: a necessary component for frost-heave (Konrad, 1999, 2005; Peppin and Style, 2013). Therefore,

given the paleoclimatic context of the region, our alternative hypothesis was that the portion of SE1 group located outside of karst context may have originated from the thawing of frost mounds during the Last Permafrost Maximum. In such a case, we should be able to detect typical solifluction features and frost susceptible sediments close to the surface (Bussière et al., 2022; Pissart, 2000).

Near the VL anticline, however, the geophysical survey specifically designed to target periglacial features showed only evidence of subsidence, and no trace of such frost-susceptible sediment (or solifluction feature typical of permafrost). This new finding contradicts the periglacial hypothesis, and therefore supports the most recent publications in demonstrating the absence of permafrost in Southwest France during the Last Glacial Maximum (Saltel et al., 2019; Stadelmaier et al., 2021). The detection of subsidence features definitely eliminates the aeolian hypothesis as well. Alternatively, this close to the anticline axis the MPQ layers become thinner, the karst layers shallower, and the faults more abundant: three elements likely to strengthen water percolation and crypto-karst. The step-like shape of the main reflector across depression #135 bears strong similarities to the sedimentary structures obtained under laboratory conditions when simulating ground collapse below a sand deposit (Geyer et al., 2006). A crypto-karst origin is therefore prone to explain the origin of the depressions near the VL anticline. Deeper geophysical investigation designed to target karst features (like gravimetry) would allow verification.

## 5.2. Validity of our morphogenetic classification: comparison with previous work

In this section we verify our rough morphogenetic classification - opposing the SE1 group of karst depressions to the SE2 group of aeolian ones - against local investigations available in the literature. As far as we know, only five academic publications provide observations at the depression scale (Figure 12). Additionally, some dating carried out by Faure and Galop (2011) help refine the chronology of their formation, if not the process involved.

The depressions in Garein and Brocas municipalities used by Sitzia (2014) as example of deflation hollows are correctly included in the SE2 group.

Likewise, the depressions identified as karstic by Becheler (2015) in Budos and by Saltel (2008) in Saint-Magne all belong to the SE1 group. Specifically, Saltel (2008) mapped the bathymetry of eight large depressions and highlighted their common step-like profile in cross-section, which he interpreted as a succession of karst collapses over time. The infill of one of them (labelled "Les Anguilleyrons") was dated 2,312-2,107 BP by Faure and Galop (2011) at 170 cm depth. Another one ("La Hucau") was exceptionally deep (10 m), reinforcing the interpretation as a doline.

Legigan and Marambat (1993) investigated one depression belonging to the SE2 group ("Bordelounque", in Canenx-et-Réhaut) - a subcircular depression of 150 m in diameter and 1.8 m deep. According to the authors, it displays a slightly asymmetric profile and a sandy surrounding
<figure><img src="image_12.png" /><figcaption>**FIGURE 11** Location of previously studied depressions. Those not included in the inventory or which have disappeared are indicated by arrows. Local names: **1** and **2**: Les Anguilleyrons; **3**: L'Aouarsey; **4**: Eaux belles; **5**: Le Cam; **6**: La Hucau; **7**: Du Bois; **8**: La Honteyre; **11**: Roudrigue; **12**: Bordelounque; **13**: Matalin; **14**: La Hubla. [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

by Faure & Galop (2011) help refine the chronology of their formation, if not the process involved.

The depressions in Garein and Brocas municipalities, identified by Sitzia et al. (2015) as deflation hollows, are included in the SE2 group, supporting the interpretation of an aeolian origin for this group.

Likewise, the depressions identified as karstic by Becheler (2015) in Budos and by Saltel (2008) in Saint-Magne all belong to the SE1 group, reinforcing the interpretation of a karst origin for this subset. Specifically, Saltel (2008) identified step-like profiles in the bathymetry of eight large depressions, which they interpreted as a succession

of karst collapses over time. The infill of one of them (labelled “Les Anguilleyrons”) was dated 2,312–2,107 BP by Faure & Galop (2011) at 170 cm depth. Another one (“La Hucau”) was exceptionally deep (10 m), reinforcing the interpretation as a doline.

Legigan & Marambat (1993) investigated one depression belonging to the SE2 group (“Bordelounque”, in Canenx-et-Réhaut) - a subcircular depression of 150 m in diameter and 1.8 m deep. According to the authors, it displayed a slightly asymmetric profile and a sandy surrounding rim. These elements are again consistent with an aeolian origin, with the surrounding rim possibly corresponding to a
<figure><img src="image_15.png" /></figure>

**Figure 12:** Location of previously studied depressions. Those not included in the inventory or which have disappeared are indicated by arrows. Local names: **1** and **2**: Les Anguilleyrons; **3**: L'Aouarsey; **4**: Eaux belles; **5**: Le Cam; **6**: La Hucau; **7**: Du Bois; **8**: La Honteyre; **11**: Roudrigue; **12**: Bordelounque; **13**: Matalin; **14**: La Hubla.

rim. These elements are again consistent with an aeolian origin, with the surrounding rim possibly corresponding to a deflation lobe. Faure and Galop (2011) later dated its peaty filling, that stretches from 8,391 BP (at 61.5 cm) to 2,470 BP (at 25 cm depth). This implies a much earlier onset compared to the “Anguilleyrons” depression (SE1).

Finally, among the 4 depressions selected for stratigraphic cross-section analysis by Texier (2011), only 2 were included in the present inventory. The first one (“Roudrigue”, in Canenx-et-Réhaut) belongs to SE2. According to the author, it was formed at the expense of an aeolian sand sheet. The cross-section presented no deformation (soft or brittle), but some gley horizons embedded in the sand deposits. Those were interpreted as the result of temporary flooding between low amplitude dunes. Therefore, the aeolian origin of the SE2 group is still consistent. The second one (“du Bois”, in Hostens) has been classified as “unknown”. It is surrounded by another unassigned depression (350 m away), as well as depressions from both the SE1 (700 m away) and SE2 (< 2,500 m away) groups. This reflects a particularly complex context and explains the indecision of the clustering algorithm. However, Texier (2011) clearly identified multiple small collapses and shears along the stratigraphic cross-section, that they interpreted as karst.

The other two depressions studied by Texier (2011), have disappeared due to human activity, but their stratigraphy can provide insights into the processes that occurred in their area, where other depressions still remain. One (“Matalin”, in Canenx-et-Réhaut) presented a surrounding rim, but no stratigraphic deformation. This implies the absence of karst

processes. This depression was located in an area containing dozens of SE2 and “unknown” depressions. It was notably 900 m away from the “Bourdelounque” depression (SE2) - which also presented a sandy rim - and less than 2 km away from the “La Hubla” one, dated by Faure and Galop (2011) around 7,687-7,587 BP at 46 cm depth (before it also disappeared).

The second one (“La Honteyre”, in Hostens) was also located close to depressions currently classified in SE2 group or as “unknown”. However, the stratigraphy displayed numerous deformations and significant change in layer thickness, interpreted as karst features Texier (2011). But the author also noted the presence of a surrounding rim, that we would rather associate with a deflation lobe. The peat filling was dated between 1,100 BP (at 29.5 cm depth) and 4,455 BP (at 98.5 m) (Faure and Galop, 2011). Again, the complexity of the context in consistent with the variability of classes (SE2 and “unknown”) found in the vicinity. It reveals that karst occurred in this area, whereas depressions of SE1 are only found further north, about 3 km away.

Overall, these local-scale investigations demonstrate that the morphogenetic distinction between SE1 and SE2 groups is mainly accurate, even though fine tuning may be necessary in some areas - notably where most of the depressions with “unknown” affiliations are found, or where the geological context is complex (e.g., in Saint-Magne). Additionally, the last example suggests that some old depressions may have been affected by both karstic and aeolian processes, at different periods of time.

## 5.3. Critical thinking on geospatial data limitation and outlooks

The version 2 of BDLISA, which we used here, only provides a qualitative description of the vertical succession of geological layers, but gives no information about the thickness or the depth of each layer. Therefore, the vertical distance from a given depression to the closest potential karst layers beneath could not be evaluated at the regional scale. Instead, we combined the qualitative stratigraphic information with the lateral distance to the closest mapped karst features, as an insight of the karst context. This approximation relies on a local version of the Walter’s law: in the absence of break in the sedimentary sequence, the lateral sequence of facies gives an insight into the vertical one. However, the uncertainty associated to this approximation is very high: it has probably led to an underestimation of the karst potential in some areas, particularly at the edges of the sand cover (which encompass the whole area of the SE1 depressions subset, as well as the border between the SE1 and SE2 subsets). In these areas, the stratigraphy needs to be refined on a case-by-case basis, using borehole stratigraphy. The next generation of BDLISA 3 will soon allow such a refinement, by providing layer thickness estimation.

The precision in spatial extent of the layers is also critical in the karst assessment methodology presented here. The BDLISA has been built at the regional scale on the basis of
deflation lobe. Faure & Galop (2011) later dated its peaty filling, which stretches from 8,391 BP (at 61.5 cm) to 2,470 BP (at 25 cm depth). This implies a much earlier onset compared to the “Anguilleyrons” depression (SE1).

Finally, among the 4 depressions selected for stratigraphic cross-section analysis by Texier (2011), only 2 were included in the present inventory. The first one (“Roudrigue”, in Canenx-et-Réhaut) belongs to SE2. According to the author, it was formed at the expense of an aeolian sand sheet. The cross-section presented no deformation (soft or brittle), but some gley horizons were embedded in the sand deposits. Those were interpreted as the result of temporary flooding between low-amplitude dunes. Therefore, the aeolian origin of the SE2 group is still consistent. The second one (“du Bois”, in Hostens) has been classified as “unknown”. It is surrounded by another unassigned depression (350 m away), as well as depressions from both the SE1 (700 m away) and SE2 (<2,500 m away) groups. This reflects a particularly complex context and explains the indecision of the clustering algorithm. However, Texier (2011) clearly identified multiple small collapses and shears along the stratigraphic cross-section, which they interpreted as karst.

The other two depressions studied by Texier (2011) have disappeared due to human activity, but their stratigraphy can provide insights into the processes that occurred in their area, where other depressions still remain. One (“Matalin”, in Canenx-et-Réhaut) presented a surrounding rim, but no stratigraphic deformation. This implies the absence of karst processes. This depression was located in an area containing dozens of SE2 and “unknown” depressions. It was notably 900 m away from the “Bourdelouique” depression (SE2) - which also presented a sandy rim - and less than 2 km away from the “La Hubla” one, dated by Faure & Galop (2011) around 7,687–7,587 BP at 46 cm depth (before it also disappeared).

The second one (“La Hontyre”, in Hostens) was also located close to depressions currently classified in SE2 group or as “unknown”. However, the stratigraphy displayed numerous deformations and significant changes in layer thickness, interpreted as karst features (Texier, 2011). But the author also noted the presence of a surrounding rim, which we would rather associate with a deflation lobe. The peat filling was dated between 1,100 BP (at 29.5 cm depth) and 4,455 BP (at 98.5 m) (Faure & Galop, 2011). Again, the complexity of the context is consistent with the variability of classes (SE2 and “unknown”) found in the vicinity. It reveals that karst occurred in this area, whereas depressions of SE1 are only found further north, about 3 km away.

Overall, these local-scale investigations demonstrate that the morphogenetic distinction between SE1 and SE2 groups is mainly accurate, even though fine tuning may be necessary in some areas - notably where most of the depressions with “unknown” affiliations are found, or where the geological context is complex (e.g., in Saint-Magne). Additionally, the last example suggests that some old depressions may have been affected by both karstic and aeolian processes at different periods of time.

## 5.3 | Geospatial data limitation and outlooks

Version 2 of BDLISA, which we used here, only provides a qualitative description of the vertical succession of geological layers, but gives no

information about the thickness or the depth of each layer. Therefore, the vertical distance from a given depression to the nearest potential karst layers beneath could not be evaluated at the regional scale. Instead, we combined the qualitative stratigraphic information with the lateral distance to the nearest mapped karst features, as an insight of the karst context. This approximation relies on a local version of Walter’s law: in the absence of a break in the sedimentary sequence, the lateral sequence of facies gives an insight into the vertical one. However, the uncertainty associated to this approximation is very high: it has probably led to an underestimation of the karst potential in some areas, particularly at the edges of the sand cover (which encompass the whole area of the SE1 depressions subset, as well as the border between the SE1 and SE2 subsets). In these areas, the stratigraphy needs to be refined on a case-by-case basis, using borehole stratigraphy. The next generation of BDLISA 3 will soon allow such a refinement by providing layer thickness estimation.

The precision in spatial extent of the layers is also critical in the karst assessment methodology presented here. The BDLISA has been built at the regional scale on the basis of spatial correlation between geological soundings that are unevenly distributed across the Landes Triangle. Despite best efforts, the mapping is inevitably subject to interpretation. As an illustration, if we extend the search for impermeable layers to a radius of 200 m around the depressions when using the decision tree (Figure 3), 6% of the depressions previously classified as in a karst context switch to a non-karst context. Conversely, some depressions located within the extent of an impermeable layer but close to its edges could actually be outside it and subject to crypto-karst, as it seems to be the case for depression #135. This highlights the need to pay special attention to depressions at the boundaries between karst and non-karst contexts.

## 6 | CONCLUSION

At the regional scale, we showed that the Landes Triangle depressions can be classified into two subsets. The first one (SE1, with 1,366 depressions) is consistent with a crypto-karst origin, whereas the second one (SE2, 976 depressions) is likely to be of aeolian origin. Ninety depressions remain unassigned because they had no marked resemblance to either of the subsets. Near the VL anticline, the vast majority belongs to SE1, with 34% are not obviously in a karst context (according to the available data). To explore an alternative periglacial origin, we conducted a geophysical survey specifically designed to target periglacial features beneath the sand layers, such as solifluction lobes in frost-susceptible silt. Nevertheless, the geophysical images revealed only evidence of subsidence, with no indication of frost-susceptible sediment or thermokarst features. This aligns with recent paleoclimatic results, which show that climatic conditions in the Landes Triangle during the LGM were too mild for permafrost to develop. We conclude that a karst origin is the most plausible explanation for the formation of most depressions near the VL anticline, and that further geological investigation is needed to capture the extent of karst and map the related water pathways in that area.

More broadly, our methodology offers a robust framework for distinguishing between geomorphological processes. By integrating statistical analyses with geophysical imaging, it revealed the presence
spatial correlation between geological sounding that are unevenly distributed across the Landes Triangle. Despite best efforts, the mapping is inevitably subject to interpretation. As an illustration, if we extend the search for impermeable layers to a radius of 200 m around the depressions when using the decision tree (Figure 3), 6% of the depressions previously classified as in a karst context switch to a non-karst context. Conversely, some depressions located within the extent of an impermeable layer but close to its edges could actually be outside it and subject to crypto-karst, as it seems to be the case for depression #135. This highlights the need to pay special attention to depressions at the boundaries between karst and non-karst contexts.

## 6. Conclusion

At the regional scale, we showed that the Land Triangle depressions can be classified in two subsets. The first one (SE1, with 1,366 depressions) is consistent with a crypto-karst origin, whereas the second one (SE2, 976 depressions) is likely to be of aeolian origin. Ninety depressions have been left unassigned, because they bore no marked resemblance to either of the subsets. Near the VL anticline, the vast majority belongs to SE1, and 34% are not obviously in a karst context (according to the available data). We thus tested the alternative hypothesis of a periglacial origin. Nevertheless, the geophysical survey specifically designed to target periglacial features in this context showed only evidence of subsidence, and none of frost-susceptible sediment or thermokarst feature. This is consistent with the most recent paleoclimatic results, which show that climatic conditions in the Landes Triangle during the LGM were too mild for permafrost to develop. We conclude that a karst origin is prone to explain the formation of most of the depressions near the VL anticline, and that further geological investigation is needed to capture the extent of karst and map the related water pathways in that area. The existence of such pathways is indicative of a rapid connection between surface water and groundwater. As a consequence, the groundwater recharge is likely to be vulnerable to changes in surface hydroclimatic context. Lastly, our overall methodology - which combines statistical analysis of conscientiously chosen variables with geophysical measurements - can be extended to the study of any pond-like ecosystem in order to assess its involvement in groundwater recharge.

## Acknowledgments

The research presented herein was supported by the Initiative d'Excellence program of Université de Bordeaux and by the Sentinel North program of Université Laval. The authors are also grateful to UMR 7619 METIS for the loan of their GPR and for their pertinent advice; to the CBNSA, CEN, CD40 and PNRLG for sharing the inventories of depressions; to J. Paris (UMR 5805 EPOC) for taking part in GPR acquisition; and to P.Bertran (UMR 5805 EPOC) for his careful reading and pertinent comments on the manuscript.

## Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## Conflict of interest statement

The authors declare that they have no conflict of interest regarding the work presented here.

## CRediT authorship contribution statement

**Léa Bussière:** Conceptualization, Methodology, Investigation, Formal analysis, Visualization, Writing - original draft. **Myriam Schmutz:** Conceptualization, Methodology, Funding acquisition, Resources, Supervision, Writing - review & editing. **Alain Dupuy:** Funding acquisition, Project administration, Supervision, Writing - review & editing.

## References

Andrieux, E., Bertran, P., Antoine, P., Deschodt, L., Lenoble, A., Coutard, S., Ajas, A., Borderie, Q., Coutard, J.P., Didierjean, F., et al., 2016a. Database of pleistocene periglacial features in France: description of the online version. Quaternaire. Revue de l'Association française pour l'étude du Quaternaire 27, 329-339. doi:10.4000/quaternaire.7717.

Andrieux, E., Bertran, P., Saito, K., 2016b. Spatial analysis of the French Pleistocene permafrost by a GIS database. Permafrost and periglacial processes 27, 17-30. doi:10.1002/ppp.1856.

Annan, A.P., 2005. GPR methods for hydrogeological studies, in: Rubin, Y., Hubbard, S. (Eds.), Hydrogeophysics. Springer, Dordrecht. chapter 7, pp. 185-213. doi:10.1007/1-4020-3102-5_7.

Becheler, P., 2014. L'origine tectono-karstique des lagunes dans la région de Villagrains-Landiras. (In French). L'échos des faluns 35-36, 11-12.

Becheler, P., 2015. L'origine tectono-karstique des lagunes dans la région de Villagrains-Landiras (suite). (In French). L'échos des faluns 37-38, 17.

Bertran, P., Andrieux, E., Antoine, P., Coutard, S., Deschodt, L., Gardère, P., Hernandez, M., Legentil, C., Lenoble, A., Liard, M., et al., 2014. Distribution and chronology of Pleistocene permafrost features in France: Database and first results. Boreas 43, 699-711. doi:10.1111/bor.12025.

Bertran, P., Andrieux, E., Antoine, P., Deschodt, L., Font, M., Sicilia, D., 2017. Pleistocene involutions and patterned ground in France: examples and analysis using a GIS database. Permafrost and periglacial processes 28, 710-725. doi:10.1002/ppp.1957.

Bertran, P., Andrieux, E., Bateman, M.D., Fuchs, M., Klinge, M., Marembert, F., 2020. Mapping and chronology of coversands and dunes from the Aquitaine basin, southwest France. Aeolian research 47, 100628. doi:10.1016/j.aeolia.2020.100628.

Binley, A., 2015. Tools and techniques: Electrical methods, in: Schubert, G. (Ed.), Treatise on Geophysics. 2nd ed.. Elsevier, Oxford. volume 11. chapter 8, pp. 233-259. doi:10.1016/B978-0-444-53802-4.00192-5.

Binley, A., Kemna, A., 2005. DC resistivity and induced polarization methods, in: Rubin, Y., Hubbard, S. (Eds.), Hydrogeophysics. Springer, Dordrecht. chapter 5, pp. 129-156. doi:10.1007/1-4020-3102-5_5.

Blott, S.J., Pye, K., 2001. GRADISTAT: a grain size distribution and statistics package for the analysis of unconsolidated sediments. Earth surface processes and Landforms 26, 1237-1248. doi:10.1002/esp.261.

Bourbon, P., 2019. Gestion des eaux souterraines en Région Aquitaine - Carte géologique numérique à 1/250 000 de la région Aquitaine. Notice technique - Module 1.5 - Année 2. (In French). Final report. BRGM/RP-68881-FR.
of karst depressions in a key recharge area of a regional aquifer. This implies the presence of rapid surface-to-groundwater connections and emphasizes the vulnerability of the groundwater resource to climate change. As this approach relies on commonly available GIS data and rapid, non-invasive geophysical methods, it can be extended to assess the role of any pond-like ecosystems in groundwater recharge.

## ACKNOWLEDGEMENTS AND FUNDING

The research presented herein was supported by the Initiative d'Excellence program of Université de Bordeaux and by the Sentinel North program of Université Laval. The authors are also grateful to UMR 7619 METIS for the loan of their GPR and for their pertinent advice; to the CBNSA, CEN, CD40 and PNRLG for sharing the inventories of depressions; to J. Paris (UMR 5805 EPOC) for taking part in GPR acquisition; and to P. Bertran (UMR 5199 PACEA and INRAP) for his careful reading and pertinent comments on the early version of the manuscript. We thank two anonymous reviewers for their detailed comments and suggestions for improving this study.

## CONFLICT OF INTEREST STATEMENT

The authors declare that they have no conflict of interest regarding the work presented here.

## ACKNOWLEDGEMENTS

Open Access funding enabled and organized by Projekt DEAL.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ORCID

Léa Bussière | https://orcid.org/0000-0002-2321-9449

## REFERENCES

Al-Halbouni, D., Holohan, E.P., Taheri, A., Schöpfer, M.P., Emam, S. & Dahm, T. (2018) Geomechanical modelling of sinkhole development using distinct elements: model verification for a single void space and application to the Dead Sea area. *Solid Earth*, 9(6), 1341–1373. Available from: https://doi.org/10.5194/se-9-1341-2018.

Andrieux, E., Bertran, P., Antoine, P., Deschodt, L., Lenoble, A., Coutard, S., et al. (2016) Database of pleistocene periglacial features in France: description of the online version. *Quaternaire. Revue de L'association française Pour L'étude du Quaternaire*, 27(4), 329–339. Available from: https://doi.org/10.4000/quaternaire.7717.

Andrieux, E., Bertran, P. & Saito, K. (2016) Spatial analysis of the French Pleistocene permafrost by a GIS database. *Permafrost and Periglacial Processes*, 27(1), 17–30. Available from: https://doi.org/10.1002/ppp.1856.

Arlen-Pouliot, Y. & Bhiry, N. (2005) Palaeoecology of a palsa and a filled thermokarst pond in a permafrost peatland, subarctic Québec, Canada. *The Holocene*, 15(3), 408–419. Available from: https://doi.org/10.1191/0959683605hl818rp.

Ballantyne, C.K. (2018) *Periglacial geomorphology*. Chichester, UK: Wiley Blackwell, p. 454 https://doi.org/10.1002/ppp.2108.

Becheler, P. (2014) L'origine tectono-karstique des lagunes dans la région de Villagrains-Landiras. (In French). *L'échos Des Faluns*, 35, 11–12.

Becheler, P. (2015) L'origine tectono-karstique des lagunes dans la région de Villagrains-Landiras (suite). (In French). *L'échos Des Faluns*, 37, 17.

Bertran, P., Andrieux, E., Antoine, P., Coutard, S., Deschodt, L., Gardère, P., et al. (2014) Distribution and chronology of pleistocene permafrost features in France: database and first results. *Boreas*, 43(3), 699–711. Available from: https://doi.org/10.1111/bor.12025.

Bertran, P., Andrieux, E., Antoine, P., Deschodt, L., Font, M. & Sicilia, D. (2017) Pleistocene involutions and patterned ground in France: examples and analysis using a GIS database. *Permafrost and Periglacial Processes*, 28(4), 710–725. Available from: https://doi.org/10.1002/ppp.1957.

Bertran, P., Andrieux, E., Bateman, M.D., Fuchs, M., Klinge, M. & Marembert, F. (2020) Mapping and chronology of coversands and dunes from the Aquitaine basin, Southwest France. *Aeolian Research*, 47, 100628. Available from: https://doi.org/10.1016/j.aeolia.2020.100628.

Binley, A. (2015) Tools and Techniques: Electrical Methods. In: Schubert, G. (Ed.) *Treatise on geophysics*, 2nd edition, Vol. 11. Oxford: Elsevier, pp. 233–259 https://doi.org/10.1016/B978-0-444-53802-4.00192-5.

Binley, A. & Kemna, A. (2005) DC resistivity and induced polarization methods. In: Rubin, Y. & Hubbard, S.S. (Eds.) *Hydrogeophysics. Dordecht: Springer*, pp. 129–156 https://doi.org/10.1007/1-4020-3102-5-5.

Blott, S.J. & Pye, K. (2001) GRADISTAT: a grain size distribution and statistics package for the analysis of unconsolidated sediments. *Earth Surface Processes and Landforms*, 26(11), 1237–1248. Available from: https://doi.org/10.1002/esp.261.

Bourbon P. (2019) Gestion des eaux souterraines en Région Aquitaine - Carte géologique numérique à 1/250 000 de la région Aquitaine. Notice technique - Module 1.5 - Année 2. (In French). Final report. BRGM/RP-68881- FR

Brosius, L.S., Anthony, K.W., Treat, C.C., Lenz, J., Jones, M.C., Bret-Harte, M.S., et al. (2021) Spatiotemporal patterns of northern lake formation since the last glacial maximum. *Quaternary Science Reviews*, 253, 106773. Available from: https://doi.org/10.1016/j.quascirev.2020.106773.

Bussière, L., Schmutz, M., Fortier, R., Lemieux, J.M. & Dupuy, A. (2022) Near-surface geophysical imaging of a thermokarst pond in the discontinuous permafrost zone in Nunavik (Québec), Canada. *Permafrost and Periglacial Processes*, 33(4), 353–369. Available from: https://doi.org/10.1002/ppp.2166.

Cabaret O., Husson E., Ladouche B., Dewandel B., Baudement C. (2020) Synthèse hydrogéologique de la structure anticlinale d'Audignon. (In French). Final report. BRGM/RG-69926-FR

Calmels, F., Allard, M. & Delisle, G. (2008) Development and decay of a lithalsa in northern Quebec: a geomorphological history. *Geomorphology*, 97(3-4), 287–299. Available from: https://doi.org/10.1016/j.geomorph.2007.08.013.

Cooper, W.S. (1958) Coastal sand dunes of Oregon and Washington. *Geological Society of America Memoirs*, 72, 169. Available from: https://doi.org/10.1130/MEM72-p1.

Cooper, W.S. (1967) Coastal sand dunes of California. *Geological Society of America Memoirs*, 101, 131. Available from: https://doi.org/10.1130/MEM104-p1.

Corbier, P., & Karnay, G. with the collaboration of Bourbine, B. & Slatel, M. (2010) Gestion des eaux souterraines en région Aquitaine - Reconnaissance des potentialités aquifères du Mio-Plio-Quaternaire des Landes de Gascogne et du Médoc en relation avec le SAGE - Module 7 - Année 1. (In French). Final report. BRGM/RP-57813-FR.

Corbier, P. with the collaboration of Bonneau, C. A. (2015) Caractéristiques géologiques et hydrogéologiques des lagunes landaises. (In French). Final report. BRGM/RP-64633-FR.

Doctor, D.H. & Doctor, K.Z. (2012) Spatial analysis of geologic and hydrologic features relating to sinkhole occurrence in Jefferson County, West Virginia. *Carbonates and Evaporites*, 27(2), 143–152. Available from: https://doi.org/10.1007/s13146-012-0098-1.

Enjalbert, H. (1960) Les pays aquitains. Le modelé et les sols. (In French). Bordeaux.

Etienne, D., Ruffaldi, P., Goepp, S., Ritz, F., Georges-Leroy, M., Pollier, B., et al. (2011) The origin of closed depressions in northeastern France: a new assessment. *Geomorphology*, 126(1-2), 121–131. Available from: https://doi.org/10.1016/j geomorph.2010.10.036.

Farrant, A.R. & Cooper, A.H. (2008) Karst geohazards in the UK: the use of digital data for hazard management. *Quarterly Journal of Engineering*
Bussière, L., Schmutz, M., Fortier, R., Lemieux, J.M., Dupuy, A., 2022. Near-surface geophysical imaging of a thermokarst pond in the discontinuous permafrost zone in Nunavik (Québec), Canada. Permafrost and Periglacial Processes 33, 353–369. doi:10.1002/ppp.2166.

Cabaret, O., Husson, E., Ladouche, B., Dewandel, B., Baudement, C., 2020. Synthèse hydrogéologique de la structure anticlinale d'Audignon. (In French). Final report. BRGM/RG-69926-FR.

Calmels, F., Allard, M., Delisle, G., 2008. Development and decay of a lithalsa in Northern Quebec: a geomorphological history. Geomorphology 97, 287–299. doi:10.1016/j geomorph.2007.08.013.

Cooper, W.S., 1958. Coastal sand dunes of Oregon and Washington. Geological Society of America - Memoir 72. doi:10.1130/MEM72-p1.

Cooper, W.S., 1967. Coastal sand dunes of California. Geological Society of America - Memoir 101. doi:10.1130/MEM104-p1.

Corbier, P., Karnay, G., with the collaboration of Bourgine, B., & Slatel, M., 2010. Gestion des eaux souterraines en région Aquitaine - Reconnaisance des potentialités aquifères du Mio-Plio-Quaternaire des Landes de Gascogne et du Médoc en relation avec le SAGE - Module 7 - Année 1. (In French). Final report. BRGM/RP-57813-FR.

Corbier, P., with the collaboration of Bonneau, C. A., 2015. Caractéristiques géologiques et hydrogéologiques des lagunes landaises. (In French). Final report. BRGM/RP-64633-FR.

Doctor, D.H., Doctor, K.Z., 2012. Spatial analysis of geologic and hydrologic features relating to sinkhole occurrence in Jefferson County, West Virginia. Carbonates and Evaporites 27, 143–152. doi:10.1007/s13146-012-0098-1.

Enjalbert, H., 1960. Les pays aquitains. Le modelé et les sols. (In French).

Everett, M.E., 2013. Ground-penetrating radar, in: Near-surface applied geophysics. Cambridge University Press, chapter 9, pp. 239–278. doi:10.1017/CBO9781139088435.

Farrant, A., Cooper, A., 2008. Karst geohazards in the UK: the use of digital data for hazard management. Quarterly Journal of Engineering Geology and Hydrogeology 41, 339–356. doi:10.1144/1470-9236/07-201.

Faure, E., Galop, D., 2011. La fin du paradigme du désert landais : histoire de la végétation et de l'anthropisation à partir de l'étude palynologique de quelques lagunes de la Grande-Lande. (In French). Aquitania 24, 43–59.

Fortier, R., Aubé-Maurice, B., 2008. Fast permafrost degradation near Umiujaq in Nunavik (Canada) since 1957 assessed from time-lapse aerial and satellite photographs, in: Proceedings of Ninth International Conference on Permafrost, pp. 457–462.

Galve, J.P., Gutiérrez, F., Remondo, J., Bonachea, J., Lucha, P., Cendrero, A., 2009. Evaluating and comparing methods of sinkhole susceptibility mapping in the Ebro Valley evaporite karst (NE Spain). Geomorphology 111, 160–172. doi:10.1016/j geomorph.2009.04.017.

Geyer, A., Folch, A., Martí, J., 2006. Relationship between caldera collapse and magma chamber withdrawal: An experimental approach. Journal of Volcanology and Geothermal Research 157, 375–386. doi:10.1016/j.jvolgeores.2006.05.001.

González-Villanueva, R., Costas, S., Duarte, H., Pérez-Arlucea, M., Alejo, I., 2011. Blowout evolution in a coastal dune: using GPR, aerial imagery and core records, in: Journal of Coastal Research, Poland. pp. 278–282. ISSN 0749-0208.

Hesp, P., 2002. Foredunes and blowouts: initiation, geomorphology and dynamics. Geomorphology 48, 245–268. doi:10.1016/S0169-555X(02)00184-8.

James, G., Witten, D., Hastie, T., Tibshirani, R., 2021. Unsupervised learning, in: An Introduction to Statistical Learning with Applications in R. 2 ed.. Springer Texts in Statistics, New-York. chapter 12, pp. 497–552. doi:10.1007/978-1-0716-1418-1_12.

Konrad, J.M., 1999. Frost susceptibility related to soil index properties. Canadian Geotechnical Journal 36, 403–417. doi:10.1139/t99-008.

Konrad, J.M., 2005. Estimation of the segregation potential of fine-grained soils using the frost heave response of two reference soils. Canadian Geotechnical Journal 42, 38–50. doi:https://doi.org/10.1139/t04-080.

Labat, C., 2021. Fonctionnement hydrogéologique d'un système aquifère multicouche aux abords d'une struture anticlinale : le cas de l'anticlinal de Villagrains-Landiras en Gironde : méthodologie et apport d'une

approche pluridisciplinaire. (In French). Ph.D. thesis. Université Michel de Montaigne - Bordeaux III. URL: https://theses.hal.science/tel-03549437. French repository of doctoral dissertations.

Labat, C., Larroque, F., de Grissac, B., Dupuy, A., Saltel, M., Bourbon, P., 2021. Influence of an anticline structure on hydrogeological functioning and aquifer interactions in a multilayered aquifer system: the case of Villagrains-Landiras anticline (Gironde, France). Hydrogeology Journal 29, 1711–1732. doi:10.1007/s10040-021-02333-z.

Lamelas, M., Marinoni, O., Hoppe, A., De La Riva, J., 2008. Doline probability map using logistic regression and GIS technology in the central Ebro Basin (Spain). Environmental Geology 54, 963–977. doi:10.1007/s00254-007-0895-3.

Legigan, P., 1979. Elaboration de la formation du Sable des Landes : dépôt résiduel de l'environnement sédimentaire pliocène-pléistocène centre aquitaine. (In French). Ph.D. thesis. Université Bordeaux 1.

Legigan, P., Marambat, L., 1993. Age de la formation d'une lagune landaise : premières données palynologiques et radiométriques. (In French). Bulletin de la Société de Borda 118, 433–443.

Lindgren, A., Hugelius, G., Kuhry, P., Christensen, T.R., Vandenberghe, J., 2016. GIS-based maps and area estimates of Northern Hemisphere permafrost extent during the last glacial maximum. Permafrost and Periglacial Processes 27, 6–16. doi:10.1002/ppp.1851.

Loke, M., Barker, R., 1995. Least-squares deconvolution of apparent resistivity pseudosections. Geophysics 60, 1682–1690. doi:10.1190/1.1443900.

Loke, M.H., Barker, R.D., 1996. Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-Newton method. Geophysical prospecting 44, 131–152. doi:10.1111/j.1365-2478.1996.tb00142.x.

Mariani, M.C., Tweneboah, O.K., Beccar-Varela, M.P., 2021a. Discriminant and Cluster analysis, in: Data science in theory and practice: techniques for Big Data analytics and complex data sets. John Wiley & Sons. chapter 12, p. 165–178.

Mariani, M.C., Tweneboah, O.K., Beccar-Varela, M.P., 2021b. Principal Component Analysis, in: Data science in theory and practice: techniques for Big Data analytics and complex data sets. John Wiley & Sons. chapter 11, p. 151–163.

Neal, A., Roberts, C.L., 2001. Internal structure of a trough blowout, determined from migrated ground-penetrating radar profiles. Sedimentology 48, 791–810. doi:10.1046/j.1365-3091.2001.00382.x.

Oldenburg, D.W., Li, Y., 1999. Estimating depth of investigation in DC resistivity and IP surveys. Geophysics 64, 403–416. doi:10.1190/1.1444545.

Orndorff, R.C., Weary, D.J., Lagueux, K.M., 2000. Geographic information systems analysis of geologic controls on the distribution on dolines in the Ozarks of south-central Missouri, USA. Acta Carsologica 29. doi:10.3986/ac.v29i2.456.

Peppin, S.S., Style, R.W., 2013. The physics of frost heave and ice-lens growth. Vadose Zone Journal 12, 1–12. doi:10.2136/vzj2012.0049.

Perrin, J., Noury, G., Cartannaz, C., ALEA-karst : vers une approche multicritères de l'aléa lié aux cavités karstiques en contexte calcaire. (In French). Volume 1: Trame méthodologique. Report. BRGM/RP-63771-FR, 77 pages.

Pissart, A., 2000. Remnants of lithalsas of the Hautes Fagnes, Belgium: a summary of present-day knowledge. Permafrost and periglacial Processes 11, 327–355. doi:10.1002/1099-1530(200012)11:4%3C327::AID-PPP370%3E3.0.CO;2-Q.

Platel, J.P., Dubreuilh, J., 1986. Le bassin d'Aquitaine, in: Lajoinie, J.P., Platel, r. (Eds.), Grandes étapes de l'histoire géologique de quatre régions françaises. (In French). BRGM Report 86-SGN-424-GEO, 18 p., 12 fig.

Saltel, M., 2008. Impact de structures géologiques sur l'alimentation de systèmes aquifères profonds - Fonctionnement hydrogéologique des antiformes du sud de Bordeaux. (In French). Ph.D. thesis. Université Michel de Montaigne - Bordeaux III. URL: https://theses.fr/2008BOR30062. French repository of doctoral dissertations.

Saltel, M., Rebeix, R., Thomas, B., Franceschi, M., Lavielle, B., Bertran, P., 2019. Paleoclimate variations and impact on groundwater recharge in multi-layer aquifer systems using a multi-tracer approach (northern
*Geology and Hydrogeology*, 41(3), 339–356. Available from: https://doi.org/10.1144/1470-9236/07-201.

Faure, E. & Galop, D. (2011) La fin du paradigme du désert landais: histoire de la végétation et de l'anthropisation à partir de l'étude palynologique de quelques lagunes de la Grande-Lande (In French). *Aquitania*, 24, 43–59.

Ford, D. & Williams, P.D. (2007) Karst Hydrogeology. In: Ford & Williams. (Ed.) *Karst hydrogeology and geomorphology*. Chichester, UK: John Wiley & Sons, pp. 103–144 https://doi.org/10.1002/9781118684986.ch5.

Fortier, R. & Aubé-Maurice, B. (2008) Fast permafrost degradation near Umiujaq in Nunavik (Canada) since 1957 assessed from time-lapse aerial and satellite photographs. *Proceedings of Ninth International Conference on Permafrost*, 9, 457–462.

French, H.M. & Demitroff, M. (2001) Cold-climate origin of the enclosed depressions and wetlands ('spungs') of the pine barrens, southern New Jersey, USA. *Permafrost and Periglacial Processes*, 12(4), 337–350. Available from: https://doi.org/10.1002/ppp.401.

Galve, J.P., Gutiérrez, F., Remondo, J., Bonachea, J., Lucha, P. & Cendrero, A. (2009) Evaluating and comparing methods of sinkhole susceptibility mapping in the Ebro Valley evaporite karst (NE Spain). *Geomorphology*, 111(3-4), 160–172. Available from: https://doi.org/10.1016/j geomorph.2009.04.017.

Geyer, A., Folch, A. & Martí, J. (2006) Relationship between caldera collapse and magma chamber withdrawal: an experimental approach. *Journal of Volcanology and Geothermal Research*, 157(4), 375–386. Available from: https://doi.org/10.1016/j.jvolgeores.2006.05.001.

Gökkaya, E., Gutiérrez, F., Ferk, M. & Görüm, T. (2021) Sinkhole development in the Sivas gypsum karst, Turkey. *Geomorphology*, 386, 107746. Available from: https://doi.org/10.1016/j geomorph.2021.107746.

González-Villanueva, R., Costas, S., Duarte, H., Pérez-Arlucea, M. & Alejo, I. (2011) Blowout evolution in a coastal dune: using GPR, aerial imagery and core records. *Journal of Coastal Research*, Special Issue 64: *Proceedings of the 11th International Coastal Symposium ICS2011*, 278–282. ISSN 0749-0208.

Goudie, A.S. & Wells, G.L. (1995) The nature, distribution and formation of pans in arid zones. *Earth-Science Reviews*, 38(1), 1–69. Available from: https://doi.org/10.1016/0012-8252(94)00066-6.

Grosse, G., Jones, B., Arp, C. 2013. Thermokarst lakes, drainage, and drained basins. In: Schroder, J. (Editor in Chief), Giardino, R., Harbor, J. (Eds.), *Treatise on geomorphology*. Academic Press, San Diego, CA, vol. 8, Glacial and Periglacial Geomorphology, pp. 325–353. https://doi.org/10.1016/B978-0-12-374739-6.00216-5

Harris, S.A., Brouchkov, A. & Guodong, C. (2017) *Geomorphology: characteristics and use of frozen ground and permafrost landforms*. London, UK: CRC press https://doi.org/10.4324/9781315166988.

Hesp, P. (2002) Foredunes and blowouts: initiation, geomorphology and dynamics. *Geomorphology*, 48(1–3), 245–268. Available from: https://doi.org/10.1016/S0169-555X(02)00184-8.

Hugenholtz, C.H. & Wolfe, S.A. (2006) Morphodynamics and climate controls of two aeolian blowouts on the northern great plains, Canada. *Earth Surface Processes and Landforms*, 31(12), 1540–1557. Available from: https://doi.org/10.1002/esp.1367.

James, G., Witten, D., Hastie, T. & Tibshirani, R. (2021) Unsupervised Learning. In: James, G. et al. (Eds.) *An introduction to statistical learning with applications in R*, Second edition. New-York: Springer Texts in Statistics, pp. 497–552 https://doi.org/10.1007/978-1-0716-1418-1_12.

Jasechko, S. (2016) Partitioning young and old groundwater with geochemical tracers. *Chemical Geology*, 427, 35–42. Available from: https://doi.org/10.1016/j.chemgeo.2016.02.012.

Jiráková, H., Huneau, F., Celle-Jeanton, H., Hrkal, Z. & Le Coustumer, P. (2011) Insights into palaeorecharge conditions for European deep aquifers. *Hydrogeology Journal*, 19(8), 1545–1562. Available from: https://doi.org/10.1007/s10040-011-0765-7.

Keskin, I. & Yilmaz, I. (2016) Morphometric and geological features of karstic depressions in gypsum (Sivas, Turkey). *Environmental Earth Sciences*, 75(12), 1040. Available from: https://doi.org/10.1007/s12665-016-5845-5.

Kolodyńska-Gawrysiak, R. & Poesen, J. (2017) Closed depressions in the European loess belt-natural or anthropogenic origin? *Geomorphology*, 288, 111–128. Available from: https://doi.org/10.1016/j.geomorph.2017.02.004.

Konrad, J.M. (1999) Frost susceptibility related to soil index properties. *Canadian Geotechnical Journal*, 36(3), 403–417. Available from: https://doi.org/10.1139/t99-008.

Konrad, J.M. (2005) Estimation of the segregation potential of fine-grained soils using the frost heave response of two reference soils. *Canadian Geotechnical Journal*, 42(1), 38–50. Available from: https://doi.org/10.1139/t04-080.

Kranjc, A. (2013). Classification of closed depressions in carbonate karst. *Shroder, J. (Editor in Chief), Frumkin, A. (Eds.), Treatise on geomorphology*. Academic Press, San Diego, CA, 6, Karst Geomorphology, 104–111. https://doi.org/10.1016/B978-0-12-374739-6.00125-1

Labat C. (2021) Fonctionnement hydrogéologique d'un système aquifère multicouche aux abords d'une structure anticlinale: le cas de l'anticlinal de Villagrains-Landiras en Gironde: méthodologie et apport d'une approche pluridisciplinaire. (In French). [Doctoral dissertation, Université Michel de Montaigne - Bordeaux III]. French repository of doctoral dissertations. https://theses.hal.science/tel-03549437

Labat, C., Larroque, F., de Grissac, B., Dupuy, A., Saltel, M. & Bourbon, P. (2021) Influence of an anticline structure on hydrogeological functioning and aquifer interactions in a multilayered aquifer system: the case of Villagrains-Landiras anticline (Gironde, France). *Hydrogeology Journal*, 29(5), 1–22. Available from: https://doi.org/10.1007/s10040-021-02333-z.

Lamelas, M.T., Marinoni, O., Hoppe, A. & De La Riva, J. (2008) Doline probability map using logistic regression and GIS technology in the Central Ebro Basin (Spain). *Environmental Geology*, 54(5), 963–977. Available from: https://doi.org/10.1007/s00254-007-0895-3.

Legigan P. (1979) Elaboration de la formation du Sable des Landes: dépôt résiduel de l'environnement sédimentaire, pliocène-pléistocène centre aquitaine. (In French). [Doctoral dissertation, Université Bordeaux 1]

Legigan, P. & Marambat, L. (1993) Age de la formation d'une lagune landaise: premières données palynologiques et radiométriques (In French). *Bulletin de la Société de Borda*, 118(432), 433–443.

Letz, O., Siebner, H., Avrahamov, N., Egozi, R., Eshel, G. & Dahan, O. (2021) The impact of geomorphology on groundwater recharge in a semi-arid mountainous area. *Journal of Hydrology*, 603, 127029. Available from: https://doi.org/10.1016/j.jhydrol.2021.127029.

Lindgren, A., Hugelius, G., Kuhry, P., Christensen, T.R. & Vandenberghe, J. (2016) GIS-based maps and area estimates of northern hemisphere permafrost extent during the last glacial maximum. *Permafrost and Periglacial Processes*, 27(1), 6–16. Available from: https://doi.org/10.1002/ppp.1851.

Loke, M.H. & Barker, R.D. (1995) Least-squares deconvolution of apparent resistivity pseudosections. *Geophysics*, 60(6), 1682–1690. Available from: https://doi.org/10.1190/1.1443900.

Loke, M.H. & Barker, R.D. (1996) Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-Newton method. *Geophysical Prospecting*, 44(1), 131–152. Available from: https://doi.org/10.1111/j.1365-2478.1996.tb00142.x.

Majecka, A., Forysiak, J., Marks, L. & Tołoczko-Pasek, A. (2019) Lithological diversity of the deposits of closed depressions in Central Poland as a result of their origin conditions. *Quaternary International*, 501, 208–218. Available from: https://doi.org/10.1016/j.quaint.2018.02.018.

Mariani, M.C., Tweneboah, O.K. & Beccar-Varela, M.P. (2021) Discriminant and Cluster Analysis. In: Mariani, M.C. et al. (Eds.) *Data science in theory and practice: techniques for big data analytics and complex data sets*. Hoboken, USA: John Wiley & Sons, pp. 165–178.

Neal, A. & Roberts, C.L. (2001) Internal structure of a trough blowout, determined from migrated ground-penetrating radar profiles. *Sedimentology*, 48(4), 791–810. Available from: https://doi.org/10.1046/j.1365-3091.2001.00382.x.

Oldenburg, D.W. & Li, Y. (1999) Estimating depth of investigation in dc resistivity and IP surveys. *Geophysics*, 64(2), 403–416. Available from: https://doi.org/10.1190/1.1444545.
Aquitaine basin, France). Hydrogeology journal 27, 1439-1457. doi:10.1007/s10040-019-01944-x.

Sinclair, W.C., Stewart, J.W., 1985. Sinkhole type, development, and distribution in Florida. Published online. URL: http://ufdc.ufl.edu/UF90000362/00001/1x.

Sitzia, L., 2014. Chronostratigraphie et distribution spatiale des dépôts éoliens quaternaires du Bassin Aquitain. (In French). Ph.D. thesis. Université de Bordeaux. URL: https://theses.hal.science/tel-01197192. French repository of doctoral dissertations.

Sitzia, L., Bertran, P., Sima, A., Chery, P., Queffelec, A., Rousseau, D.D., 2017. Dynamics and sources of last glacial aeolian deposition in southwest France derived from dune patterns, grain-size gradients and geochemistry, and reconstruction of efficient wind directions. Quaternary Science Reviews 170, 250-268. doi:10.1016/j.quascirev.2017.06.029.

Smith, H., 1960. Physiography and photo interpretation of coastal sand dunes. Final Report Contract NONR-2242 (00), Office of Naval Research, Geographical Branch.

Stadelmaier, K.H., Ludwig, P., Bertran, P., Antoine, P., Shi, X., Lohmann, G., Pinto, J.G., 2021. A new perspective on permafrost boundaries in france during the Last Glacial Maximum. Climate of the Past 17, 2559-2576. doi:10.5194/cp-17-2559-2021.

Texier, J., 2011. Genèse des lagunes landaises: un point sur la question, in: Merlet, J., Bost, J. (Eds.), De la lagune à l'airial. Le peuplement de la Grande-Lande. (In French). Pessac, p. 23-42. Aquitania suppl. 24.

Trask, P.D., 1930. Mechanical analyses of sediments by centrifuge. Economic Geology 25, 581-599. doi:10.2113/gsecongeo.25.6.581.

Vandenberghe, J., French, H.M., Gorbunov, A., Marchenko, S., Velichko, A.A., Jin, H., Cui, Z., Zhang, T., Wan, X., 2014. The Last Permafrost Maximum (LPM) map of the Northern Hemisphere: permafrost extent and mean annual air temperatures, 25-17 ka BP. Boreas 43, 652-666. doi:10.1111/bor.12070.

Waltham, A., Fookes, P., 2003. Engineering classification of karst ground conditions. Quarterly Journal of Engineering Geology and Hydrogeology 36, 101-118. doi:10.1144/1470-9236/2002-33.
Orndorff, R.C., Weary, D.J. & Lagueux, K.M. (2000) Geographic information systems analysis of geologic controls on the distribution on dolines in the Ozarks of south-Central Missouri, USA. *Acta Carsologica*, 29(2), 161-175. Available from: https://doi.org/10.3986/ac.v29i2.456.

Pardo-Igúzquiza, E. & Dowd, P.A. (2021) The mapping of closed depressions and its contribution to the geodiversity inventory. *International Journal of Geoheritage and Parks*, 9(4), 480-495. Available from: https://doi.org/10.1016/j.ijgeop.2021.11.007.

Peppin, S.S. & Style, R.W. (2013) The physics of frost heave and ice-lens growth. *Vadose Zone Journal*, 12(1), 1-12. Available from: https://doi.org/10.2136/vzj2012.0049.

Perrin, J., Cartannaz, C., Noury, G. & Vanoudheusden, E. (2015) A multicriteria approach to karst subsidence hazard mapping supported by weights-of-evidence analysis. *Engineering Geology*, 197, 296-305. Available from: https://doi.org/10.1016/j.enggeo.2015.09.001.

Perrin, J., Noury, G. and Cartannaz, C.2014. ALEA-Karst: vers une approche multicritères de l'aléa lié aux cavités karstiques en contexte calcaire. (In French). Volume 1: Trame méthodologique. Report. BRGM/RP-63771-FR, 77 pages.

Pissart, A. (2000) Remnants of lithalsas of the Hautes Fagnes, Belgium: a summary of present-day knowledge. *Permafrost and Periglacial Processes*, 11(4), 327-355. Available from: https://doi.org/10.1002/1099-1530(200012)11:4%3C327::AID-PPP370%3E3.0.CO;2-Q.

Pissart, A. (2002) Palsas, lithalsas and remnants of these periglacial mounds. A progress report. *Progress in Physical Geography*, 26(4), 605-621. Available from: https://doi.org/10.1191/0309133302pp354ra.

Pissart, A. (2003) The remnants of younger Dryas lithalsas on the Hautes Fagnes plateau in Belgium and elsewhere in the world. *Geomorphology*, 52(1-2), 5-38. Available from: https://doi.org/10.1016/S0169-555X(02)00246-5.

Platel, J.P. & Dubreuilh, J. (1986) Le bassin d'Aquitaine (In French). In: Lajoinie, J.P., Platel, J.P., et al. (Eds.) *Grandes étapes de l'histoire géologique de quatre régions françaises*. Report number: 86-SGN-424-GEO. Paris, France: BRGM. 18 p., 12 fig.

Priestley, S.C., Meredith, K.T., Treble, P.C., Cendón, D.I., Griffiths, A.D., Hollins, S.E., et al. (2020) A 35 ka record of groundwater recharge in south-West Australia using stable water isotopes. *Science of the Total Environment*, 717, 135105. Available from: https://doi.org/10.1016/j.scitotenv.2019.135105.

Prowse, T.D., Wrona, F.J., Reist, J.D., Hobbie, J.E., Lévesque, L.M. & Vincent, W.F. (2006) General features of the arctic relevant to climate change in freshwater ecosystems. *Ambio*, 35(7), 330-338. Available from: https://doi.org/10.1579/0044-7447(2006)35[330:GFOTAR]2.0.CO;2.

Rains, M.C. (2011) Water sources and hydrodynamics of closed-basin depressions, cook inlet region, Alaska. *Wetlands*, 31(2), 377-387. Available from: https://doi.org/10.1007/s13157-011-0147-x.

Saltel, M. (2008). Impact de structures géologiques sur l'alimentation de systèmes aquifères profonds - Fonctionnement hydrogéologique des antiformes du sud de Bordeaux. (In French). [Doctoral dissertation, Université Michel de Montaigne - Bordeaux III]. French repository of doctoral dissertations. https://theses.fr/2008BOR30062

Saltel, M., Rebeix, R., Thomas, B., Franceschi, M., Lavielle, B. & Bertran, P. (2019) Paleoclimate variations and impact on groundwater recharge in multi-layer aquifer systems using a multi-tracer approach (northern Aquitaine basin, France). *Hydrogeology Journal*, 27(4), 1439-1457. Available from: https://doi.org/10.1007/s10040-019-01944-x.

Sassolas-Serrayet, T., Cattin, R. & Ferry, M. (2018) The shape of watersheds. *Nature Communications*, 9(1), 3791. Available from: https://doi.org/10.1038/s41467-018-06210-4.

Sauro, U. (2019) Closed depressions in karst areas. In: White, W.B. et al. (Eds.) *Encyclopedia of caves*, Third edition. Amsterdam: Academic Press, pp. 285-300 https://doi.org/10.1016/B978-0-12-814124-3.00032-7.

Seppälä, M. (2011) Synthesis of studies of palsa formation underlining the importance of local environmental and physical characteristics. *Quaternary Research*, 75(2), 366-370. Available from: https://doi.org/10.1016/j.yqres.2010.09.007.

Sitzia L. (2014) Chronostratigraphie et distribution spatiale des dépôts éoliens quaternaires du Bassin Aquitain. (In French). [Doctoral dissertation, Université de Bordeaux]. French repository of doctoral dissertations https://theses.hal.science/tel-01197192

Sitzia, L., Bertran, P., Bahain, J.J., Bateman, M.D., Hernandez, M., Garon, H., et al. (2015) The quaternary coversands of Southwest France. *Quaternary Science Reviews*, 124, 84-105. Available from: https://doi.org/10.1016/j.quascirev.2015.06.019.

Sitzia, L., Bertran, P., Sima, A., Chery, P., Queffelec, A. & Rousseau, D.D. (2017) Dynamics and sources of last glacial aeolian deposition in Southwest France derived from dune patterns, grain-size gradients and geochemistry, and reconstruction of efficient wind directions. *Quaternary Science Reviews*, 170, 250-268. Available from: https://doi.org/10.1016/j.quascirev.2017.06.029.

Smith H. T. U. (1960) Physiography and photo interpretation of coastal sand dunes. Final report contract NONR-2242 (00), Office of Naval Research, Geographical Branch.

Stadelmaier, K.H., Ludwig, P., Bertran, P., Antoine, P., Shi, X., Lohmann, G., et al. (2021) A new perspective on permafrost boundaries in France during the last glacial maximum. *Climate of the Past*, 17(6), 2559-2576. Available from: https://doi.org/10.5194/cp-17-2559-2021.

Texier, J. (2011) Genèse des lagunes landaises: un point sur la question. (In French). In: Merlet, J. & Bost, J. (Eds.) *De la lagune à l'airial: le peuplement de la Grande Lande*. Pessac, France: Fédération Aquitania, pp. 23-42.

Trask, P.D. (1930) Mechanical analyses of sediments by centrifuge. *Economic Geology*, 25(6), 581-599. Available from: https://doi.org/10.2113/gsecongeo.25.6.581.

Vandenberghe, J., French, H.M., Gorbunov, A., Marchenko, S., Velichko, A.A., Jin, H., et al. (2014) The last permafrost maximum (LPM) map of the northern hemisphere: permafrost extent and mean annual air temperatures, 25-17 ka BP. *Boreas*, 43(3), 652-666. Available from: https://doi.org/10.1111/bor.12070.

Waltham, A.C. & Fookes, P.G. (2003) Engineering classification of karst ground conditions. *Quarterly Journal of Engineering Geology and Hydrogeology*, 36(2), 101-118. Available from: https://doi.org/10.1144/1470-9236/2002-33.

Walvoord, M.A. & Kurylyk, B.L. (2016) Hydrologic impacts of thawing permafrost—a review. *Vadose Zone Journal*, 15(6), vzj2016-01. Available from: https://doi.org/10.2136/vzj2016.01.0010.

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

**How to cite this article:** Bussière, L., Schmutz, M. & Dupuy, A. (2025) A karst rather than periglacial origin for small enclosed depressions of the Landes Triangle, southwest France. *Earth Surface Processes and Landforms*, 50(6), e70091. Available from: https://doi.org/10.1002/esp.70091
Supplementary material S1

List of impervious and karst layers taken into account to assess the vulnerability to karst processes. The * mark indicates that the unit does not belong to the BDLISA list of karst units.

<table><thead><tr><th>Type</th><th>Code BDLISA</th><th>Description</th></tr></thead><tbody><tr><td rowspan="7">Impervious layers</td><td>308AC02</td><td>Brach's clay</td></tr><tr><td>308AC04</td><td>Clay</td></tr><tr><td>312AA01</td><td>Clay</td></tr><tr><td>312AA03</td><td>Loam and pebble clay</td></tr><tr><td>318AA01</td><td>Clay</td></tr><tr><td>322AA01</td><td>Clay</td></tr><tr><td>326AA01</td><td>Deep-sea marl</td></tr><tr><td rowspan="20">Potential karst layers</td><td>324AA01</td><td>Starfish limestone, faluns and sandstone</td></tr><tr><td>346AA03</td><td>Chalk-marl limestone and marl</td></tr><tr><td>346AA05</td><td>Chalk-marl limestone and marl</td></tr><tr><td>348AA01</td><td>Limestone, sandstone, and marl</td></tr><tr><td>348AA02</td><td>Marly limestone</td></tr><tr><td>348AA03</td><td>Limestone, sandstone, and sand</td></tr><tr><td>350AA01</td><td>Limestone and sandstone</td></tr><tr><td>352AA01</td><td>Limestone and dolomite</td></tr><tr><td>356AB01</td><td>Limestone</td></tr><tr><td>358AE03</td><td>Micritic and bioclastic limestone</td></tr><tr><td>358AE07</td><td>Limestone and dolomite</td></tr><tr><td>362AA01</td><td>Limestone and sandstone</td></tr><tr><td>362AA03</td><td>Limestone, dolomite and anhydrite</td></tr><tr><td>362AG01</td><td>Sandstone and dolomite</td></tr><tr><td>366AA01</td><td>Sandstone and dolomite</td></tr><tr><td>320AA01*</td><td>Faluns, sandstone, and limestone</td></tr></tbody></table>
Enclosed depressions in SW France

<figure><img src="image_16.png" /></figure>

**Supplementary material S2:** Results of the geophysical investigation along the survey lines L1 to L7. From left to right: (1) GPR at 250 MHz, (2) GPR at 100 MHz, and (3) ERT models. All the graphs are at the same scale (assuming a wave speed of 0.06 m.ns⁻¹ for the radargrams), with vertical exaggeration by a factor of two. The displayed water tables are those measured specifically on the day of the corresponding survey in well PQ-08515X0053 (250 m away).
<figure><img src="image_17.png" /></figure>

**Supplementary material S3:** Results of the geophysical investigation along the survey lines L1 to L7. From left to right: (1) GPR at 250 MHz, (2) GPR at 100 MHz, and (3) ERT models. All the graphs are at the same scale (assuming a wave speed of 0.06 m.ns⁻¹ for the radargrams), with vertical exaggeration by a factor of two. The displayed water tables are those measured specifically on the day of the corresponding survey in well PQ-08515X0053 (250 m away).