domain: sdv shs

text:
Boris Brasseur
email: boris.brasseur@u-picardie.fr

Théo Allalou

Laurent Chalumeau

Emilie Gallet-Moron

Jérôme Buridant

A chronology of marl pits in the Paris Basin determined from the infilling of closed depressions: do marling practices have their origin in the Neolithic period?

Keywords: Lidar, topographic depression, limestone quarry, calcaric, technosol, soil alkalinisation

Marling (limestone alkaline amendment) agrarian practices have a plurimillennial influence on soil pH and on soil-associated ecosystems. Although the earliest written records in Europe date back to antiquity, the origin of this agrarian practice is not well known. In order to trace the evolution of this practice in the early agrarian societies of Western Europe, we searched for topographic anomalies of ancient limestone quarries that may have been used for marling.

Several hundred circular closed topographic depressions (CCDs) were identified under two areas of ancient forest in northern France with Lidar digital elevation models (DEMs). Soils, morphology and archaeological context were studied and filling sediments were dated in several dozen CCDs.

The gaps in the carbonate rocks, soils with calcaric properties (scattered limestone fragments), unstratified fills, the absence of buried paleosols, and dates (OSL, C-14) later than the first local agricultural societies made it possible to recognise the CCDs that may have served as limestone quarries. This was the case for most of the 30 CCDs studied, leaving a minority of natural origin (suffosion, depressions of thermokarstic origins). Amongst the quarries

Introduction

Many plant species in food crops require near-neutral soil pH conditions for optimal growth and yields. Marling or liming of soils is practised in agriculture to neutralise soil acidity and improve nutrient availability [START_REF] Läuchli | Plant stress under non-optimal soil pH[END_REF][START_REF] Matthew | Marling in British agriculture: a case of partial identity[END_REF]. This practice is particularly relevant in areas where soils are subject to natural acidification when rainfall exceeds evapotranspiration [START_REF] Slessarev | Water balance creates a threshold in soil pH at the global scale[END_REF]. The most common materials used for marling are various limestones in granular (friable or fragmented) form and their derivatives such as lime [START_REF] Goulding | Soil acidification and the importance of liming agricultural soils with particular reference to the United Kingdom[END_REF]. However, in the absence of these materials, farmers may also use sandy deposits of calcareous seaweed [START_REF] Clout | Fertilisants minéraux en France au XIXe siècle[END_REF] or the alkaline plant ash [START_REF] Juo | Chemical dynamics in slash-and-burn agriculture[END_REF] during slash and burn practices. For all soils undergoing acidolysis and acido-complexolysis processes, the practice of marling will constitute a revolution in the dynamics of functioning processes. Due to the increase of pH and the This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d loess quarries. The authors seriously considered the hypothesis that the quarries used to serve for calcareous/alkaline amendments for crop improvement.

In the Merdaal forest (central Belgium), [START_REF] Vanwalleghem | Origin and evolution of closed depressions in central Belgium, European loess belt[END_REF] dated the colluvium infill of a CD between the first century BC and the fourth century AD. The CD was contemporary with local agricultural activities, inferred from gully infills dated from the 18th century BC to the third century AD. The authors attributed the excavation of the calcareous loess in the local CDs to the practice of marling. Near the same forest, [START_REF] Rommens | Reconstruction of late-Holocene slope and dry valley sediment dynamics in a Belgian loess environment[END_REF] dated the colluvium infill of one of these anthropogenic CDs to the Middle Bronze Age (1120-900 BC). [START_REF] Meylemans | The evaluation of archaeological sites using LIDAR and erosion/sedimentation modelling[END_REF] also reported similar CDs, without dating the infill, in the vicinity of the Middle Neolithic site of Ottenburg (Central Belgium).

In western Germany [START_REF] Baumewerd-Schmidt | Von Restfundstellen und Scheinfundstellen-Ergebnisse einer Grabenbetreuung in der Lösslandschaft[END_REF] and in Britain [START_REF] Prince | Some reflections on the origin of hollows in Norfolk compared with those in the Paris region[END_REF], other CDs have been recognised as marl pits, but no further information on their age is available.

Other CDs of anthropogenic origin had multiple functions compatible with a possible initial use as marling. Thus, in the North-East of France [START_REF] Etienne | The origin of closed depressions in Northeastern France: A new assessment[END_REF], CDs on marl were probably later used for watering herds and/or for retting flax. The filling of the CDs studied dates from the Iron Age to the Roman period. Some of these CDs could also be of natural origin following the dissolution of diapirs or the subsidence of the soil. These findings were complemented by other work in Luxembourg [START_REF] Slotboom | Pollen records of mardel deposits: The effects of climatic oscillations and land management on soil erosion in Gutland, Luxembourg[END_REF][START_REF] Van Mourik | Natural versus anthropogenic genesis of mardels (closed depressions) on the Gutland plateau (Luxembourg); archaeometrical and palynological evidence of Roman clay excavation from mardels[END_REF], where an anthropogenic origin was found in some cases, with possible clay quarries for ceramics. The age of the fills in these quarries ranges from Roman times to the Middle Ages. The absence of buried paleosols under/inside the colluvium infill is the sign of a CD of anthropogenic origin. [START_REF] Kołodyńska ; Gawrysiak | Closed depressions in the European loess belt -Natural or anthropogenic origin?[END_REF] carried out an overview of the work on CDs in the loess belt in Europe. The analysis of Eastern European CDs mainly suggests natural landforms (from loess deflation, dissolution, thermokarst processes). According to 14 C datings, the CDs in eastern Poland were formed between the Late Glacial and the mid-Holocene period. The morphometric characteristics, ages and

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d stratigraphy of the infills clearly indicate that natural thermokarstic depressions were filled by colluvial deposits due to agricultural land use in the second half of the Holocene. The presence of mid-Holocene mature paleosols (Chernozems, Phaeozems, Luvisols or Albeluvisols) buried under the colluvium infilling also points to an ancient natural origin of the CDs (Kołodyńska-Gawrysiak, 2019). In Western Europe, a dominant anthropogenic origin (i.e. loess and limestone quarries, underground quarry collapses, bomb and mine craters) is more common for CDs. But natural landforms such as sinkholes, salt lenses and limestone dissolution dolines have also been reported [START_REF] Delafosse | De l'origine des mardelles de Lorraine[END_REF], especially in the Paris region [START_REF] Pissart | Les dépressions fermées de la région parisienne. Le problème de leur origine[END_REF].

To sum up, the similar surface forms resulting from this wide variety of CDs constitutes a possible confusing factor in the interpretation of these structures. Systematic investigation of the fillings and their dating is therefore necessary. Cereals, originating from the Near East, are cultivars adapted to local neutral-alkaline soils and must therefore have required soils with neutralised acidity from the beginning. Our hypothesis is that the practice of marling is much older than the antiquity period because the first farmers in Northwest Europe had to adapt the chemical qualities of naturally acidic soils to the plants they cultivated. We believe that in former agrarian territories the anthropic impact of marling on the soils is often underestimated and that this practice has strongly and durably structured the characteristics of soilscapes where it has developed.

Our objective was to date a set of filled marl pits in the north-central part of the Paris basin in order to estimate the age of the soil marling practices. To identify a maximum number of marl pits, we targeted i) former agrarian territories where CDs are abundant, ii) depressions under forest cover where the ancient agrarian microtopography is better preserved [START_REF] Lenoir | Unveil the unseen: Using LiDAR to capture time-lag dynamics in the herbaceous layer of European temperate forests[END_REF], and iii) areas where the subsoil would have been favourable for marling.

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d In order to identify the origin of the CDs, we studied their infilling (stratigraphy), their pedological contexts and the contact with the local limestone substratum. To date the selected CDs, we decided to analyse the lowermost infillings with Optically Stimulated Luminescence (OSL) and carbon 14 ( 14 C) when rare charcoal remains were discovered. We selected two large woodland areas, Compiègne and Saint-Gobain forest massifs, exhibiting a high diversity of soil types over contrasted geological substrates, making available the targeted limestone properties. The geological strata in these forests are relatively thin and have a lateral variability in the facies that makes it easier to discern the strategies used by former inhabitants for the acquisition of carbonate rocks and the characteristics they sought. Both forests were agricultural areas during the Roman period and probably even a few centuries before. High-resolution lidar DEMs enabled us to locate and then determine the geomorphological and fossil agrarian context of the observed CDs.

Finally, it should be noted that karstic networks are very rare in the limestone areas in the northern Paris Basin [START_REF] Rodet | Les karsts de la craie et du Cénozoïque du bassin de Paris, in: Grottes et Karsts de France, Karstologia-Mémoires[END_REF], which limits the risk of natural sinkholes forming on the surface.

Material and methods

Study sites, past and present contexts

The forests of Compiègne and Saint-Gobain are located in the north-central part of the Paris basin (Fig. 1). The first farming populations settled around 4500 BC in the Oise valley bordering the current territory of the two forests [START_REF] Leroyer | L'IMPACT DES GROUPES NÉOLITHIQUES DU BASSIN PARISIEN sur le milieu végétal: évolution et approche territoriale d'après les données polliniques[END_REF][START_REF] Leroyer | Nouvelles données sur le paysage végétalde la moyenne vallée de l'Oise au subboréal :l'étude de la séquence organique de Thourotte(Oise, France). Quaternaire[END_REF]. Field surveys, complemented by a remote sensing survey with Lidar technology [START_REF] David | Rapport de synthèse des prospections archéologiques d'après les données lidar[END_REF], carried out on the two forests revealed hundreds of Roman sites, mainly rural settlements, within the forests [START_REF] Pichon | Carte archéologique de la Gaule 02 : L'Aisne, Carte Archéologique de la Gaule[END_REF]Thuillier, 2017).

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d The Gallo-Roman settlements were abandoned rather abruptly towards the end of the 4 th century AD, with archaeological evidence of destruction and fires as early as the 3 rd century AD (Fremont andWoimant, 1976, 1975;[START_REF] Tuffreau-Libre | La céramique commune gallo-romaine de la forêt de Compiègne (Oise) au Musée des Antiquités nationales[END_REF]. On the territory of the Compiègne forest, a woodland is mentioned for the first time in the 6 th century by Grégoire de Tours (de Tours, 1836); the forest then became a royal hunting reserve until the French revolution in 1789 AD. The Saint-Gobain forest was part of the Voas forest in the early Middle Ages (Acta Sanctorum, Vita Sancti Gobani). The study of the This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d local Lidar DEM reveals numerous Gallo-roman agricultural parcels. This forest therefore probably developed between the final period of the Western Roman Empire and the early Middle Ages. Today, the local climate is sub-oceanic with an average annual temperature of 10.5°C and an annual rainfall of 650 mm (Climate-data.org, n.d.). Acidolysis and acidocomplexolysis are the two main geochemical soil-forming processes. The geological substratum is mainly composed of marine tertiary tabular sedimentary rocks (Thanetian to Barthonian), sometimes covered by more recent continental deposits (Pleistocene to Holocene). In both forests, the three dominant soil groups are Cambisol, Luvisol and Podzol. In the 20th century, a 1200-ha plateau area south of the Compiegne forest, mostly covered with fields, was supplied with limestone by three marl pits [START_REF] Guérin | Carte Géologique de la France à 1/50 000[END_REF]. These pits exploited Middle-Upper Lutetian limestone and were all closed towards the end of the 20th century. This recent situation indicates a rarefication of marl pits over two centuries [START_REF] Clout | Fertilisants minéraux en France au XIXe siècle[END_REF][START_REF] Matthew | Marling in British agriculture: a case of partial identity[END_REF].

Field survey, laboratory analyses and data analyses

Around 150 and 250 circular closed depressions (CCDs) have been characterized within Compiègne (14,357 ha) and Saint-Gobain (8,470 ha) forests, respectively, thanks to airborne Lidar DEMs (hillshade models with two mutually perpendicular directions). The realisation of the shaded DEM and the calculation of depth, area and diameter of CCDs were carried out with the ArcGis Pro 2.8 software (ESRI).

In the Compiègne forest (FDC), among the 150 CCDs recorded on the DEM, a high density is observed on the southern plateau (Fig. 1). Finally, only the central and eastern parts of the forest, with upper Pleistocene blown sand formations, have few CCDs. It is also on these siliceous sand dunes that the absence of agrarian parcelling and other archaeological structures had previously been noted [START_REF] Horen | Relation entre les structures archéologiques identifiées sur un levé LIDAR et la typologie des sols du massif forestier de Compiègne (Nord de la France)[END_REF]. In the Saint-Gobain massif (SGB), where such thick dune formations are absent, the 250 CCDs are more equally distributed.

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d We selected CCDs that were sub-circular, over 10m in diameter, outside the areas of ancient settlements and did not include bomb holes (characterised by smaller, rounder and sharper holes). We then studied the infillings of 30 of these CCDs in various geological contexts, ranging from plains (minimum altitude 48m) to plateaus (maximum altitude 170m).

For each CCD, hand auger cores were taken both inside and outside the depression i) until the first hard rock is encountered irrespective of the depth; ii) to 1.5-metre-deep inside in situ carbonate soft rocks; or iii) to 6m deep in other cases. Outside the CCDs, we determined the "local reference" for soil and substrate rocks through (at least) one auger survey. Within the CCDs, two boreholes (auger and/or soil pit) were drilled or dug to characterise the infills. We removed the CCDs with no carbonate rocks outside and inside the depression in the first five metres and those where local carbonated rocks appeared to be intact in the CCD footprint. We assumed the remaining 25 CCDs were potential marl pits and therefore retained them for further investigation.

Based on this survey and with the help of the DEM, several morphological characteristics were collected for the 25 CCDs: average elevation at the centre and outside of the CCDs; average diameter and circularity (Perimeter^2 / (4 * pi * Area)); thickness of the infill in the CCDs; thickness of the limestone gap; minimum depth of the carbonate fragments in the fill; and relative infill ratio (Infill thickness/(Delta topographic + Infill thickness)). The infill thickness of the CCDs was defined as the depth required to reach the unreworked carbonate rock at the centre of the depression. Outside the CCDs, we also measured the depth of the local soil and the thickness of any underlying soft limestone. Seventeen were chosen for dating purposes; a trench was dug inside these CCDs to 1.5/2.2m in depth or to the first hard rock substrate. A description of the pedogenetic, sedimentary and stratigraphic characteristics of the CCD infilling was then carried out on the soil profiles obtained. Soil horizons and sediment units were characterised in the field (based on colour, texture, structure, spatial relationships) and following the World Reference Base for Soil Resources (IUSS Working Group WRB,

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d 2015). In order to obtain a complete nomenclature of the horizons adapted to the description of the diversity of the underlying source rocks, we also described the sequences of horizons according to the French pedological reference system (AFES 2009). We also carried out pH analyses on the soils (5, 25 and 50cm depth) inside and outside the CCDs and measured CaCO 3 levels with a Bernard calcimeter (volumetric method) on the same soil samples and on the calcareous substrates in the CCDs. Soil pH was measured with two glass electrodes in distilled water (ISO 10390:2004) on two subsamples, then averaged. We performed a Wilcoxon pairwise comparison rank test to evaluate differences in pH water values in the soils inside and outside the CCDs soils. All statistical analyses were performed within the R software environment (R Development Core Team, 2016), with the "ggplot2", "psych" and "corrplot" packages.

Of the 17 CCDs, 16 were dated by OSL (n=14) or/and 14 C (n=4). OSL dates were established by the Luminescence Dating Laboratory of the Silesian University of Technology, Gliwice. Multigrain measurement of quartz grains (125-200µm) was realised following the single aliquot regenerative method for 16 aliquots and a "Central Age Model" was established [START_REF] Galbraith | Optical Dating of Single and Multiple Grains of Quartz from Jinmium Rock Shelter, Northern Australia: Part I, Experimental Design and Statistical Models*[END_REF]. 14 C dating was carried out following the AMS method at the Poznań Radiocarbon Laboratory. The calibration was done on the OxCal software v4.4.2 [START_REF] Ramsey | Bayesian Analysis of Radiocarbon Dates[END_REF] and with atmospheric data from Reimer et al [START_REF] Reimer | The IntCal20 Northern Hemisphere Radiocarbon Age Calibration Curve (0-55 cal kBP)[END_REF].

Results

Geological context & carbonate rocks gaps

Within the 25 CCDs, gaps (compared to their surroundings) were observed for several carbonate rock facies: loess (2), Lutetian limestone (15), Ypresian coarse limestone and "faluns" (4), and Thanetian limestone and marl (4) (Fig. 2). Loess has a relatively low carbonate content (8-15%) whereas Lutetian This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d limestones can reach 97% CaCO 3 . Of the 25 CCDs, 18 were dug at the expense of soft (granular to friable) limestone material. In the case of the CCDs present on the Lutetian limestone, 9 CCDs (FDC-CCD1, 2, 3, 4, 5, 9, 10, 12; SGB-CCD4, CCD8) have a flat bottom (sometimes with a slight slope from the centre to the edges) composed of massive hard limestone bedrock. Fig. 2: Stratigraphic log of the global sequence of Tertiary rock found in the Compiègne and Saint-Gobain forests. White rocks with brick-like features correspond to calcareous rocks. The 25 CCDs with limestone gap were counted for these different limestone facies on the log. Thickness data and descriptions are taken from the 1:50,000 geological maps of "Attichy" and "La Fère" (BRGM). The % of CaCO3 corresponds to the calcium carbonate content we measured with the Bernard calcimeter or was taken from published data in the local geological map records.

Closed depressions: morphology and relative infilling

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d On the Lidar DEM, the 25 candidate CCDs for marls pits have sub-circular shapes. Their circularity index ranges from 1.02 to 1.19. Regarding their average diameters, they vary from 55.5m to 28.5m and their topographic depths (relative to the average local ground elevation) from -320cm to -60cm. No positive relief of extracted soil is associated with these CCDs, so if it was extracted, the exceeding earth appears to have been spread further beyond. The thicknesses of the calcareous rocks gaps range from 65 to 335cm. Finally, the thickness of the sedimentary fill in the CCDs varies from 360 to 25cm and their relative infill ratio from 0.09 to 0.82. Analysis of the morphological features of the CCDs reveals that the sharpest and deepest topographic depressions are those with the lowest infilling ratio (e.g. FDC-CCD1 to CCD4, Fig. 3). Conversely, the shallow marked depressions have a higher infilling ratio. Thus, two opposing morphotypes are distinguished, one with a well-defined recessed relief and a relatively thin infilling, while the other has a faint relief but a much thicker filling. FDC-CCD1 and CCD3 illustrate these two opposite morphotypes (supplementary data, Fig.

A.1).

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d Fig. 3. Top: The topographic profiles of the four CCDs with red points illustrate the diversity of landforms encountered. Bottom: Graph representing the distribution of CCDs as a function of two variables, on the X-axis the difference in altitude (in centimetres) on the Lidar DEM between the centre of the CCD and the outside, and on the Y-axis the Relative infill ratio. The projected regression line with confidence intervals (shaded area) shows a negative relationship between these two variables.

The correlation matrix below (Fig. 4) also shows a positive correlation between the CaCO 3 thickness deficit and the topographic depth of the CCDs and a negative correlation between the topographic depth and the relative infill ratio. We can also observe that the CCDs will be more filled (Relat_infill_ratio & Infill_tchick) if the external soils are thick. Finally, these soils will become thinner if a significant thickness of carbonate rock has been extracted (CaCO3_gap_thick). Among the 25 marl pits, one of them (FDC-CCD1) had almost no fill, which allowed us to use it as an example to describe the morphology of the contact between bedrock and fill for the entire CCD. FDC-CCD1 is a relatively deep CCD (250cm) with a particularly shallow fill (20-30cm) and steep sloping edges. The low fill revealed the general shape of the contact with the un-removed soil and rock. In this

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d CCD, we therefore noticed a lower slope on one sector of the southern edge. A 80 cm-deep soil pit was made at this point. This revealed a steep boundary with a dense level of limestone fragments associated with few pieces of reworked soils and some charcoal under the fill composed of soil and few limestone fragments. The high degree of compaction of this level strongly resembles a trampled surface. This surface follows a gentle slope from the outside to the inside of the depression.

3.3.Infilling sequences: description and dating

We described 17 CCDs infillings by assessing excavated trench walls and a complementary augering survey (supplementary data, Fig . A.3). Particular attention was paid to breakage phenomena in the filling sequences of the CCDs, whether they were stratifications, the formation of a soil, or cross-cutting by secondary excavations. Among the 17 CCDs, only one (FDC-CCD11) had a buried soil horizon, indicating a stable pedogenesis step inside the total infilling sequence. At the base of the infilling, we found a 40cm-thick buried leptic calcaric pretic Anthrosol. We retrieved some roman archaeological material into this thin buried soil layer (ceramic plate and earthenware fragments, charcoal, flint chips). In the case of the CCDs SGB-CCD1 and SGB-CCD11, a sub-horizontal discontinuity intersected the deposits of a previously formed fill. This could indicate a local secondary excavation, possibly related to the installation of a charcoal kiln platform during recent forestry history. With the exception of FDC-CCD11, the absence of an intercalated pedogenetic layer in 15 out of 16 infill sequences indicates single-phase events. The charcoal and datable archaeological artefacts (ceramics) found in the pits all predate the infilling event. However, when these elements are present close to the soil surface (<100cm), this may indicate more recent contamination by bioturbation (biogaleries, windfalls). This was the case for four CCDs (FDC-CCD6 and CCD8, SGB-CCD1 and CCD8) where kilns had been installed on the soil surface. Some pieces of charcoal had reached depths of nearly one metre from much more abundant accumulations near the surface (<25cm). In such situations we did not consider these coals for the dating of the CCDs fill.

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d Optically stimulated luminescence thermochronometry (OSL) tube sampling was carried out towards the base of the depression fill. Multi-grain quartz dating of the samples yielded ages between 28.9 and 1.793 kyrs BP (Table 1). Six dates were between 28.900 and 12.710 kyrs BP; that is towards the end of the Late Pleistocene. The remaining eight dates were in the second half of the Holocene, from 5.850 to 1.793 kyrs BP. Where charcoal was present in a substantial number at the base of the depression fill, a sample was taken and analysed for 14 C dating. Four calibrated age intervals were obtained, with medians of 4306, 2500, 2128 and 1936 yrs BP (Table 2). All datings were within the second half of the Holocene.

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d FDC-CCD3 was dated by both OSL and 14 C, yielding 2400 yrs BP and 2500 cal. yrs BP respectively. The two methods therefore provided very consistent ages for the fill in this depression. The fill in the FDC-CCD17 depression was also doubly dated. The OSL age obtained was 5380 BP while the 14 C age was 4306 cal. yrs BP. The difference is greater but the two dates are still relatively close if we consider the distribution of all the dates obtained. 19,20,21;11) had a fill that predated the early Holocene which excludes them from being limestone quarries contemporary with agricultural societies.

3.4.Pedogenetic development inside and outside the CCDs

The most common soil type found inside the 25 CCDs (Table 3) was Calcaric Cambisol (12/25) followed by a diverse set of Luvisols (11/25). A few of the soils (FDC-CCD14/ CCD19/ CCD20/ CCD21, SGB-CCD1/ CCD5/ CCD11) had no calcareous characteristics. Of the latter, all but one (SGB-CCD1) had pre-Holocene infills. In the case of SGB-CCD1, the thick surrounding local decarbonated silt soils could explain the absence of carbonates in the fill. For Luvisols, the base of the BT horizon often coincided with fill levels that were still rich in limestone fragments. Outside the CCDs, the most common soil type was Luvisol (16/25), followed by Calcaric Cambisol (9/25). In several cases (FDC-CCD12 & CCD17, SGB-CCD1), the Luvisols encountered outside the CCDs had sharper and more clayey subsurface soil horizons (Bt).

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d Regardless of depth (5, 25, 50cm) Wilcoxon rank tests returned non-significant differences in pH water values between inside and outside the CCDs (Table 4). -50

6.29 6.40 0.925 Several CCDs (FDC-CCD11, CCD12, CCD16, SGB-CCD1) also contained pieces of ancient A, S, E and BT horizons mixed with limestone fragments in the lower part of their fill. During the digging of the soil pits, these pieces of ancient horizons were distinguishable by clear boundaries between blocks of different colours and textures. They disappeared progressively in the upper part of the fill, particularly in the 50 cm just below the surface. The frequency of bioturbation in the superficial part of the presentday soils could explain the disappearance of these very friable fragments with fine matrices through homogenisation. After this homogenisation, only the limestone fragments, if they have not been dissolved by the surface acidification of the soil, still bear witness to the diversity of origins of the infill materials.

Discussion

4.1.A minority of natural depressions with pre-Holocene infill and a majority of limestone quarry (potential marl pits) The fill of six CCDs dated back from Late Glacial Maximum (LGM) to Younger Dryas (Fig. 5) match the age of aeolian sandy (FDC-CCD19, CCD20, CCD21, SGB-CCD5) or silty (FDC-CCD14, SGB-CCD11) deposits and miss artefacts from Holocene agricultural societies as well as scattered carbonate rock fragments.

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d These depressions are likely natural. Several hypotheses can explain how they trapped sediments during the last Ice Age, but two are more plausible, namely thermokarst depressions in a context of wind-blown sands and loessic silts 19,20,; and suffosion due to the strong irregularities in the depth of the limestone/fill interface (for SGB-CCD5, located on a limestone plateau). Alternative hypotheses include various natural processes leading to CDD formation in the European loess belt (see for example [START_REF] Kołodyńska ; Gawrysiak | Closed depressions in the European loess belt -Natural or anthropogenic origin?[END_REF], like. piping erosion, irregularities in the deposition of aeolian sediments, deflation hollows [START_REF] Gillijns | On the characteristics and origin of closed depressions in loess-derived soils in Europe-a case study from central Belgium[END_REF]. According to the age of infilling, 10/16 of CCDs date back to the second half of the Holocene (Fig. 5) and correspond to potential marl pits. They span time periods from the Middle Neolithic, to the Bronze and Iron Ages, to the Roman antiquity. This is consistent with the time at which both forests started to establish on these former agricultural lands, i.e. the end of antiquity (de Tours, 1836; Fremont andWoimant, 1976, 1975;[START_REF] Pichon | Carte archéologique de la Gaule 02 : L'Aisne, Carte Archéologique de la Gaule[END_REF][START_REF] Thuilliers | Dynamique des paysages de l'Antiquité au Moyen Age[END_REF][START_REF] Tuffreau-Libre | La céramique commune gallo-romaine de la forêt de Compiègne (Oise) au Musée des Antiquités nationales[END_REF]Acta Sanctorum, Vita Sancti Gobani). The earliest potential marl pit we found date back to the Middle Neolithic, hence it is 2000 years older than marl pits dated by [START_REF] Rommens | Reconstruction of late-Holocene slope and dry valley sediment dynamics in a Belgian loess environment[END_REF] from the Middle Bronze Age.

4.2.Limestone quarry infilling and pedogenesis

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d The nine potential marl pits showing a single-phase filling suggests intentional filling immediately after reclamation of the quarry, plausibly with the aim of rehabilitating it into agricultural lands. A similar observation was reported by [START_REF] Gillijns | On the characteristics and origin of closed depressions in loess-derived soils in Europe-a case study from central Belgium[END_REF] and [START_REF] Poesen | Gullies and Closed Depressions in the Loess Belt: Scars of Human-Environment Interactions[END_REF], who identified carbonate rock quarries filled that had been more or less levelled depending on their location in the field or under forest cover in central Belgium. However, in their case the filling was not single-phased but composed of several successive layers of colluvium. Our sole exception was site FDC-CCD11 which showed a paleosol at the basis on a thin, likely natural fill. It indicates a delay between the abandonment of the quarry and the infilling operation.

In contrast, FDC-CCD1 is an abandoned quarry with no previous infilling; therefore, the thin existing fill is probably natural. This CCD allows to observe the natural functional morphology (bowl shape, steep slopes) and provides a model for the other local CCDs. The same CDD showed an well-preserved edge with a levelled slope and a buried compacted surface which might have been an access ramp to the bottom of the quarry. It is thus likely that this quarry has been abandoned at the same time as agricultural activities, before the land became afforested (4th to 5th century AD). Carbonate rock fragments dispersed in a matrix of earth were retrieved in the lower half of the filling of 0/6 natural CCDs and 9/10 of potential marl pits, consistent with their hypothesized use as quarries. Among The notable exception is site SGB-CCD1 which contained loess with 10-15% calcium carbonate, where soil acidification may have more rapidly eliminated the remnants of this fine and weak carbonate component in the filling. The absence of carbonate rock fragments near the surface can be due to soil acidification (from the surface) after the afforestation of these territories [START_REF] Brasseur | What deep-soil profiles can teach us on deep-time pH dynamics after land use change? Land Degradation & Development[END_REF]. Applying the criteria of presence of carbonate rock fragments and artefacts of ancient agrarian societies in the fill, there are out of 25 CCDs with a limestone gap (on the inner side), 19 quarries (& marl pit candidates). We conclude that in the northern part of the Paris Basin nearly 75% of the CCDs associated with a carbonate gap are former carbonate rock quarries now filled with a mixture of carbonate materials.

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d Without the help of a DTM Lidar, these advanced Technosols, with frequently Calcaric, Leptic, Cambic or Luvic diagnostic horizons can easily be considered as calcaric cambisols or a rendzic leptosols at first glance. Considering the long history of local agriculture and anthropogenic erosion, the soils surrounding these CCDs can often be classified as Anthrosols.

The thicker and more calcareous soils inside former quarries compared to soils outside (Table A.2), as well as the less pronounced horizonation (cf. rarity of clay (BT) and albic (E) horizons) can be explained by the relatively young age of the parent material and thus to less advanced depletion/enrichment processes.

In the absence of ancient marling activities and erosion due to agriculture, the soils encountered in the areas studied would have been Albic Stagnic Fragic Luvisols and Dystric Cambisols depending on the depth of the limestone susbstrate. Such soils are now only found at the bottom of the slope, buried under the thick colluvium of the first phases of cultivation. In areas covered by quartz drifted sands, Albic Podzosols are expected to dominate. Where these sands reach significant thicknesses and in the absence of nearby calcareous substrates the dunes formed have sometimes been spared by agriculture and these soils are then preserved [START_REF] Horen | Relation entre les structures archéologiques identifiées sur un levé LIDAR et la typologie des sols du massif forestier de Compiègne (Nord de la France)[END_REF].

4.3.Materials suitable for marling extracted from quarries and archeogeographic distribution

Where several carbonate rock types were locally accessible, as in our study area, extraction was prioritary done on the most friable carbonate layers with low clay content. Once those extracted, extraction often stopped on a hard limestone (FDC-CCD1, CCD2, CCD3, CCD4, CCD5, CCD9, CCD10, CCD11, CCD12, CCD16, CCD17; SGB-CCD4, CCD8) or a marl (FDC-CCD7, CCD13). This further supports the hypothesis that these former quarries were used for carbonate amendment in surrounding agricultural lands, though secondary uses cannot be fully ruled out (e.g. paving of roads and paths, mortar). Compared to the situations observed in central Belgium [START_REF] Poesen | Gullies and Closed Depressions in the Loess Belt: Scars of Human-Environment Interactions[END_REF], carbonate loess was rarely used in our study area, probably because more suitable and/or more accessible substrates (e.g. CaCO3-rich limestone rocks below a thin layer of silt) were accessible in our case.

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d We therefore favour the hypothesis that these depressions are indeed old marl pits rather than carbonate quarries of undetermined use. The spatial layout of the most recent marl pits supports this hypothesis since they were clearly integrated into the agrarian plots (see appendix Fig. A4). Another striking feature is the high density of quarries in a limited area. On the Compiègne forest southern plateau, each unit of land (from 5 to 13 ha) had its own marl pit. This illustrates the importance in the past of avoiding the cost of transporting large volumes of the local limestone. Marl pits were therefore much more abundant features in earlier agrarian landscapes.

How old is the alkaline revolution in the north-western European soilscapes?

Here we provided robust diagnostic criteria to distinguish former carbonate rocks quarries among lidar-detected CDDs: i) the presence of carbonate rocks fragments in the fill matrix; ii) dating to the second half of the Holocene; and iii) the presence of other artefacts associated with agricultural societies.

Compared to other continents, Europe has the highest proportion of carbonate-containing rocks near the surface [START_REF] Hartmann | The new global lithological map database GLiM: A representation of rock properties at the Earth surface[END_REF]. Most of North-West Europe agricultural soils are currently affected by marling or more modern alkaline amendments, which has an enormous influence on their physico-chemical and biological processes [START_REF] Goulding | Soil acidification and the importance of liming agricultural soils with particular reference to the United Kingdom[END_REF][START_REF] Griffiths | Soil biodiversity and ecosystem functions across Europe: A transect covering variations in bio-geographical zones[END_REF][START_REF] Jauzein | Flore des champs cultivés[END_REF][START_REF] Turbé | Soil biodiversity: functions, threats and tools for policy makers[END_REF][START_REF] Winiwarter | From marl to rock powder: On the history of soil fertility management by rock materials[END_REF]. This widespread practice is much more deeply rooted in our agricultural history than previously thought. It has also contributed to structure ancient agrarian landscapes for several thousand years. The legacy on pedological processes and edaphic ecosystems is still very important [START_REF] Brasseur | What deep-soil profiles can teach us on deep-time pH dynamics after land use change? Land Degradation & Development[END_REF][START_REF] Dambrine | Present Forest Biodiversity Patterns in France Related to Former Roman Agriculture[END_REF], even on territories considered to be little impacted, such as circa. 1.6 millenia old forests [START_REF] Dupouey | Irreversible impact of past land use on forest soils and biodiversity[END_REF]. More studies are needed to better measure the temporal depth of this anthropic influence in order to better model the long-term influence on the evolutionary trajectories of soils and associated ecosystems [START_REF] Ewald | The calcareous riddle: Why are there so many calciphilous species in the Central European flora?[END_REF]. Still our results support the hypothesis of an "alkaline soil revolution" at the Neolithic period, which primarily affected siliceous soils sensitive to

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d natural acidification due to humid climate conditions [START_REF] Slessarev | Water balance creates a threshold in soil pH at the global scale[END_REF] when a limestone substrate was present nearby. It is likely that this alkaline revolution deviated these acidic soils durably from their natural evolutionary trajectory.

Fig. 1 .

 1 Fig. 1. Sitemap of the Saint-Gobain and Compiègne forests in the north-central part of the Paris

Fig. 4 .

 4 Fig. 4. Six morphological features of the 25 CCDs are compared in this correlation matrix with Spearman

Fig. 5 .

 5 Fig. 5. OSL (before 1950) and Carbon-14 (cal. BP) dating results of the CCD infill integrated into a (non-

Table 1

 1 Dating results of CCD basal infillings. OSL ages stated are expressed in years before 1950.

	Forest CCD	epth	Material	Th-232	228U	40K	Total	Equivalent	Age OSL (Kyrs
		n°D (m)		(Bq/kg)	(Bq/kg)	(Bq/kg)	Dose	dose DE	Before 1950)
								rate	(Gy)
								(Gy/kyrs)
	FDC	CCD3	1.5	Quartz	27	30.1	301	2.016	4.99	2.400(+/-0.013)
				very fine
				sand
	FDC	CCD4	1.1	Quartz	31.5	45.2	288	2.24	4.18	1.793(+/-0.086)
				very fine
				sand
	FDC	CCD9	1.2	Quartz	18.3	27.9	215	1.599	5.08	3.110(+/-0.019)
				very fine
				sand
	FDC	CCD12 0.95	Quartz	13	29.5	199	1.455	3.269	2.180(+/-0.012)
				very fine
				sand
	FDC	CCD13 1.2	Quartz	10.5	22.88	96.8	1.007	2.86	2.770(+/-0.014)
				fine sand
	FDC	CCD14 1.3	Quartz	45.2	32.8	471	2.75	62.3	22.600(+/-0.012)
				fine sand
	FDC	CCD17 0.7	Quartz	38	26.5	372	2.283	12.44	5.380(+/-0.032)
				fine sand
	FDC	CCD19 0.7	Quartz	10.9	11.3	321	1.47	34.3	23.300(+/-0.018)
				fine sand
	FDC	CCD20 0.7	Quartz	11	13.9	93.9	0.878	11.29	12.790(+/-0.074)
				fine sand
	FDC	CCD21 0.75	Quartz	8.22	7.09	165	0.919	13.84	14.990(+/-0.091)
				fine sand
	SGB	CCD1	1.8	Quartz	43.5	43	532	3.05	8.87	2.840(+/-0.013)
				fine sand
	SGB	CCD5	1.1	Quartz	27.26	22.57	239	1.661	21.23	12.710(+/-0.073)
				fine sand
	SGB	CCD8	1.1	Quartz	26.8	23.22	242	1.673	9.91	5.850(+/-0.034)
				fine sand
	SGB	CCD11 1.7	Quartz	29.31	25.12	350	2.025	58.6	28.900(+/-0.014)
				fine sand

Table 2

 2 Results of 14 C dating for charcoal from CCD infillings. The calibrated 14 C dates given are intervals of calendar years, encompassing the true age of the samples with a probability of ca. 95%. The calibration was done with the OxCal software(OxCal v4.4.2 Bronk Ramsey (2009); r:5 Atmospheric data from[START_REF] Reimer | The IntCal20 Northern Hemisphere Radiocarbon Age Calibration Curve (0-55 cal kBP)[END_REF]). CCDs likely to be former limestone quarries, 16 have been dated. Of the 16 CCDs dated, six

	Forest CCD	epth(m) Material 14 C Age	Calibrated age (OxCal 4.4)	Median age cal. BP
		n°D		BP	Probability intervals
	FDC	CCD3	1.4	Charcoal	2440 ±35	594BC (59.6%) 407BC	2500
	FDC	CCD11 1.9	Charcoal	1985 ±30	44BC (86.4%) 85AD	1936
	FDC	CCD16 0.9	Charcoal	2135 ±30	206BC (81.2%) 51BC	2128
	FDC	CCD17 0.9	Charcoal	3845 ±30	2410BC (86.1%) 2203BC	4306
	Of the 25

Table 3

 3 Classification (IUSS Working Group WRB, 2015) of soils and underlying rocks, inside and outside the 25 surveyed CCDs (candidate CCDs for marls pits). Light orange indicates CCDs without the diagnostic features of a limestone quarry.

	Forest	N°	Soil (& Horizons) inside CCD	Soil (& Horizons) outside CCD
		CCD
	FDC	1	Rendzic Leptosol (Aca, Rca)	Endoleptic Luvisol (A, E, BT, Mca, Rca)
	FDC	2	Luvic Endocalcaric Cambisol (A, BT, Cca,	Endoleptic Luvisol (A, E, BT, Mca, Rca)
			Rca)
	FDC	3	Luvic Calcaric Cambisol (A, S/BT, Cca, Rca)	Endoleptic Luvisol (A, E, BT, Mca, Rca)
	FDC	4	Calcaric Cambisol (Aci, Sca, Cca, Rca)	Endoleptic Luvisol (A, S/E, BT, Mca, Rca)
	FDC	5	Leptic Calcaric Cambisol (Aci, Sca, Cca, Rca)	Epileptic (truncated) Luvisol (A, BT, Mca,
				Rca)
	FDC	6	Endoleptic Calcaric Cambisol (A, Sci, Sca,	Epileptic Calcaric Cambisol (Aci, Sca, Mca,
			Cca, Mma)	Mma)
	FDC	7	Calcaric Cambisol (A, Sca, Mma)	Ferric Cambisol (A, S, S à sequioxydes;
				buried Calcaric Paleogley A, Bgca, Mca)
	FDC	8	Calcaric Cambisol (Aca, Bca, Mca)	Epileptic Calcaric Cambisol (Aci, Sca, Cca,
				Mca)
	FDC	9	Calcaric Cambisol (A, Sca, Rca)	Endoleptic Calcaric Cambisol (Aci, Sca, Cca,
				Mca)
	FDC	10	Endoleptic Calcaric Cambisol (A, Sca, Rca)	Endoleptic Calcaric Cambisol (Epiclayic) (A,
				Sca, BTca, Rca)
	FDC	11	Luvisol (Siltic) (A, E, BT, C) over buried	Epileptic (truncated) Luvisol (A, BT, Mca,
			Pretic Anthrosol (Leptic Calcaric) (Aca, Rca)	Rca)
	FDC	12	Endoleptic Endocalcaric Luvisol (A, E, BTca,	Endoleptic Luvisol (A, E, BT, Mca, Rca)
			Rca)
	FDC	13	Calcaric Cambisol (A, S, Sca, Mma)	Luvisol (A, E1, E2, BT, Mca, Mma)
	FDC	14	Luvisol (Siltic) (A, E, BT1, BT2, C, R)	Luvisol (Siltic) (A, E, BT1, C, BT2, R)
	FDC	16	Calcaric Cambisol (A, Sca, Rca)	Epileptic Calcaric Cambisol (Aca, Sca, Mca,
				Rca)
	FDC	17	Endoleptic Luvisol (A, E, BT, Rca), rounded	Luvisol (A, E, BT, Mca, Rca)
			and isolated limestone blocks inside BT
	FDC	18	Calcaric Cambisol (Aca, Sca, Cca, Rca)	Epileptic Calcaric Cambisol (Aca, Sca, Mca,
				Rca)
	FDC	19	Albic Protospodic Dystric Arenic Luvisol (A,	Albic Dystric Arenic Luvisol (A, E, BT, Mma,
			E, BP, BT, Msi)	Mca, Msi)
	FDC	20	Endoleptic Arenic Luvisol (Pretic) (A, E, BT,	Endoleptic Arenic Luvisol (A, E, BT, Mca)
			Mca)
	FDC	21	Arenic Dystric Luvisol (A, E1, E2, BT,	Arenic Dystric Luvisol (A, S, BT, Mca,
			Msi/Mca)	Msi/Mca)
	SGB	1	Gleyic Luvisol (Siltic Pretic) (A, Bi, Btg1,	Luvisol (Siltic) (A, E, BT, C, Cca)
			Btg2, C, Cca)
	SGB	4	Epileptic Calcaric Cambisol (Aca, Sca, Rca)	Epileptic Calcaric Cambisol (Aci, Sca, Rca)
	SGB	5	Chromic Luvisol (

A, E, BT, Cca/C, Rca) Endoleptique Luvisol (A, E,

BT, Mca, Rca)

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d

Table 4

 4 Comparative results of soil pH water analyses inside and outside the 25 CCDs.

	Soils	pH water	pH water	Wilcox.test
	Depths	inside	outside
	(cm)	CCDs	CCDs
	-5	5.70	5.36	0.2348
	-25	5.91	5.94	0.5416

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4259820 P r e p r i n t n o t p e e r r e v i e w e d

Acknowledgements

The authors would like to thank the French Ministère de la Culture which funded the ARPEGE project. We also thank the Office National des Forêts for its support of the project, in particular for making the Lidar available on the two forests. Thanks to Guillaume Decocq for his careful proofreading and Yves Le Bechennec for his help in the chronocultural determination of pottery fragments.


Yann Jobic
email: yann.jobic@univ-amu.fr

Marc Médale
email: marc.medale@univ-amu.fr

Frédéric Topin
email: frederic.topin@univ-amu.fr

Flow stability and regime transitions on periodic open foams

Keywords: Incompressible Navier-Stokes equations, Darcy's law, Forchheimer's law, Asymptotic Numerical Method, Linear Stability Analysis, Steady-state and Hopf bifurcations

This study aims to link critical Reynolds numbers associated with either steadystate or temporal bifurcations to fluid flow regimes described by macroscopic laws (Darcy, Forchheimer) in a periodic 3D Kelvin foam, using direct numerical simulations. We first identify the permeability and inertial coefficient of the Darcy's law and the Forchheimer's one. We explicit the different flow regimes that are accounted for in the Forchheimer's framework (Darcian, weak inertia and strong inertia regimes, respectively). We present an original systematic way to determine the critical Reynolds number associated with both regime transitions. We calculate them over a wide range of porosities and present two power-law correlations that locate these regime transitions for engineering purposes. We have performed pore-scale resolved calculations for various porosities with the Asymptotic Numerical Method (ANM) to find steady-state bifurcations, if any, along with a Linear Stability Analysis (LSA) to find temporal (Hopf) bifurcations from steady-state base flows. All the computed bifurcations occur at Reynolds numbers in the vicinity of the transition from weak to strong inertia regimes, where a change of behavior takes place. On the other hand, no bifurcation has been found in the transition between Darcian and weak inertia regimes.

Introduction

Predicting flow structures in porous media is a key task in several environmental and industrial applications (hydro-geology, geothermal energy, chemical engineering, etc., to quote a few). This rich field has been widely studied for the last two centuries. The pioneering approach is due to Darcy [START_REF] Darcy | Les fontaines publiques de la ville de dijon[END_REF], who identified a macroscopic law for creeping flows in a packed-sand beds. Specific macroscopic properties associated with the porous media were at that time only determined from experiments. It has been verified countless times and it linearly relates the macroscopic gradient of the fluid phase average pressure ∂ < p f > /∂x (Pa.m -1 ) to the filtration (or Darcian, superficial) velocity u (m.s -1 ), as follows:

- ∂ < p > f ∂x = µK -1 D u, (1)
with µ the dynamic viscosity of the fluid (Pa.s) and K D the permeability (m 2 ) of the porous media. Also from experimental observations, Forchheimer ( [START_REF] Forchheimer | Wasserbewegung durch boden[END_REF]) extended the Darcy's law to accommodate a wider range of velocities by adding a quadratic term into Eq. 1, which accounts for the fluid-flow kinetic energy. So, the Forchheimer's law reads:

- ∂ < p > f ∂x = au + bu m , ( 2
)
where a is a constant related to the first order polynomial in u, b is the inertia coefficient, that depends on the geometry of the porous media (e.g., strut shape, tortuosity, etc.), and the power coefficient m is close to 2. We have written eq. 2 in a general way to describe the different methods used for the identification of macroscopic parameters (a,b) that enter the law. They are related to the way one determines a and b by fitting experimental or from computations at pore-scale data-points. The uniqueness of the solution is not guaranteed and strongly depends on the working range of Reynolds numbers ( [START_REF] Kumar | State-of-the-Art of Pressure Drop in Open-Cell Porous Foams: Review of Experiments and Correlations[END_REF]).

There are three main ways to determine (a,b). The first one is to identify a from the permeability defined from the Darcy's equation, eq. 1. In a second step, one looks for b, keeping a fixed ( [START_REF] Kumar | Predicting pressure drop in open-cell foams by adopting forchheimer number[END_REF][START_REF] Kumar | Investigation of fluid flow properties in open cell foams: Darcy and weak inertia regimes[END_REF]). This method is mainly used with numerical simulations, as it's very difficult to obtain accurate experimental data in the Darcian regime ( [START_REF] Kumar | Investigation of fluid flow properties in open cell foams: Darcy and weak inertia regimes[END_REF])for such media. A second method determines b first as one intends to fit the high inertia regime (high Reynolds numbers). Then, a is determined with b fixed. Finally, the third, mostly used, method is to determine (a,b) altogether, by fitting directly the curve of the pressure drop as a function of the filtration velocity. The differences between those methods can be found in [START_REF] Kumar | Micro-structural impact of different strut shapes and porosity on hydraulic properties of kelvin-like metal foams[END_REF], which shows a difference of 20-97% for the permeability, and up to 6% for the inertia coefficient. In other words, this method while often used for engineering applications leads to coefficient values valid only for the range of identification and may leads to significant errors in the prediction of pressure drop outside of this range (e.g. for creeping flow conditions that are often extrapolated in industrial situations)

However, as the macroscopic flow law strongly depends on the flow regime, there is a need to correctly characterize the validity domain of the previous macroscopic laws, together with transition regions between them. Three zones are usually identified in the related literature [START_REF] Firdaouss | Nonlinear corrections to darcy's law at low reynolds numbers[END_REF][START_REF] Muskat | The flow of homogeneous fluids through porous media[END_REF]. The first one is called the Darcian zone, where viscous effect are dominant with respect to inertia terms, and corresponds at the pore scale to a Stokes flow. In this region, the Darcy's law is valid and accurate and predicts that the filtration velocity linearly depends on the macroscopic pressure drop. The third zone corresponds to the range of high Reynolds numbers, where turbulence occurs, characterized by unsteady and chaotic flows. In this zone, non-linear effects are dominant. The second zone is the transition zone in between the two previous ones. The inertia effects are of the same order of magnitude than the viscous ones. It is well established [START_REF] Mei | The effect of weak inertia on flow through a porous medium[END_REF][START_REF] Pauthenet | Inertial sensitivity of porous microstructures[END_REF] that we can divide this zone in two secondary sub-zones. The first one defines the limit from Darcy to what we call the weak inertia regime. In this sub-zone, the higher the filtration velocity, the greater the pressure drop deviates from the Darcy law (eq. 1). There is therefore a transition between the Darcy zone and the weak inertia one. Moreover, as the Reynolds number increases, the non-linear effects become dominant, while the flow remains laminar at the pore scale. We call this zone the strong inertia regime, which is different from turbulent flows. There is also a transition from weak inertia to strong inertia regimes or zones, to be defined. All those transitions may be related to the Reynolds number (Re = ρuL µ ), which obviously depends on the choice of the reference length-scale L. Note that the Forchheimer equation models at the macroscopic scale all of those zones with a reasonable accuracy.

Aside the two famous macroscopic laws (Darcy and Forchheimer equations), some other works have been devoted to recover macroscopic laws from theoretical developments. Different techniques were used, e.g., from dimensional analysis [START_REF] Ward | Turbulent flow in porous media[END_REF], capillary models [START_REF] Azzam | Flow rate-pressure gradient measurements in periodically nonuniform capillary tubes[END_REF], averaging methods [START_REF] Ruth | On the derivation of the forchheimer equation by means of the averaging theorem[END_REF], asymptotic expansions [START_REF] Giorgi | Derivation of the forchheimer law via matched asymptotic expansions[END_REF], volume averaging [START_REF] Whitaker | The forchheimer equation: A theoretical development[END_REF], etc. Those developments were correlated to the need to understand the physical meaning of the above macroscopic laws. However, the origin of the deviation from Darcy's law is not yet completely understood. Indeed, it was first attributed to turbulence ( [START_REF] Muskat | The flow of homogeneous fluids through porous media[END_REF], p.62), even for low Reynolds numbers (of values ranging from 1 to 10 in that reference). However, as the flow exhibits non-linear effects without turbulence at even lower Reynolds numbers, other explanations should be advanced. Giorgi ([14]) suggested that it could originate from the inertial microscopic forces and lead to a dependence of the inertia coefficient to the filtration velocity of the Forchheimer's law.

The increase of computational power enables nowadays to compute the Navier-Stokes equations at pore scale in realistic 3D geometries. Indeed, they quantitatively describe the fluid flow behavior for all Reynolds numbers, provided one assigns suitable spatial resolution to the considered hydrodynamic regime. Moreover, CFD computations at pore scale give access to a large amount of information, such as streamlines, fluid accelerations, flow variation around the mean one, drag forces, etc. So, by upscaling the computed flow fields, one can identify the macroscopic parameters.

Various numerical methods have been used: Finite Element Method ( [START_REF] Yazdchi | Towards unified drag laws for inertial flow through fibrous materials[END_REF]), Lattice Boltzmann Method [START_REF] Narváez | From creeping to inertial flow in porous media: a lattice boltzmann-finite element study[END_REF][START_REF] Jobic | Determining permeability tensors of porous media: A novel 'vector kinetic' numerical approach[END_REF], along with commercial softwares ( [START_REF] Kumar | Predicting pressure drop in open-cell foams by adopting forchheimer number[END_REF][START_REF] Jobic | Transport properties of solid foams having circular strut cross section using pore scale numerical simulations[END_REF]), etc. The main idea is to find laws describing the dependency of the macroscopic properties of the porous media to intrinsic geometric parameters (such as porosity, tortuosity, etc.). We can cite the famous Ergun's law ( [START_REF] Ergun | Fluid flow through packed columns[END_REF], which links the pressure drop to the porosity and the grain size. This work leads to many variations (experimental and numerical ones) depending of the corresponding field (a review can be found in [START_REF] Yazdchi | Towards unified drag laws for inertial flow through fibrous materials[END_REF], table 2). Moreover this formulation developed for packed bed of spheres has been applied to many other media and it has been shown that the macroscopic flow laws depends also on other geometric quantities.

The bridge between computations at the pore scale and macroscopic laws is the existence of elementary representative volume from some periodic geometry that could be found in the porous media. Among idealized elementary representative volumes is the Kelvin cell, for which a comprehensive survey of macroscopic flow laws is presented in [START_REF] Kumar | State-of-the-Art of Pressure Drop in Open-Cell Porous Foams: Review of Experiments and Correlations[END_REF]. Among all the pressure drop studies and correlations presented in this review paper, critical remarks are summarized as the follows:

The permeability and inertia coefficients are very sensitive to the porosity range (mainly low-porosity range) and are strictly strut shape dependent; Separate determination of flow law characteristics (K D and C For in viscous and inertia regimes) is needed. Extraction of permeability from second-order polynomial pressure drop data would lead to wrong values of Darcian permeability in the viscous regime; Values of Ergun's parameters vary clearly with porosity and strut shape and do not possess constant numerical values; Common definition of characteristic length scale (and associated measurements) is missing to reduce the dispersion in friction factor data.

Thermal instabilities in foams are first studied at the macroscopic level, using Darcy's laws for the flow (see for example [START_REF] Ryland | A bifurcation study of convective heat transfer in porous media. Part II: Effect of tilt on stationary and nonstationary solutions*[END_REF][START_REF] Riley | Time-periodic convection in porous media: the evolution of hopf bifurcations with aspect ratio[END_REF][START_REF] Weinitschke | A bifurcation study of convective heat transfer in porous media[END_REF]). Those studies aim to examine the convective heat transfer in a homogeneous porous duct of rectangular cross section. They link the stability to the Rayleigh number, and find for those configurations some Hopf bifurcations. Another classical approach, also at macroscopic level, is to study the stability of the flow overlaying a porous layer, also modeled at the macroscopic level (see for e.g. [START_REF] Samanta | Linear stability of a plane couette-poiseuille flow overlying a porous layer[END_REF][START_REF] Liu | Instabilities of a liquid film flowing down an inclined porous plane[END_REF][START_REF] Sengupta | Stability of poiseuille flow of a bingham fluid overlying an anisotropic and inhomogeneous porous layer[END_REF]). They are interested in the interaction between the flow and the porous layer, in term of production or dissipation of kinetic energy, linking it to instabilities. However, the study of fluid-flow instabilities at pore level in 3D foam is a new topic. To our knowledge, only one paper reported 2D bifurcations into periodic ordered and disordered non-deformable porous structures, composed of squares or circles ( [START_REF] Agnaou | From steady to unsteady laminar flow in model porous structures: an investigation of the first hopf bifurcation[END_REF]). For each porous medium, a systematic study of increasing mass-flow rate was performed for unsteady Navier-Stokes equations, giving the critical Reynolds numbers of the Hopf bifurcations. Instabilities of 2D fluidic pinball has also been studies [START_REF] Deng | Coinciding local bifurcations in the navier-stokes equations[END_REF]. We can note that 3D flow and/or instabilities have no reason to be identical to 2D ones.

For example, critical Reynolds number of the Hopf bifurcation for the flow over a 2D circular cylinder arise at about Re c =43, for the 3D cylinder Re c =188, and for the 3D sphere, Re c =277. Therefore, critical Reynolds number of the same shape family leads to very different values in 2D and 3D. In the present work we are interested in intrinsic 3D geometries, presenting their own complexities.

So, to address some of the above remarks, the present work aims to compute the incompressible Navier-Stokes equations at pore level in a Kelvin's like foam to provide reference values for the identification of macroscopic laws. Moreover, we also have computed the instabilities that take place and identified critical Reynolds numbers at which they occur. Indeed, we have computed both steady-state and Hopf bifurcations, depending on the fluid-flow regime. The outcome of this work is to relate the computed instabilities to the transitions in macroscopic flow regime predicted by macroscopic laws or correlations.

Numerical methods

We are interested in the study of real foams, however, they are too complex as we can not control individually each parameters (e.g., strut shape, pore elongation, different local porosities, etc). Moreover, the individual manufacturing process of each producer (even each batch) lead to specific relations between the geometrical parameters of the commercially available samples. Therefore, the Kelvin foam structure, a widely recognized idealized structure of real foams, has been chosen. Its symmetries enable to restrict the computational domain to the periodic pattern and thus strongly reduce computational cost. The chosen strut shape is triangular, with various cross-section to change the porosity. The present study is deliberately focused on the high porosity range (68% < ε < 97%), which is of industrial application relevance.

At pore scale, we consider an incompressible fluid-flow in laminar regime, which is governed by the Incompressible Navier-Stokes equations:

∂v ∂t + (v • ∇)v = -∇p + 1 Re ∇ 2 v, ∇ • v = 0, (3)
with v the velocity vector at the x position, p is the pressure, Re is the Reynolds number based on the relevant length-scale L (Re = ρU L/µ) and t is time.

Periodical boundary conditions are imposed along the streamwise direction (Z direction), meanwhile symmetry conditions are applied along the span-wise horizontal directions (X and Y directions). The periodic condition is implemented by changing the connectivity of the elements While the symmetry conditions are obtained by cancelling the transverse gradient on the appropriate surfaces As we consider the streamwise periodicity of the flow, the fluid flow is driven by a macroscopic body force along the Z direction.

The governing equations are discretized in space by the Galerkin Finite Element Method [START_REF] Dhatt | Methode des éléments finis -une présentation[END_REF][START_REF] Carey | Penalty finite element method for the navier-stokes equations[END_REF], with tri-quadratic Lagrange polynomial approximations for the velocity and a linear one for the pressure to satisfy the Ladyzhenskaya-Babuška-Brezzi (LBB) condition. Three complementary algorithms have been used in the present work: a steady-state continuation algorithm, based on the Asymptotic Numerical Method (ANM, [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Medale | High performance computations of steady-state bifurcations in 3d incompressible fluid flows by asymptotic numerical method[END_REF]) to compute the branches of steady-state solutions; a Linear Stability Analysis algorithm to compute the generalized eigenvalue problem associated with the linear stability of steady-state solutions ( [START_REF] Nicolis | Introduction to Nonlinear Science[END_REF][START_REF] Medale | A one-dimensional full-range twophase model to efficiently compute bifurcation diagrams in sub-cooled boiling flows in vertical heated tube[END_REF]); a time integration algorithm of the Incompressible Navier-Stokes equations, based on an unconditionally stable projection method ( [START_REF] Guermond | A new class of truly consistent splitting schemes for incompressible flows[END_REF]) to compute unsteady fluid flows.

When looking for a steady solution, when it exists, we have used a penalized formulation of the governing equations which enables direct access to the solution. However, at high Reynolds numbers this process leads to solve a highly non-linear problem. Therefore, one uses the ANM to compute this non-linear solution with a continuation algorithm( [START_REF] Cochelin | Méthode asymptotique numérique[END_REF]). Moreover this continuation algorithm computes the various branches of steady-state solutions if present, i.e. when it exists steady-state bifurcations. The ANM is based on the combination of a perturbation method and a continuation one. The perturbation method introduces high-order Taylor series representation of the unknowns along with the continuation parameter. So, it transforms a highly non-linear and implicit parametric problem into a series of linear ones. As a continuation algorithm, it is based on the implicit function theorem that enables to extrapolate a regular (continuous and derivable) initial solution to its vicinity along the loading parameter for a regular solution. It is therefore a high-order predictor/corrector algorithm in which the correction stage is performed with a Newton-Raphson algorithm, if needed.

To overcome the very poor conditioning of the penalized formulation, the resulting algebraic systems are solved with the direct parallel MUMPS (LU parallel) solver ( [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF][START_REF] Amestoy | Performance and Scalability of the Block Low-Rank Multifrontal Factorization on Multicore Architectures[END_REF]). This formulation can significantly reduce the computational cost of steadystate problems, compared to any time-marching algorithm. On the other hand, this method is costly in terms of RAM, especially for 3D computations.

The implementation of the steady-state solver with the ANM continuation algorithm is performed into the massively parallel PETSc framework library ( [START_REF] Balay | PETSc users manual[END_REF][START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF][START_REF] Balay | PETSc Web page[END_REF]).

The ANM algorithm enables to compute all steady-state solutions (stable and unstable ones), but also to detect and compute steady-state bifurcations, if any. Therefore it has to be complemented by a Linear Stability Analysis (LSA) in order to find out other bifuraction types, such as Hopf bifurcations.

In the present work, the LSA of steady-state solutions relies on the normal mode analysis. This method consists in the linearization of the Navier-Stokes equations by applying a perturbation of exponential form in time. The solution can this way be decomposed into normal modes. The stability analysis of the steady-state solution is achieved by solving a generalized eigenvalue problem owing to the Finite Element discretization of the dynamical terms (the mass matrix). We apply this method to each steady-state solution computed by the ANM. The implementation of this linear analysis aglorithm is performed with the SLEPc library ( [START_REF] Hernandez | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF]), which relies on the PETSc one.

For the unsteady calculations, the use of a projection algorithm leads to a "segregated" approach, which consists in successively compute velocity-pressure subproblems. This choice has two main consequences from a numerical and computational point of view. The first consequence is that the size of each sub-system is reduced, and more importantly that the conditioning of each algebraic sub-system is significantly improved compared to the original global problem. This allows us to efficiently solve these systems with iterative methods (of preconditioned Krylov subspace type), and this on massively parallel architectures. Regarding the solution of the incompressible Navier-Stokes equations, we have implemented the inconditionnaly stable algorithm of Guermond-Shen ( [START_REF] Guermond | A new class of truly consistent splitting schemes for incompressible flows[END_REF]). The time integration is performed with a first order Backward Euler scheme with appropriate time steps.

Results and discussion

The quality of the numerical solution depends strongly on the quality of the mesh, specially for the presented numerical tools. Indeed, the ANM relies on high order power series in a recursive way. The numerical errors could then adds up in a dramatic way. We then want hexahedral meshes as it's the one that offers accurate solutions compared to others (tetrahedral, pyramids, ...), provided they are not too much distorted (invertible geometric transformation from parametric to actual elements). However, the mesh should correctly describe the geometry (negligible geometric discretization errors). Note this structure is an idealized Kelvin Cell (see Fig. 1) as it is composed of straight struts of nearly equilateral triangle cross-section (often named Kelvin Like Cell) while the original Kelvin Cell is composed of sligthly curved struts as it is the dual of stacked interpenetrating spheres. Gyroid surfaces are also used to create Kelvin cell structure. There are many ways to create these geometric structures. The basic idea is to start with a truncated octahedron. Then a first option is to create strut cross section centered along the edge of the face and extrude them; this method is simple, fast as there is many symmetry planes, allow to parametrically create various strut shape and cross section and associated unstructured meshes [START_REF] Kumar | The geometric and thermohydraulic characterization of ceramic foams: An analytical approach[END_REF]. On the other hand, the subdivision of the poral space (usually obtained via boolean subtraction of the assembled struts from a bounding cube) in hexahedral patches if needed is tedious and complicated. We thus choose to construct directly the poral space. The domain was obtained by extruding (outward, normal) several faces of the truncated octahedron and taking advantage of the symmetries. The resulting structure is the poral space of a triangular strut Kelvin cell formed by the empty spaces left by the extrusion process [START_REF] Klostermann | Meshing of porous foam structures on the micro-scale[END_REF]. This space is then easily subdivided in hexahedral shapes that are paved with structured hex meshes. The resulting mesh is isotropic and periodic with only high-quality hexahedral elements, as presented Fig. 2. This conformal mesh possesses a complicated periodicity, consisting in one cell plus two half-cells in two directions, and one in the last one. It's however periodic in all the directions. We can generate different porosities by changing the extrude length (which correspond to the strut cross section edge length), which also change the size of the Kelvin cell. The inner volume of the cell (the original truncated octahedron) stay the same, but the distance between struts centerline change and thus also the size of the bounding cuboid enclosing the Kelvin cell.

This domain discretization quality is of primary importance for numerical precision at higher Reynolds numbers, as the inertia term comes into play. A previous attempt has been made to use a regular structured grid generated by a marching cube algorithm applied to the Kelvin cell. In this case, the unstructured mesh is perfect in the numerical way (small cubes), but it is aligned with the Cartesian coordinate system and does not follow correctly the geometry, creating a staircase effect on the struts-fluid interfaces. Refining the mesh creates different porous media (the geometrical discretisation errors change with the mesh size), with different properties (mostly porosity, specific surface area, and the interface roughness that is artificially added). Considering these bias we prefer a initially more complicated but far more efficient mesh creation procedure to ensure that mesh convergence studies were carried out on the exact same geometry. Note : The dependency of the permeability ( D 2 bc ), inertia coefficient ( D -1 bc ) as well as the geometrical parameters of the Kelvin cell with the size of the Kelvin cell (here the dimension of the bounding cube -D bc -) is well known and could be removed from all further analysis [START_REF] Kanaun | Effective conductive properties of open-cell foams[END_REF]. Figure 3 depict the results of a pore scale simulation at Re Dh = 183.7, using the hydraulic diameter Dh as the length scale (see below). We can observe the complicated flow path going around the struts,with the streamlines highlighting some recirculation patterns. That dynamic is complex, and is difficult to extract simple -scalar-indicators to correlate to the shape of the struts.

Coefficient identification for macroscopic laws

Identifying flow properties in porous media is critical for distinguishing regime transitions, if any. As highlighted in the introduction, there are several methods to identify them. We have chosen the one that makes the best use of the our numerical tools. We have chosen the following form of the Forchheimer's equation:

- ∆p L = µK -1 D u + βρu 2 , (4)
with β the inertial coefficient (m -1 ) and ρ is the fluid density (Kg.m -3 ). The identification method we have used is as follows. First, we use the computed steady-state solutions in the range of very small Reynolds numbers (Re ∈ [10 -5 , 1]) to identify the Darcian permeability K D from eq. 4 in which we drop the inertial term (see Figure 4a for the case of porosity ε = 0.68). There are at least seven computations in this range that enable us to perform a highly accurate identification of the K D coefficient (regression coefficient: 0.9999534). By doing so, we take advantage of the FE method that is particularly efficient for elliptic problems. In a second step, we take the deviation from the Darcy's law as follows:

(-∆p/L)/u -K D u. This quantity is a linear function of u 2 , according to eq. 4. We plot this quantity and perform the linear fit with results from higher values in velocity. From this fit, we can group the pair (u, ∂ < p > /∂x) until the deviation to the linear fit does influence the regression coefficient. We have shown the case of porosity ε = 0.68 in Figure 4b.

We have at least 50 continuation points in this range and get regression coefficient of 0.9999. Finally, we extract the inertia coefficient β, which is the slope of this linear fit.

Regarding the case of open-cell foams as the one considered in this work, one has to select one among the various characteristic length scales L c (m), to build up dimensionless quantities, such as:

The Reynlods number:

Re = ρuL c µ . (5)
The Darcy friction coefficient:

f D = ∆P L c 2L c ρu 2 . ( 6
)
There is no consensus in the literature for the relevant choice of L c . The hydraulic diameter, a 3D extension to the cylinder diameter historically used in fluid studies, defined as D h = 4V f /S p , with V f the volume of the fluid phase (m 3 , called the total wetted volume), and S p the total wetted surface area (m 2 ) are commonly used.

Other reference quantities are also found in the literature, as the strut length L s , the periodic length of one cell foam L p , √ K D only for Darcian flows, or βK D for Forchheimer's ones, etc (e.g. [START_REF] Zami-Pierre | Transition in the flow of power-law fluids through isotropic porous media[END_REF][START_REF] Mills | The wet kelvin model for air flow through open-cell polyurethane foams[END_REF]), reported in table 1. The first three (D h , L p , L s ) are purely geometrical defined quantities. D h and L p are macroscopic quantities whereas L s is a 10 times smaller quantity than the two previous ones. Finally the fourth one, √ K D combines geometrical quantity (porosity, ...) and macroscopic integration of the pressure drop over the cell.

Table 1: Permeabilities, inertia coefficients with length scales for a Kelvin cell of various porosities ε.

Macroscopic properties

Length scales ). Previous results were obtained with different numerical tools such as Starccm+ (commercial software) and a in-house code using a kinetic method, close to Lattice Boltzmann ( [START_REF] Jobic | Determining permeability tensors of porous media: A novel 'vector kinetic' numerical approach[END_REF]). The values computed in the present work (Kelvin cell with triangular struts, numerical model described in §2) are reported in Table 1 and also reported in Figure 5. There is a very good agreement between all methods, whatever being the strut shape.

Porosity K D β D h Lp Ls √ K D (%) (m 2 ) (m -1 ) (m) (m) (m) (m

Flow regime transitions

Using the dimensionless pressure drop in the form of friction coefficient as a function of the Reynolds number leads to the Moody's diagram. It's the usual way to qualitatively identify transition regimes. Using Fig. 8, we recover three regimes: i) the Darcian one in the linear part with the minus one slope; ii) a first transition departing from the previous linear zone and corresponding to the weak inertia zone; iii) the strong inertia zone. The curves are ordered from top to bottom for increasing porosity (68.3% < ε < 97.6%). Indeed, at a given Reynolds number, the friction coefficient increases as the porosity decreases (as it is proportional to the specific surface area).

When increasing the Reynolds number above the Darcian region (of minus one slope in Fig. 8), the departure from the Darcy's law grows as Re increases. So, one computes the distance to the Darcy's law and defines a user defined threshold to identify the critical Reynolds number at which the flow gets out from the Darcian regime. The pressure-drop relationship to the filtration velocity is nearly linear in this zone, and changing this threshold only shifts the critical Reynolds number. For example, taking 1% over 2.5% changes the critical Reynolds of about 25%. In this work, we have chosen this threshold to be at 2.5%. Fitting the values reported in Table 2, one finds that the Reynolds number evolves as a function of the compacity (1-ε) according to the power law: Re Dh = 13.35(1 -ε) 0.456 . [START_REF] Firdaouss | Nonlinear corrections to darcy's law at low reynolds numbers[END_REF] Moving to the transition between weak and strong inertia regions, another method has been used. We represent the flow as the deviation to Darcy's law as a function of u 2 (see Fig. 9). This curve has the weak and strong zones. For simplicity, we can model their respective evolution by a polynomial, for each zone. We found that there is a unique intersection point of those polynomials, representing the point where the transition occurs. With this method, we do not need to define a threshold. The critical Reynolds number related to that transition follow another power law, as follows:

Re Dh = 2154.4(1 -ε) -0.823 . (8)
We have compared the values of transitional critical Reynolds number based on different length scales, in Table 2 and3, summarised in Figure 6. Obviously, the slope's curves depend of the variation of the considered reference length scale. The zone defining the weak inertia flow zone is always of two orders of magnitude in Reynolds number, but is shifted from 0.01 for the strut size considered as the length scale, to 10 for the hydraulic diameter. The shift is so huge that the critical Reynolds transition from weak to strong inertia based on the strut size is at the same order of magnitude of the end of the Darcian regime based on the periodic length.

We can make a mechanical analysis to chose a particular length scale as the characteristic one. Two factors play a high role in modifying the pressure drop: the change of cross section, and the change of direction due to the path taken by the fluid in the porous media. The first one is related to the divergence free velocity field, changing the flow velocity between each section that the fluid cross. The second one is related to the tortuosity of the porous media. Thus, the friction length is related to both parameter. From all the possible length scales cited above, the hydraulic diameter could be seen as a friction length, and therefore could be more accurate to express the changes at which the macroscopic regimes occur. We use from now on the hydraulic diameter as the characteristic length scale. The first one is the Hopf bifurcation (or dynamic one), and the second one is the steady-state bifurcation (or static one).

We have performed a spatial convergence study (see Table 4) to check the mesh sensitivity on the critical Reynolds number associated with the first encountered bifurcation. We have used a smaller (quarter) representative domain, which enables to get a finer mesh. One effectively gets a convergence to the same critical Reynolds as the mesh gets finer. Applying the LSA to all steady-state solutions computed by the ANM, we find the critical Reynolds numbers associated with the first Hopf bifurcations, reported in Table 5. A steady-state bifurcation study has been done for various porosities and is summarized in Table 6 and Figure 7. The higher the porosity the higher is the Reynolds number for a given pressure drop. For such high porosities the first encountered bifurcation is a Hopf one. This means that the steady-state flow looses at that point its stability towards a periodical limit cycle, telling us that inertia changes the dynamics of the system. On the other hand, for low porosity cases, the Reynolds number is much lower for the same given pressure drop and the first encountered bifurcation is a steady-state one, whose nature is reported in table 6.

The Reynolds numbers associated with such steady-state bifurcations evolve as a function of the porosity according to the following power law:

Re Dh = 4778.4(1 -ε) -1.03 (9)
Two main points are noteworthy. The first one is that in the range of spanned flow rates we found two Hopf bifurcations for the Kelvin foam at porosity 86.4%. It's not uncommon, but we only have this behavior for this porosity. The second remark is that in the range of spanned flow rates we didn't find any Hopf bifurcation for the porosity of 80.1%. In this case, the main branch along the control parameter's path computed by the ANM is unstable. We were able to take the first branch at the steady-state bifurcation point, which is a pitchfork, and take the stable branch. But this bifurcation act like a strong attractor, making it difficult to explore this branch, as this would need further algorithmic developments. Moreover, from a more practical point of view, the ANM algorithm may sometimes detect some strange bifurcations, specially at low Reynolds number. To check the physical or spurious origin of the computed bifurcations, the kernel size of the algebraic system was also computed. Indeed, at a bifurcation point, the matrix associated to this algebraic problem is rank-deficient. That means that there is different orders in magnitude between the smallest eigenvalues, some at zero to the numeric precision. For each bifurcation candidate, we have checked this property, and we are thus aware of the physical relevance of the steady-state bifurcation. The spurious bifurcations arise when the truncation order of the power series computed by the ANM is too high compared to the spatial discretization. In this case, the errors coming from the spatial discretization grow as the terms of the series expansion get smaller. The bifurcation detection is perturbed by this artefact. 3.4 Discussion: pore-scale instabilities and impact on flow law at macroscopic-scale

After identifying the transition between the various flow regimes (creeping, weak and strong intertia), and found the critical Reynolds numbers associated with the encountered bifurcations (Hopf or steady-state), we can compare their respective position in a same plot (Fig. 10). The first comment is that all bifurcations are located in the vicinity of the transition between the weak and strong inertia regimes. They occur at moderate Reynolds numbers that depend on the foam porosity, (expressed in term . of compacity). The Hopf bifurcations are mainly located just before this transition, whereas the steady-state ones are located just beyond it. Moreover, as the porosity decreases, the steady-state bifurcations occur much closer to the transition between the weak and strong inertia regimes. Both steady-state and Hopf bifurcations have at least one eigenvalue that has a null real part. The hopf bifurcation has a non null imaginary part, whereas the steadystate one is null. They both indicate a stability change, where any non-infinitesimal perturbation results in a change in the fluid-flow pattern or feature. It thus makes sense to gather them and represent their evolution using the same fitting curve, which has been shown in Fig. 10 (solid black line). This curve is almost confounded with the transition from weak to strong inertia regimes, which has been clearly established. The bifurcations which are fluid-flow instabilities at pore-scale, are then located in the area where the inertia terms are more significant as compared to the viscous forces. They happen at a relatively moderate Reynolds number and clearly depend on the porosity. No more detailed influence of foam geometry (e.g. strut shape) has yet been obtained.

Conclusion

The incompressible Navier-Stokes equations have been computed at the pore-scale in a Kelvin foam in order to get all knowledge necessary to identify a macroscopic behavior in such periodic geometry. We have identified macroscopic properties such as the permeability and the inertia term of a complex idealized 3D foams. To do so, we considered the first order polynomial of the Forchheimer's equation to be the Darcian permeability and identified in a subsequent step the inertia term.

We have introduced a systematic way to define flow transitions from three different regimes (creeping, weak and strong inertia), in which the flow has different macroscopic behaviors. They all follow a power law equation, whose fitting parameters have been given.

To our knowledge, it's the first time that the ANM and LSA method were used to automatically detect bifurcations in such complex geometrical fluid-flow configurations. We found some steady-state and Hopf bifurcations, for Kelvin cells with various porosities. We have correlated the critical Reynolds numbers related to those bifurcations to the flow regimes and showed that they appear in the vicinity of the transition from weak to strong inertia regimes. This strong result validates our methodology in defining flow transitions, where flow behavior changes. We may also emphasize that the used numerical strategies are computationally efficient, as they rely on steady-state solutions to perform linear stability analyses.

We partially used the field information given by the direct numerical solutions in such spatially resolved configurations. Some other valuable quantities can be exploited in a future work, such that the eigenvector field associated with the most unstable mode given by the LSA, which provides the most dangerous perturbation field along with its magnitude. It will help to better understand the flow behavior, linking it to the foam geometry.

  (a) A Kelvin cell. (b) Exploded view of the same cell.

Fig. 1 :

 1 Fig. 1: Solid phase of a Kelvin cell with triangular struts (porosity ε = 68.3%).

Fig. 2 :

 2 Fig. 2: Left-half computational fluid domain, right-half solid-fluid interface. Porosity: ε = 92.5%, mesh composed of 290'000 hexahedral elements.

Fig. 3 :

 3 Fig. 3: Computed results at Re Dh = 183.7 for porosity ε = 86.4%. Ribbons (streamlines of the flow) and velocity vectors (both colorized by velocity magnitude).

3. 3

 3 Pore-scale analysis: fluid flow bifurcations Using the numerical tools presented in Section 2, two types of bifurcations have been encountered while spanning the fluid-flow rate in the range of Re Dh ∈ [0 -2500].

  Determination of the permeability K D . Determination of the inertia coefficient β.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Identification of the macroscopic flow properties of a Kelvin cell for the porosity ε = 0.68.

Fig. 6 :Fig. 7 :

 67 Fig.6: Critical Reynolds number for the transition from Darcy to weak inertia regimes (dotted lines) and for the weak to strong inertia regimes (solid lines). The square marks refer to the length scale based on the hydraulic diameter, the diamond marks for the periodic length, the triangular marks for the strut size and the circle marks for √ K D .

Fig. 8 : 4 Fig. 9 :

 849 Fig. 8: Moody's diagram.

Fig. 10 :

 10 Fig. 10: Validity domain of macroscopic laws along with the critical Reynolds numbers associated with the computed bifurcations. Transition between regimes are represented by dashed lines. Steady-state bifurcations are represented by a cross, Hopf bifurcations by a squares. The fit of all bifurcations is represented by a black line.

  However, it's possible to compare the permeability obtained from previous work with comparable Kelvin cells, since it varies proportionally to the square of the cell size (i.e., K D ∝ d cell 2

						)
	97.6	0.1190	0.1441	5.17 2.54	0.1	0.3444
	92.5	0.0838	0.5989	3.18 2.77	0.2	0.2895
	86.4	0.0646	0.9974	2.52 3.00	0.3	0.2542
	80.1	0.0519	1.3677	2.19 3.23	0.4	0.2278
	68.3	0.0353	2.2202	1.86 3.69	0.6	0.1878

Table 2 :

 2 Reynolds number of the right bound of the Darcian validity domain, based on several length scales.

			Rec (Darcy) based on
	porosity (%) Compacity (%)	D h	Lp	Ls	√	K D
	97.6	2.4	8.59 4.21 0.17	0.57
	92.5	7.5	5.61 4.88 0.35	0.51
	86.4	13.6	4.22 5.02 0.50	0.43
	80.1	19.9	3.45 5.10 0.63	0.36
	68.3	31.7	2.61 5.19 0.84	0.26

Table 3 :

 3 Critical Reynolds number of the transition from weak to strong inertia regimes, based on several length scales.

			Rec (strong inertia) based on
	porosity (%) Compacity (%)	D h	Lp	Ls	√	K D
	97.6	2.4	1099.2 539.3 21.2	73.2
	92.5	7.5	467.0	406.2 29.3	42.4
	86.4	13.6	228.8	272.2 27.2	23.1
	80.1	19.9	184.4	272.1 33.7	19.2
	68.3	31.7	135.4	269.1 43.7	13.7

Table 4 :

 4 Spatial convergence study for the critical Reynolds number Re c (based on the hydraulic diameter) associated with the first Hopf bifurcation. h is the mean edge size of all the elements of the domain.

	porosity = 97.6% mesh #1 mesh #2	mesh #3	mesh #4
	mean h	0.0601	0.0451	0.0361	0.0285
	Number of dof	262 466	611 750	1 182 650 2 382 666
	Rec Hopf	1376	1390	1399	1410
	Relative error	2.41%	1.38%	0.76%	-

Table 5

 5

	:	Critical
	Reynolds number asso-
	ciated with the first
	Hopf bifurcation.
	porosity (%) Rec Hopf
	97.6	710.7
	92.5	294.1
	86.4	317.5
	80.1	-
	68.3	154.0

Table 6 :

 6 Critical Reynolds associated with the first steadystate bifurcation, with the bifurcation type

	porosity (%) Rec for steady-state bifurcation bifurcation type
	97.6	2077.5	transcritical
	92.5	485.6	pitchfork
	86.4	385.7	pitchfork
	80.1	209.7	pitchfork
	68.3	137.6	transcritical

Acknowledgments. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources.