# From Outcrop to Groundwater Flow: The Impact of Overlooked Heterogeneity on Advective Transport in Lower Triassic Sandstones, Eastern France

Lucas Bofill, Gerhard Schäfer, Guilherme Bozetti, Jean-François Ghienne, Mathieu Schuster

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

## ► To cite this version:

Lucas Bofill, Gerhard Schäfer, Guilherme Bozetti, Jean-François Ghienne, Mathieu Schuster. From Outcrop to Groundwater Flow: The Impact of Overlooked Heterogeneity on Advective Transport in Lower Triassic Sandstones, Eastern France. Geothermics, 2025, 136, pp.103581. 10.1016/j.geothermics.2025.103581. hal-05478872

HAL Id: hal-05478872

https://univoak.hal.science/hal-05478872v1

Submitted on 27 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

# From Outcrop to Groundwater Flow: The Impact of Overlooked Heterogeneity on Advective Transport in Lower Triassic Sandstones, Eastern France

Lucas Bofill<sup>a,b,*</sup> <sup>□</sup>, Gerhard Schäfer<sup>a</sup> <sup>□</sup>, Guilherme Bozetti<sup>a</sup> <sup>□</sup>, Jean-François Ghienne<sup>a</sup> <sup>□</sup>, Mathieu Schuster<sup>a</sup> <sup>□</sup>

## ARTICLE INFO

### Keywords:

- Sedimentary heterogeneity

- Advective transport

- Buntsandstein

- Upscaled anisotropy

- Breakthrough curve

- Groundwater flow

## ABSTRACT

Sandstone aquifers lacking mudstone interbeds are often considered homogeneous. However, depositional processes and diagenesis can generate significant permeability contrasts without notable variations in granulometry. Capturing this heterogeneity in flow models is essential but challenging. Outcrop studies provide valuable analogues for subsurface conditions, revealing sedimentary architectures and facies distributions not observable from borehole data alone. This study focuses on the Lower Grès Vosgien formation in eastern France, an important aquifer hosting lithium-rich geothermal brines in the Upper Rhine Graben, deposited by braided fluvial and aeolian systems during the Lower Triassic. This work aims to evaluate the influence of decimetre- to metre-scale sandstone features on advective transport modelling, through a two-dimensional conceptual model developed using a digital outcrop model and sedimentary facies analysis. Six progressively simplified hydrostratigraphic models are tested to assess how heterogeneity and its representation affect particle residence times, breakthrough curves, longitudinal macrodispersivity, and upscaled anisotropy. Results indicate that distinguishing fluvial- and aeolian-related sandstones significantly influences both horizontal and vertical advective flow. Failure to discretise individual fluvial facies and reliance on deterministic hydraulic conductivity values led to the underestimation of preferential flow pathways, delaying the prediction of first particle arrivals. Moreover, the discretisation between fluvial and aeolian-related deposits significantly affects macrodispersion results and yields considerable anisotropy when the conceptual model is upscaled. The strong coupling between advective flow with heat transfer and solute transport underscores the critical role of the observed sedimentary heterogeneity on the accurate understanding of lithium-rich geothermal brines circulation in the matrix porosity of the Lower Triassic sedimentary successions in the Upper Rhine Graben reservoirs.

## 1. Introduction

Accurately characterising subsurface decimetre-scale heterogeneities in sedimentary aquifers is practically impossible due to the complex spatial variability of sedimentary deposits (Miall, 1988; Høimyr et al., 1993; Zheng and Gorelick, 2003). Nonetheless, these characteristics have direct implications for subsurface fluid flow, heat-transfer, and solute transport (Gelhar et al., 1992; Zheng and Gorelick, 2003; Ferguson, 2007; Ramanathan et al., 2010; Possemiers et al., 2015; Ren et al., 2022). Laboratory columns or tank experiments and field tracer tests are commonly used to investigate the effects of such heterogeneities on

solute transport (Herr et al., 1989; Schäfer, 1991; Gelhar et al., 1992; Ma et al., 2025). Nevertheless, despite the valuable insights obtained by these approaches, both present limitations. The former is limited in space and, therefore, not able to replicate the scale effect of field heterogeneity (Gelhar et al., 1979). The latter is time-consuming, presents infrastructure costs, and does not allow direct observations to link solute behaviour and subsurface heterogeneity.

Geostatistical approaches based on variograms and covariances are commonly used to predict unknown sedimentary heterogeneity architecture. However, field studies have shown that these methods do not fully capture the natural complexity of sedimentary systems and often

<footnote><sup>a</sup> Institut Terre et Environnement de Strasbourg (UMR7063), CNRS, Université de Strasbourg, Strasbourg, France</footnote>

<footnote><sup>b</sup> Azerbaijan State Oil and Industry University, French Azerbaijani University, Baku, Azerbaijan</footnote>

<footnote>* Corresponding author at: 33 Uzeyir Hajibeyov, Baku, Azerbaijan.</footnote>

<footnote>E-mail address: medeirosbofill@unistra.fr (L. Bofill).</footnote>
fail to reproduce the connectivity of preferential flow paths (Webb and Anderson, 1996; Fogg et al., 1998; Klise et al., 2009). Forward stratigraphic modelling approaches help bridge this gap by linking sedimentary processes to their resulting deposits and improving the prediction of petrophysical heterogeneity (e.g. Snider et al., 2021; Yan et al., 2024). Nevertheless, there are currently no models available that can accurately replicate the heterogeneities at the scale and with the specific characteristics of the depositional system studied here.

Outcrop studies offer valuable analogues for subsurface conditions (e.g., Jussel, 1989; Alexander, 1993; Aigner et al., 1996; Rogiers et al., 2013; Howell et al., 2014; Weissmann et al., 2015; Thompson et al., 2019), providing insights into the geometry and distribution of sedimentary facies that are difficult to obtain from borehole data alone. Although hydraulic properties measured from outcrops may not exhibit similar absolute values compared to subsurface conditions (e.g. Koltermann and Gorelick, 1996; Rogiers et al., 2013; Busch et al., 2022; Bofill et al., 2025), outcrop measurements can be useful in a relative sense (Koltermann and Gorelick, 1996; Bofill et al., 2025).

Sedimentary facies analysis using outcrops has been widely used to enhance the efficiency of hydrocarbon exploration (Walker, 1992; Høimyr et al., 1993); however, its application remains underutilised in groundwater studies (Crooijmans et al., 2016; Åberg et al., 2021), despite the need of more comprehensive data collection and analysis of sedimentary features, for instance, to better understand their influence on solute transport behaviour (Anderson et al., 1999; Ma et al., 2025) and heat-transfer in geothermal reservoirs (Aghaei et al., 2024). Integrating quantitative facies analysis with field measurements of permeability from outcrops enhances the realism of hydrogeological models and improves our understanding of subsurface flow dynamics (Alexander, 1993; Theel et al., 2020).

A sedimentary facies consists of a body of rock characterised by a particular combination of lithology with physical and/or biological structures, which is distinct from the surrounding rock bodies (Walker, 1992). Facies can be combined into facies associations, which refer to a group of facies genetically related to one another, recording some environmental significance (Collinson, 1969), and into architectural elements, which are morphological subdivisions characterised by distinct geometries and depositional processes (Walker, 1992). Identifying facies, facies associations, and architectural elements provides a framework for understanding the depositional system that formed the aquifer (Walker, 1992) and, consequently, improves the predictability of subsurface heterogeneity framework (Posamentier et al., 1988; Davis et al., 1993).

The influence of sediment grain size on the permeability of granular aquifers and reservoirs has been widely discussed (e.g., Høimyr et al., 1993; Rogiers et al., 2013; Scorgie et al., 2021; Wadsworth et al., 2023), and grain size is commonly used as the main criterion to differentiate hydrostratigraphic units (Koltermann and Gorelick, 1996), particularly when notable differences in grain sizes are present, such as between mud (e.g. floodplain), sand, and gravel (e.g. Fogg et al., 1998; Klise et al., 2009; Zhang et al., 2013; Bianchi and Zheng, 2016; Soltanian et al., 2017; Yin et al., 2023; Aghaei et al., 2024). Consequently, slight variations between sand-grade facies are often overlooked, even though studies have shown that reservoir quality differences can arise due to factors such as sediment sorting, composition (e.g., Bjørlykke, 2014; Makeen et al., 2021; Castro et al., 2023; Bofill et al., 2025) and the geometry of stratification (e.g., Chandler et al., 1989; Fischer et al., 2007; Bofill et al., 2025).

Such factors are common in sand-dominated fluvio-aeolian systems, where significant permeability variations related to decimetre- to metre-scale sedimentary architecture exist, even in the absence of conglomerates or mudstone layers (e.g., Chandler et al., 1989; Fischer et al., 2007; Newell et al., 2019; Thompson et al., 2019; Bofill et al., 2025). This is the case for the upper ~70 m of the Lower Grès Vosgien formation (LGV) in eastern France (Théobald et al., 1967) (Fig. 1).

The LGV is the thickest formation of the Buntsandstein Group (Gp).

(Duringinger et al., 2019), which is a major aquifer known for hosting lithium-rich geothermal brines (Aquilina et al., 1997; Kölbel et al., 2020; Fries et al., 2022) in the Upper Rhine graben (URG) (Fig. 1), considered as one of the most favourable lithium reserves in Europe (Sanjuan et al., 2022). Additionally, to the west of the Vosges Mountains in eastern France, it is an important shallow aquifer for supplying public drinking water (Vaute et al., 2007; 2013) (Fig. 1B).

The interval studied here (Fig. 1B and 1C) is often referred to as a single, homogeneous, “fluvio-aeolian” hydrostratigraphic unit. However, Bofill et al. (2024; 2025) revealed, based on outcrops along the Vosges Mountains and core samples from the geothermal well EPS1 in the URG (Fig. 1A), several laterally extensive decimetre-scale low-permeability intervals (Fig. 1C). These, together with other heterogeneities controlled by sedimentary processes and diagenesis, occur without significant changes in clast granulometry. Such a lack of correlation between reservoir quality and grain size has also been observed in other Lower Triassic formations within the Upper Rhine Graben (e.g. Busch et al., 2025), highlighting the potential preferential flow pathways and flow baffles that have not yet been considered in flow predictions.

These findings underscore the presence of *unknown unknowns* (sensu Hunt and Welter, 2010) in hydrogeological conceptual models in the URG, as decimetre- to metre-scale heterogeneity can, for instance, considerably affect the heat recovered from an aquifer thermal energy storage (Ferguson, 2007; Possemiers et al., 2015), and the solute transport behaviour in porous aquifers (e.g. Zheng and Gorelick, 2003; Ramanathan et al., 2010; Ren et al., 2022).

The circulation of geothermal brines in the Buntsandstein Gp. is controlled by a matrix-fracture porosity system (Haffen et al., 2013; Glaas et al., 2025), with the highest permeabilities and currently economically viable geothermal loops associated with fractures (Vidal and Genter, 2018). Nevertheless, away from major fault zones, where fracture density is low (Haffen et al., 2013), fluid circulation is expected to occur primarily through sedimentary beds (Haffen et al., 2013; Heap et al., 2017).

Matrix permeability within the Buntsandstein Gp., measured from core samples obtained from geothermal wells, is highly variable, ranging from ~$10^{-4}$ mD to ~$10^3$ mD (Haffen et al., 2013; Heap et al., 2017; Bofill et al., 2025), and can be strongly controlled by sedimentary features, such as facies and facies associations (Bofill et al., 2025). Considering that (a) even a few high-transmissivity zones can be responsible for most water circulation (e.g. Jung, 2013); (b) thermal dispersivity is directly influenced by the subsurface heterogeneity architecture (e.g. Aghaei et al., 2024); and (c) in geothermal doublets, injected cold water migrates along preferential pathways, reducing thermal breakthrough time and affecting predictions of geothermal reservoir lifetime (e.g. Sommer et al., 2013; Di Dato et al., 2022; Molnár et al., 2025), it is vital to determine the location and structure of these zones and to understand how they influence fluid flow, solute transport, and heat transfer within the reservoir. This is particularly true for the LGV, which has been characterised as a favourable zone for fluid circulation within the Buntsandstein Gp. matrix due to its relatively high permeability and porosity (Haffen et al., 2013).

Understanding these flow and transport mechanisms is important not only for geothermal energy production but also for emerging applications such as lithium exploration. Comprehending the circulation of lithium-rich geothermal brines within both the matrix porosity and fractures is essential for assessing reservoir volume and estimating the available lithium resource. Lithium production from geothermal brines involves reinjecting lithium-depleted water into the reservoir and extracting lithium-rich brine through a production well. This process gradually lowers lithium concentrations in the reservoir over time (Dugamin et al., 2021) and requires accurate prediction of thermal and chemical breakthroughs (O'Sullivan et al., 2023). The greater the reservoir heterogeneity, the longer it takes to extract the maximum lithium content (Ghergut et al., 2023). Although studies suggest that,
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** (A) Location and geological map showing the studied outcrop situated in the Vosges Mountains, eastern France. It also shows the proximity between the outcrop and the geothermal well EPS1. (B) Schematic geological cross-section and stratigraphic column of the Upper Rhine Graben (URG), indicating zones where the Buntsandstein Group hosts lithium-rich geothermal brine in deep aquifers, and where it serves as a source of public water supply from shallow aquifers; (C) Conceptual model for the upper 70 m of the Lower Grès Vosgien formation, based on the geothermal well EPS1 in eastern France. The model shows an upward increase in the frequency of wind- and water-laid facies association (WWL) intervals. St = trough cross-bedded sandstone; Ss = sigmoidal cross-bedded sandstone; Sl = low-angle cross-bedded sandstone; HSS = hybrid sand sheet; f.a. = facies association. Adapted from Bouly and Demassieux (1996); Guillocheau et al. (2000); Herivaux and Marechal (2019); Duringer et al. (2019) and Bofill et al. (2025).</figcaption></figure>
once a tracer test has been conducted, the reservoir can be treated as a “black box” and a response function for lithium production can be established (Behrens et al., 2022; Kölbel et al., 2023), assumptions about the solute transport and reservoir structure remain necessary in the early stages of tracer testing (Ghergut et al., 2023).

This study integrates outcrop data from sandstones deposited in fluvio-aeolian settings into a numerical groundwater flow model to evaluate the influence of the permeability heterogeneity described by Bofill et al. (2025) on groundwater flow and transport at the site scale. Specifically, it addresses the following research questions: (i) To what extent do facies and facies association control flow behaviour in the

LGV? (ii) How does the simplification of such sedimentary heterogeneities in conceptual modelling affect the numerical model results, particularly concerning (1) particle residence time distribution within the model domain; (2) longitudinal macrodispersivity; and (3) anisotropy induced during model upscaling?

## 2. Materials and methods

### 2.1. Study area

This study is conducted using one outcrop situated at the western

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Main facies in the Lower Grès Vosgien formation observed in outcrops and core samples from the well EPS1. (A) Facies trough cross-bedded sandstone (St), and low-angle cross-bedded sandstone (Sl), from the fluvial channel facies association (FC); (B) sigmoidal cross-bedded sandstone (Ss), from the FC; (C) Hybrid sand sheet (HSS), from the wind- and water-laid facies association, overlying facies St, from the FC. Black-dashed lines mark facies internal stratifications. Core width = 8 cm.</figcaption></figure>
shoulder of the Upper Rhine Graben (URG), on the Vosges Mountains, in eastern France (Fig. 1). The region exposes several outcrops of the Buntsandstein Gp., which consists of a sequence of approximately 400 m thick sedimentary deposits composed predominantly of sandstones deposited in continental environment, mainly fluvial and aeolian, during the Lower Triassic (Dachroth, 1985).

The Lower Grès Vosgien formation (LGV) corresponds to the thickest formation within the Buntsandstein Gp., with approximately 200 m thickness measured in a geothermal well in eastern France (Aichholzer, 2019) (Fig. 1B). The LGV is built up dominantly by deposits of braided fluvial systems, intercalated with levels of aeolian sediments (Ghienne et al., 1994; Bourquin et al., 2009; Bofill et al., 2024), and has been described as preserving one of the most permeable zones for fluid migration in the Buntsandstein Gp. (Haffen et al., 2013).

## 2.2. Quantitative sedimentological context

The chosen outcrop is representative of the LGV, as it displays similar characteristics to those described and quantified by Bofill et al. (2024). The authors identified, in outcrops and in core samples, a dominance of facies from the fluvial channel facies association (FC), which accounts for 93 % of the total measured stratigraphic thickness. These are intercalated with decimetre-scale, laterally extensive tabular deposits of the wind- and water-laid facies association (WWL), comprising the remaining 7 %.

The FC is primarily composed of trough cross-bedded sandstone (St; 73 %; Fig. 2A, C), low-angle cross-bedded sandstone (Sl; 13 %; Fig. 2A), and sigmoidal cross-bedded sandstone (Ss; 11 %; Fig. 2B). These facies are predominantly medium-grained sandstones, followed by coarse-grained sandstones (Fig. 3), with scattered pebbles and granules, forming decimetre-scale individual sets, interpreted as deposited under varying flow regime conditions within braided fluvial channels (Bofill et al., 2024).

The WWL includes two element types: aeolian dunes and hybrid sand sheets (HSS). Although these elements are classified as architectural elements (sensu Miall, 1988), they are referred to as facies in this study. Because the studied outcrop exposes only a limited occurrence of aeolian dunes, only the HSS facies is considered here (Fig. 2C), representing 3 % of the LGV formation. The HSS facies consists mainly of fine-grained sandstones, with subordinate medium-grained sandstones (Fig. 3), arranged in tabular layers with an average thickness of ~30 cm (Figs. 2C and 3), and a lateral extent exceeding the maximum width of the individual outcrops (> 100 m). These deposits were formed during periods when fluvial activity was absent, and sediment transport and deposition were controlled by wind action, fluctuations in the water-table level, and ephemeral flood events (Bofill et al., 2024).

The WWL is readily identifiable in both outcrop and core samples of the LGV because, in addition to the characteristics described above, it typically displays a horizontal to subhorizontal lower boundary marked by a lag deposit composed of pebbles and granules and exhibits a bleached appearance (Fig. 2C). More detailed data on granulometry and individual set thickness by facies are presented in Fig. 3.

## 2.3. Permeability and hydraulic conductivity data

This study uses permeability data measured by Bofill et al. (2025) from outcrops using the TinyPerm3® air permeameter. Although less precise than typical laboratory methods for determining absolute permeability, air permeameters enable rapid, autonomous field measurements (Farquharson et al., 2016). This makes them highly effective for generating large datasets used to assess spatial variability in aquifer heterogeneity (e.g. Huysmans et al., 2008; Possemiers et al., 2012; Heidsiek et al., 2020). Since this study aims to evaluate the effect of relative differences in hydraulic conductivity (K) between sedimentary features in groundwater flow models, rather than to determine absolute K values, the TinyPerm 3® is well-suited for this application.

Permeability data from outcrops revealed that within the fluvial channel facies association (FC), the St, Ss, and Sl facies recorded median permeability values of 2223 mD, 1425 mD, and 554 mD, respectively, even though they generally display similar granulometry, with a slightly higher occurrence of coarse-grained sandstones in facies St (Fig. 3). Within the wind- and water-laid facies association (WWL), the hybrid sand sheet (HSS) recorded a median permeability of 56 mD (Fig. 3).

Core samples record permeability values between two and three orders of magnitude lower than those from outcrops (Fig. 3), primarily due to the higher degree of compaction, resulting from the additional overburden within the URG since the Eocene (Busch et al., 2022), and the presence of carbonate cements in the core, which were likely dissolved from the outcrops by meteoric water (Busch et al., 2022). Nevertheless, similar facies controls on permeability are observed, with progressive decreases from facies St to Ss, Sl, and HSS, with median values of 17.7 mD, 6.8 mD, 0.8 mD, and 0.11 mD, respectively (Fig. 3). Furthermore, greater permeability contrasts between facies in core samples (Fig. 3) suggest that the sedimentological controls on permeability heterogeneity are not only present in the subsurface but may also be amplified.

The permeability dataset used in this study comprises only outcrop measurements, as these provide the more extensive dataset tied to sedimentological characteristics, and offer a greater statistical representation than the core measurements. In total, 1070 measurements were used, with 972 corresponding to the fluvial channel facies association (FC) and 98 to the wind- and water-laid facies association (WWL). The compiled dataset from Bofill et al. (2025) is presented in Fig. 3.

The permeability differences are attributed to factors other than grain size and sediment sorting. Variations in the primary composition of the facies, such as the presence of ductile metamorphic rock fragments and mechanically infiltrated clay, directly influenced the diagenetic processes they underwent, particularly the degree of compaction.

Permeability (k) [m²] is fluid-independent and reflects the properties of the porous media (Loll et al., 1999). Therefore, to simulate groundwater flow, the data were converted to hydraulic conductivity (K) [mS⁻¹] using the empirical equation proposed by Iversen et al. (2003), which transforms air permeameter-based data (k)[m²] from the field into K values for saturated granular aquifers:

$$K = \frac{10^{1.22} \times \log_{10}(k \times 9.86923e^{-16}) + 13.93}{86400} \qquad (\text{Eq. 1})$$

## 2.4. From outcrop to discretised conceptual hydraulic model

An orthomosaic image of an outcrop measuring 100 m in length and 20.6 m in height—the Falkenstein Castle ruin—was created by importing drone images (with approximately 70 % overlap) into Agisoft Metashape® (Fig. 4A). The orthomosaic was used to manually delineate polygons representing zones dominated by distinct facies and facies associations (Fig. 4B). The polygons were then imported into FEFLOW 8.0®, enabling a realistic 2D vertical (cross-section) representation of the aquifer, incorporating sedimentary heterogeneities at a decimetre scale (Fig. 4). Each polygon was classified according to its corresponding facies, which FEFLOW 8.0® assigned to the mesh elements.

A triangular mesh was generated, consisting of 230,406 elements in total. Local mesh refinement was applied at the boundaries of most facies, particularly those with complex geometries or with thicknesses of only a few decimetres (Fig. 4C). To improve mesh quality, a built-in mesh smoothing tool in FEFLOW® was used. Higher mesh quality corresponds to a narrower distribution of internal triangle angles, which enhances numerical stability. Optimal mesh quality is achieved when the maximum internal angles approach 60° (Diersch, 2014).

Due to the limited number of porosity measurements from thin sections and petrophysical laboratory analyses reported by Bofill et al. (2025) and the lack of a clear understanding of porosity distribution among facies, a single representative value of 15 %, based on their data,
<figure><img src="image_6.png" /><figcaption>Fig. 3. Sedimentary data regarding facies association and their respective facies, adapted from Bofill et al. (2024); and permeability data per facies from Bofill et al. (2025). St = trough cross-bedded sandstone; Ss = sigmoidal cross-bedded sandstone; Sl = low-angle cross-bedded sandstone; HSS = hybrid sand sheet. (1mD=9.869233 × 10⁻¹⁶ m²).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** (A) Orthomosaic of the outcrop at Falkenstein Ruin. (B) Interpretation of zones dominated by individual facies, as described by Bofill et al. (2024): St = trough cross-bedded sandstone; Ss = sigmoidal cross-bedded sandstone; Sl = low-angle cross-bedded sandstone; HSS = hybrid sand sheet. (C) Finite element triangular mesh generated using FEFLOW 8.0®, with boundary conditions assigned to induce horizontal (x) and vertical (z) flow.</figcaption></figure>

was applied throughout the model domain.

## 2.5. Scenarios of heterogeneity representation

To assess the influence of sedimentary features on transport model results, six modelling scenarios were developed, each differing in sedimentary discretisation and hydraulic conductivity (K) assignment method (Table 1). These scenarios range from a fully heterogeneous model with facies-based stochastic K assignment (Scenario I, considered the reference model since it is the most realistic representation of the aquifer; Fig. 5A) to a completely homogeneous model (Scenario VI), in which a single, deterministic, weighted median K value is applied across the entire model domain. The intermediate scenarios (II to V) progressively reduce the level of sedimentary heterogeneity detail, either by merging facies or facies associations into undifferentiated units or by assigning deterministic K values based on median or weighted statistics derived from field data. The weighted K values are calculated using facies proportions reported by Bofill et al. (2024) (Fig. 3; Table 2). Fig. 5 illustrates all scenarios of heterogeneity, except for the homogeneous

model (Scenario VI).

### 2.5.1. Stochastic K assignment

The K dataset was log-transformed and then distributed stochastically using the Monte Carlo sampling method. This approach involved randomly sampling a specified number of data points (i.e. the number of mesh elements within each facies) based on a given mean and standard deviation. To avoid sampling extreme values not observed in the field data, minimum and maximum bounds were applied using a truncated normal distribution (Burkardt, 2023), implemented via the `truncnorm` module in Python's `scipy.stats` library.

Despite the log transformation, the original field K data remained negatively skewed, whereas the Monte Carlo method produced normally distributed datasets. This discrepancy led to median and mean values in the sampled data that were consistently lower than those of the original dataset. To address this issue and ensure that the median values from the Monte Carlo sampling were not significantly lower than those of the original dataset, a trial-and-error approach was used. Specifically, extremely low values were systematically removed from the original
<table><caption>**Table 1**
Scenarios of heterogeneity representation by geological discretisation scale and hydraulic conductivity (K) assignment.</caption><thead><tr><th>Scenario</th><th>Sedimentary Discretization</th><th>K Assignment Method</th><th>Heterogeneity Level</th></tr></thead><tbody><tr><td>I</td><td>Facies geometries</td><td>Stochastic distribution using field K data within facies</td><td>High (between facies)</td></tr><tr><td>II</td><td>Facies geometries</td><td>Single median K per facies</td><td>Medium (between facies)</td></tr><tr><td>III</td><td>Facies associations (FC/ WWL)</td><td>Stochastic distribution, weighted by facies proportions per association</td><td>Medium-low (between associations)</td></tr><tr><td>IV</td><td>Facies associations (FC/ WWL)</td><td>Median K weighted by facies proportion per association</td><td>Low (between associations)</td></tr><tr><td>V</td><td>No sedimentary geometry</td><td>Stochastic, weighted by facies proportions</td><td>Very low (global variability only)</td></tr><tr><td>VI</td><td>No sedimentary geometry</td><td>Single median K for the entire domain weighted by facies proportions</td><td>None (homogeneous)</td></tr></tbody></table>

dataset before sampling. This adjustment aimed to achieve a minimum ratio of 0.9 between the median of the Monte Carlo sample and that of the original dataset (e.g. see Tables 3 and 4).

For facies HSS, which exhibited the highest variance of K among all facies (Table 3) and maintained substantial deviation from a normal distribution even after log transformation, a lower ratio of 0.72 (Tables 3 and 4) was accepted. This compromise prevented excessive data exclusion solely to match the medians. Depending on the facies and the degree of skewness, this process resulted in the removal of between 5% (approximately two standard deviations below the mean) and 9% (for facies HSS) of the original data.

This approach was essential for generating a Monte Carlo-sampled K dataset that more accurately reflected field observations. It also underscores the limitations of applying log transformations without subsequent evaluation to confirm that the transformed data realistically represents the underlying field data.

### 2.5.2. Deterministic K assignment

In the deterministic scenarios (II, IV, and VI), median field K values were assigned to the mesh elements, constrained by the geometry of sedimentary features. In Scenario IV, values were weighted based on the proportion of facies within each facies association (weighted K in Table 2, and Ko in Table 4; Fig. 5D). In Scenario VI, a single K value of $4.8 \times 10^{-6}$ m/s was assigned to the entire model domain, weighted based on the proportion of facies within the entire LGV (Fig. 3). Median values are preferred when representing non-normally distributed data, as they are less affected by extreme values compared to the mean.

## 2.6. Fluid-flow and particle tracking modelling

Fluid-flow simulations were conducted in a two-dimensional vertical cross-section under saturated, steady-state flow conditions using FEFLOW 8.0®. This software employs a finite element method, capable of handling complex geometries in groundwater flow and solute transport models (Diersch, 2014). The simulations were governed by the law of conservation of mass:

$$S_s \frac{\partial h}{\partial t} = \frac{\partial}{\partial x_i} \left( K_{ij} \frac{\partial h}{\partial x_j} \right) - Q \qquad (2)$$

where, $S_s$ is the specific storage (set to zero in this study due to steady-state conditions); $h$ is the hydraulic head; $t$ is time, $K_{ij}$ is the hydraulic conductivity tensor; $x$ is the spatial coordinate; $i$ and $j$ are the coordinate directions; and $Q$ represents local sources or sinks (also set to zero in this study).

A fixed hydraulic head boundary condition was applied uniformly along both lateral boundaries to induce predominantly horizontal groundwater flow. The hydraulic gradient was set at 0.018, based on the potentiometric map of the Lower Triassic aquifer from the hydrogeological atlas of the Rhine-Meuse basin (BRGM, 2002), inducing a flow from left (inlet) to right (outlet) in the model domain represented in Fig. 4C. To simulate vertical flow, perpendicular to the sedimentary bedding, fixed hydraulic head boundaries were applied to induce flow from the top (inlet) to the bottom (outlet) of the model domain (Fig. 4C), with the same hydraulic gradient of 0.018.

Where no boundary conditions were defined, FEFLOW 8.0® assigned no-flow (impermeable) boundaries by default, meaning that no water could enter or leave through these limits (Fig. 4C).

Forward particle tracking was employed to visualise flow pathlines and estimate residence times across the model domain. This post-processing computational method traces the advective trajectories of massless particles introduced at specific starting points within a known velocity field (Diersch, 2014). Particles were introduced at all mesh elements along the upstream boundary, resulting in 369 particles for each horizontal flow simulation and 213 for the vertical flow simulations.

To ensure representative coverage of possible flow paths, multiple simulation runs of particles were conducted. These were included in the analysis until the cumulative mean residence time stabilised (see Supplementary Material, Fig. S1). This iterative approach was particularly important for scenarios with stochastically distributed K values, where complex heterogeneity creates a wide range of potential particle trajectories across the model.

## 2.7. Quantification of upscaled anisotropy

Given the relatively small dimensions of the model domain (100 m in length and 20.6 m in height), which are comparable to the size of a single cell in regional-scale groundwater models, it is essential to evaluate how permeability heterogeneity translates into upscaled anisotropy. To do this, anisotropy was assessed after each simulation by calculating upscaled hydraulic conductivity ($K_{up}$) in both horizontal and vertical directions.

During each simulation run, FEFLOW® automatically computes the average flow rate $Q$ ($m^3 \cdot d^{-1}$) across the model domain. Using this output result, $K_{up}$ was calculated based on Darcy's law:

$$K_{up} = \frac{Q \Delta L}{A \Delta h'} \qquad (3)$$

where A is the cross-sectional area ($m^2$): for horizontal flow, 1 m (width) $\times$ 20.6 m (outcrop height), and for vertical flow, 1 m (width) $\times$ 100 m (outcrop length) for vertical flow; $\Delta h$ is the hydraulic head difference between the inflow and outflow boundaries; and $\Delta L$ is the length of the model in the considered flow direction.

The upscaled anisotropy was then calculated as the ratio of vertical to horizontal $K_{up}$:

$$Upscaled anisotropy = \frac{K_{up \text{ vertical}}}{K_{up \text{ horizontal}}} \qquad (4)$$

This approach enables the comparison of how each heterogeneity scenario affects anisotropy at a scale more relevant to regional groundwater models.

## 2.8. Quantification of macrodispersivity using temporal moments analysis

In this study, the temporal moment approach is applied to analyse the residence time of water particles across various model scenarios. Temporal moment analysis is a well-established method for interpreting breakthrough curves obtained from tracer tests in porous aquifers (e.g., Kreft and Zuber, 1978; Valocchi, 1985; Schäfer, 1991). Each particle is assumed to represent a fixed, constant mass of tracer relative to a
<figure><img src="image_8.png" /><figcaption>Fig. 5. Scenarios of heterogeneity representation by geological discretisation scale and hydraulic conductivity (K) assignment.</figcaption></figure>

specific water volume, and all particles are introduced simultaneously at the model inlet. The distribution of particle arrival times at the model outlet is expressed as a frequency, which is interpreted as the tracer concentration breakthrough curve $c(x_0, t)$, where $x_0$ denotes the outlet location.

The temporal moments of a concentration distribution $c(x, t)$ are defined following Aris (1958) as:

$$M_{n,t} = \int_{0}^{\infty} t^n c(x, t) dt, \qquad (5)$$

where, $M_{n,t}$ represents the n-th temporal moment around the origin ($t = 0$). Analogous to the spatial moment method, the normalised temporal

moments $m_{n,t}$ is calculated. The first normalised moment, $m_{1,t}$, corresponds to the mean residence time of the particles:

$$m_{1,t} = \frac{M_{1,t}}{M_{0,t}} \qquad (6)$$

The zeroth moment $M_{0,t}$ refers to the total tracer mass injected at a constant flow rate at the model inlet. In this study, this corresponds to the total number of particles used in the simulation. Given that the number of particles varies across scenarios (Fig. 7), and to ensure all scenarios are normalised to the same initial concentration, $M_{0,t}$ was set as 1, representing 100% of the particles. Accordingly, $c(x, t)$ is expressed as the percentage of particles arriving at the model outlet over time.

For an instantaneous tracer injection (a Dirac impulse), the mean flow velocity ($U$) (m.d$^{-1}$) across the domain length $x_0$ is calculated as
**Table 2**

Weighted K values calculated using original field data (Ko) and facies proportion normalised from Bofill et al. (2024), for assigning K values at the facies association scale (Scenarios III and IV).

<table><thead><tr><th>Facies association</th><th>Facies</th><th>Ko (ms<sup>-1</sup>)</th><th>Normalised % of the facies association thickness (Bofill et al., 2024)</th><th>Weighted K (ms<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="3">Fluvial channel</td><td>St</td><td>5.9 × 10<sup>-6</sup></td><td>74 %</td><td rowspan="3">4.9 × 10<sup>-6</sup></td></tr><tr><td>Ss</td><td>3.4 × 10<sup>-6</sup></td><td>12 %</td></tr><tr><td>Sl</td><td>1.1 × 10<sup>-6</sup></td><td>14 %</td></tr><tr><td>Wind- and water-laid</td><td>HSS</td><td>6.6 × 10<sup>-8</sup></td><td>In the studied outcrop, it is considered 100 % HSS</td><td>6.6 × 10<sup>-8</sup></td></tr></tbody></table>

**Table 3**

Hydraulic conductivity (K) data used in Scenarios I and II, by individual facies. Original median hydraulic conductivity (Ko) measured in the field, and median hydraulic conductivity of the stochastically generated data by Monte Carlo sampling (Kmc). σ²(o) is the variance of the original dataset in the log domain; σ²(mc) is the variance of the stochastically generated data in the log domain. Kmc/Ko indicates how much the simulated data approximates the real data.

<table><thead><tr><th>Facies</th><th>Ko [ms<sup>-1</sup>]</th><th>σ<sup>2</sup>(o)[m<sup>2</sup> s<sup>-2</sup>]</th><th>Kmc [ms<sup>-1</sup>]</th><th>σ<sup>2</sup>(mc)[m<sup>2</sup> s<sup>-2</sup>]</th><th>Kmc/Ko</th></tr></thead><tbody><tr><td>St</td><td>5.9 × 10<sup>-6</sup></td><td>0.69</td><td>5.8 × 10<sup>-6</sup></td><td>0.19</td><td>0.98</td></tr><tr><td>Ss</td><td>3.4 × 10<sup>-6</sup></td><td>0.48</td><td>3.1 × 10<sup>-6</sup></td><td>0.28</td><td>0.90</td></tr><tr><td>Sl</td><td>1.1 × 10<sup>-6</sup></td><td>0.66</td><td>1.1 × 10<sup>-6</sup></td><td>0.47</td><td>0.99</td></tr><tr><td>HSS</td><td>6.6 × 10<sup>-8</sup></td><td>1.4</td><td>4.8 × 10<sup>-8</sup></td><td>0.67</td><td>0.72</td></tr></tbody></table>

**Table 4**

Original weighted median hydraulic conductivity (Ko) measured in the field, and the median hydraulic conductivity of stochastically generated data by Monte Carlo sampling (Kmc), per facies association. The Kmc/Ko ratio indicates the degree to which the simulated data approximates the field data in Scenario III.

<table><thead><tr><th>Facies association</th><th>Ko (ms<sup>-1</sup>)</th><th>Kmc (ms<sup>-1</sup>)</th><th>Kmc/Ko</th></tr></thead><tbody><tr><td>Fluvial channel</td><td>4.9 × 10<sup>-6</sup></td><td>4.9 × 10<sup>-6</sup></td><td>1.0</td></tr><tr><td>Wind- and water-laid</td><td>6.6 × 10<sup>-8</sup></td><td>4.7 × 10<sup>-8</sup></td><td>0.72</td></tr></tbody></table>

(Kreft and Zuber, 1978):

$$U = \frac{x_0}{m_{1,t}} \qquad (7)$$

Higher-order temporal moments are typically represented as central moments ($\mu_{n,t}$), computed by:

$$\mu_{n,t} = \frac{\int_{0}^{\infty} (t - m_{1,t})^{n} c(x, t) dt}{M_{0,t}} \qquad (8)$$

The second central moment ($\mu_{2,t}$), corresponds to the variance of the concentration distribution. Under Dirac impulse injection conditions, the longitudinal dispersion coefficient ($D_L$) (m²·d⁻¹) is calculated as (Kreft and Zuber, 1978):

$$D_L = \frac{\mu_{2,t} U^3}{2x_0} \qquad (9)$$

Finally, the longitudinal dispersivity ($\alpha_L$) (m) is determined from the dispersion coefficient and mean flow velocity as:

$$\alpha_L = \frac{D_L}{U} \qquad (10)$$

## 3. Results

### 3.1. Particle tracking - horizontal flow

Histograms of particle residence times for each scenario reveal variations in transport behaviour, particularly related to peaks and sub-peaks of frequency of particle residence time, skewness and kurtosis (Fig. 6). Fig. 6 and Table 5 summarise the key metrics, such as statistics of particle arrival times and upscaled horizontal K, across all six scenarios. Detailed interpretations are provided in the following subsections.

#### 3.1.1. Interpretation of particles’ arrival times

The higher permeability, dominant proportion, and inter-connectivity of facies St cause it to act as the preferential flow path, controlling the earliest particle arrivals. In contrast, the relatively thin, low-permeability, and laterally extensive facies HSS are interpreted as the primary particle retention zones, contributing to the long-tailing observed in the late particle arrivals (e.g. Schäfer, 1991; Zheng and Gorelick, 2003; Klise et al., 2009; Bianchi et al., 2011; Di Dato et al., 2022; Yin et al., 2023).

In Scenario I (Fig. 7A), where K values are distributed stochastically, the first peak is interpreted as corresponding to particles that travelled primarily through facies St, and the second peak to particles that were more affected by the lower-permeability zones of facies Sl. Facies Ss acts as an intermediate unit, where particles that travelled through it may contribute to both peaks, but are predominantly associated with the first peak. This interpretation is supported by Bofill et al. (2025), who confirm that the permeability values of St and Ss do not represent statistically distinct populations, whereas Sl constitutes a separate, statistically distinct facies.

In Scenario II (Fig. 6B), the two main peaks (modes) are sharper than the peaks observed in Scenario I, which is attributed to the deterministic assignment of median K values by facies. As observed in Scenario I, the first peak is linked to facies St with some influence from Ss. The transition between the peaks corresponds to Ss, while the second peak is associated with Sl. A third, smaller peak, occurring at 6.77 years, is interpreted as being linked to particles that were injected in the inlet in facies HSS, but relatively soon after injection migrated out of this low-permeability zone and were incorporated into the facies with higher permeability from the FC facies association.

In Scenario III (Fig. 6C), the single main peak is interpreted as resulting from the grouping of the three FC facies (St, Ss, and Sl) into a single hydrostratigraphic unit. The long, positively skewed tail is attributed to the WWL association (facies HSS), which is consistent with the absence of tailing in Scenarios V and VI, where WWL is not explicitly represented.

In Scenario IV (Fig. 7D), which represents only the two facies associations captured deterministically, a secondary peak occurs at 6.77 years, similar to that observed in Scenario II. The long tail, extending to 15.4 years of residence time, is attributed to particles strongly influenced by facies HSS. In Scenarios V and VI (Figs. 6E and 6F), where HSS is not represented, no particles arrive after 6 years. This absence reinforces the link between the second peak in Scenarios II and IV, as well as the long tail in Scenarios I to IV, and facies HSS.

The behaviour of the cumulative residence time curves (Fig. 7) reflects that of the histograms. The peaks observed in the histograms of stochastically distributed K models (Figs. 6A, 6C, and 6E) are smoothed in the cumulative curves, whereas the peaks in the deterministic K models (Figs. 6B, 6D, and 6F) are recognisable as steps in the cumulative curves. The most evident differences of cumulative residence time curves are again in Scenarios V and VI (Fig. 7), aligned with observations made by Bianchi and Pedretti (2017), where long-tailing emerges in more ordered structures (sedimentary geometry), while a disordered random distribution of K values results in less skewed residence time distribution (e.g. Fig. 6E).
<figure><img src="image_9.png" /></figure>

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

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

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

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

<figure><img src="image_14.png" /><figcaption>Fig. 6. Histograms of residence times for vertical flow particles across the different modelling scenarios.</figcaption></figure>

<table><caption>**Table 5**
Compilation of results obtained from the horizontal flow.</caption><thead><tr><th>Scenario</th><th>Sedimentary Heterogeneity</th><th>K assignment</th><th>Total of particles analysed</th><th>Mean residence time (years)</th><th>Median residence time (years)</th><th>First particle to arrive (years)</th><th>Last particle to arrive (years)</th><th>Macro-dispersivity (m)</th><th>Upscaled K (ms<sup>-1</sup>)</th></tr></thead><tbody><tr><td>I</td><td>Facies scale</td><td>Stochastic</td><td>16227</td><td>5.64</td><td>4.92</td><td>3.65</td><td>27.1</td><td>9.15</td><td>5.44 × 10<sup>-6</sup></td></tr><tr><td>II</td><td>Facies scale</td><td>Deterministic</td><td>381</td><td>6.25</td><td>5.51</td><td>4.49</td><td>15.4</td><td>7.26</td><td>4.79 × 10<sup>-6</sup></td></tr><tr><td>III</td><td>Facies association scale</td><td>Stochastic</td><td>6586</td><td>5.64</td><td>4.7</td><td>4.16</td><td>31.9</td><td>11.35</td><td>5.52 × 10<sup>-6</sup></td></tr><tr><td>IV</td><td>Facies association scale</td><td>Deterministic</td><td>337</td><td>6.38</td><td>5.48</td><td>5.37</td><td>14.8</td><td>4.78</td><td>4.65 × 10<sup>-6</sup></td></tr><tr><td>V</td><td>No sedimentary geometry</td><td>Stochastic</td><td>8277</td><td>5.15</td><td>5.06</td><td>4.63</td><td>5.9</td><td>0.03</td><td>5.26 × 10<sup>-6</sup></td></tr><tr><td>VI</td><td>No sedimentary geometry</td><td>Deterministic</td><td>368</td><td>5.57</td><td>5.49</td><td>5.48</td><td>5.49</td><td>0</td><td>4.82 × 10<sup>-6</sup></td></tr></tbody></table>

### 3.1.2. Macrodispersivity

The longitudinal dispersivities calculated based on Eq. (10) for the different model scenarios are hereafter called longitudinal macro-dispersivities as they result from the averaging process of the individual advective solute contributions of each particle to the overall mass

transport. The value obtained for reference Scenario I was 9.1 m This result aligns with the assumption proposed by Gelhar et al. (1992), who, through a compilation of numerous field experiments, suggested that longitudinal macrodispersivity scales with travel distance and is approximately one-tenth of it. Given the model length of 100 m, the
<figure><img src="image_15.png" /><figcaption>**Fig. 7.** Cumulative curves of horizontal flow particles' residence time across the distinct scenarios.</figcaption></figure>

value obtained in this study is consistent with field observations.

The most evident effect of model simplification is observed in Scenarios V and VI, where the two decimetre-scale levels of wind- and water-laid facies associations (WWL) are not represented. In these cases, longitudinal macrodispersivity due to the more homogeneous flow field remained close to zero (Table 5). Studies have shown that higher macrodispersivity is expected not only where there are greater contrasts in

facies permeability, but also where the average length of the facies is larger (Dong et al., 2018). This applies to the WWL, whose lateral extent exceeds the outcrop dimensions.

Models with deterministically distributed hydraulic conductivity (K) consistently showed lower macrodispersivity compared to scenarios with the same level of heterogeneity representation (e.g., Scenarios I vs. II, and III vs. IV; Table 5), which is expected considering that macrodispersion arises from hydraulic conductivity variability (e.g., Jussel, 1989; Schäfer, 1991).

It is important to note that the so-called asymptotic regime of macrodispersion was not reached in all heterogeneous scenarios. Breakthrough curves that exhibit a multi-step increase in concentration (Scenarios II and I; Fig. 7) and/or display long-tailing (Scenarios I to IV; Figs. 6 and 7) are indicative of a pre-asymptotic regime (Talon et al., 2023), where macrodispersivity continues to increase with time and travel distance (Gelhar et al., 1992), and cannot be described by Fickian dispersion concepts (Herr et al., 1989).

### 3.1.3. Comparison between scenarios

No notable differences are observed in the mean or median values of particle residence time across the scenarios (Figs. 6 and 8), which directly influence their similar upscaled horizontal K (Table 5). However, differences become more evident when analysing the distribution of particle arrival times at the model outlet (Figs. 6 and 8).

Scenario I shows the earliest particle arrivals among all models, at 3.7 years (Figs. 6A and 8; Table 5), interpreted as primarily due to the

<figure><img src="image_16.png" /><figcaption>**Fig. 8.** Hypothetical plume migration differences based on arrival times at the model outlet, using Scenario I as a reference.</figcaption></figure>
discretisation of facies St, which acts as a preferential flow zone because of its relatively high K and interconnectivity (Fig. 9) (e.g. Yin et al., 2023). Although facies St is also discretised in Scenario II, assigning a single deterministic K (median) removes mesh elements with K values higher than the median, thereby affecting the preferential pathways. As the complexity of sedimentary geometries is progressively simplified in subsequent scenarios, high-K mesh elements are increasingly homogenised or are lost entirely. This leads to delayed particle arrivals and a systematic increase in minimum residence times, causing failure to replicate the earliest particle arrivals observed in Scenario I (Figs. 6 and 8).

The last-arriving particles also showed differences across the scenarios (Figs. 6 and 8; Table 5). Scenarios I and III exhibited the longest tailing in particle travel times (Figs. 6A, 6C, 7, and 8), mainly due to the discretisation of the HSS facies and the stochastic assignment of K values. As with the first-arriving particles, where mesh elements with higher K values enabled faster movement to the model outlet, those passing through elements with lower K values, represented at the lower end of the negatively skewed K distribution of facies HSS (Bofill et al., 2025), experienced greater resistance to flow. Therefore, HSS facies acted as the main particle retention zones within the aquifer (Fig. 9) (e.g. Yin et al., 2023). The resulting slower movement contributed to the increased positive skewness observed in the residence time histograms for these two scenarios (Figs. 6A and 6C).

Although Scenario III recorded the highest overall residence time, only 0.2 % of its particles travelled longer than the last-arriving particle in Scenario I (27.1–31.9 years). Since both scenarios used the same statistical distribution of K to represent the low-permeability HSS facies, this minor difference is considered negligible and is not discussed further. The same is valid for their deterministic counterparts, Scenarios

II and IV.

The most notable differences in the arrival times of the last particles are observed when the model simplification does not explicitly represent the main particle retention zones, the HSS facies. In these scenarios (V and VI), the K values originally associated with the layers of HSS are redistributed throughout the entire model domain. As a result, the long tail seen in Scenarios I to IV (Fig. 6) is no longer present, and the residence times of the last-arriving particles are much closer to the mean. For example, in Scenario VI, 100 % of the particles arrived within a 5.4-day (0.01 years) interval, whereas in the reference scenario, the last particle arrived 23.4 years after the first.

Based on the first and the last particles to arrive at the model outlet, the overlap between the simplified scenarios (II to VI) and the reference Scenario I was calculated. This analysis aimed to quantify how much of the actual mass concentration ($c(x, t)$) from the most realistic case (Scenario I) would be captured within the arrival intervals predicted by the simplified scenarios. Results showed that Scenario III had the highest overlap, of 78 %, while the arrival interval predicted by Scenario VI would capture only 0.5 % of the mass concentration from Scenario I (Fig. 8).

Table 6 and Fig. 8 summarise the percentage of particle concentration ($c(x,t)$) that would have already arrived in Scenario I by the time the first particle arrives in each simplified model (denoted as F in Fig. 8). They also show the portion of $c(x,t)$ from Scenario I that arrives after the last particle in each simplified scenario (denoted as L), and the overall overlap in particle residence times between each simplified model and Scenario I (denoted as O).

<figure><img src="image_17.png" /><figcaption>**Fig. 9.** Examples of particle streamlines from FEFLOW® flow simulations in Scenario I, where facies geometries were discretised, and Scenario V, where sedimentary architecture was not used to constrain hydraulic conductivity (K) data.</figcaption></figure>
<caption>**Table 6**</caption>

Percentage of particles from Scenario I that arrived before the first particles of the other scenarios; after the last particles of the other scenarios; and the overlap of particles that arrived in the same range of residence time observed in Scenario I.

<table><thead><tr><th></th><th>% of particles from Scenario I that arrived before the 1st particle arrives</th><th>% of particles from Scenario I that arrived after the last particle arrived</th><th>Overlap of particles arriving in the range or residence time observed in Scenario I</th></tr></thead><tbody><tr><td>Scenario II</td><td>36 %</td><td>1 %</td><td>63 %</td></tr><tr><td>Scenario III</td><td>22 %</td><td>0 %</td><td>78 %</td></tr><tr><td>Scenario IV</td><td>66 %</td><td>1 %</td><td>33 %</td></tr><tr><td>Scenario V</td><td>40 %</td><td>23 %</td><td>37 %</td></tr><tr><td>Scenario VI</td><td>68.8 %</td><td>30.7 %</td><td>0.5 %</td></tr></tbody></table>

## 3.2. Particle tracking - vertical flow

The histograms of residence time in the vertical flow simulations show distributions that are more symmetric than in the horizontal flow case (Fig. 10). This is reflected in skewness values that remain close to zero for most scenarios, except for Scenario VI (homogeneous). Table 7 summarises the key metrics, including particle arrival times and upscaled vertical K, across all six scenarios. Detailed interpretations are provided in the following subsections.

### 3.2.1. Interpretation of particles’ arrival times

The relatively more symmetric residence time distributions observed in the vertical flow simulations can be explained by the fact that the flow crosses the dominant high- and low-permeability layers almost perpendicularly. Under these conditions, most particles encounter similar hydraulic contrasts and follow comparable pathways, which results in less variable and more uniform residence times. In contrast, the horizontal-flow simulations allow particles to sample heterogeneity laterally, producing broader and more skewed distributions.

<figure><img src="image_18.png" /><figcaption>**Fig. 10.** Histograms of residence times for vertical flow particles across the different modelling scenarios.</figcaption></figure>
<table><caption>**Table 7**
Summary of results obtained from the vertical flow simulations, at different scenarios.</caption><thead><tr><th>Scenario</th><th>Sedimentary Heterogeneity</th><th>K assignment</th><th>Total of particles analysed</th><th>Mean residence time (years)</th><th>Median residence time (years)</th><th>First particle to arrive (years)</th><th>Last particle to arrive (years)</th><th>Upscaled K (ms<sup>-1</sup>)</th></tr></thead><tbody><tr><td>I</td><td>Facies scale</td><td>Stochastic</td><td>6063</td><td>3.5</td><td>3.4</td><td>1.4</td><td>5.6</td><td>1.64 × 10<sup>-6</sup></td></tr><tr><td>II</td><td>Facies scale</td><td>Deterministic</td><td>209</td><td>4.7</td><td>4.6</td><td>2.2</td><td>6</td><td>1.26 × 10<sup>-6</sup></td></tr><tr><td>III</td><td>Facies association scale</td><td>Stochastic</td><td>3074</td><td>3.2</td><td>3.1</td><td>1.4</td><td>5.6</td><td>1.84 × 10<sup>-6</sup></td></tr><tr><td>IV</td><td>Facies association scale</td><td>Deterministic</td><td>211</td><td>4.3</td><td>4.2</td><td>2</td><td>6</td><td>1.39 × 10<sup>-6</sup></td></tr><tr><td>V</td><td>No sedimentary geometry</td><td>Stochastic</td><td>2054</td><td>1.1</td><td>1</td><td>0.8</td><td>1.3</td><td>5.33 × 10<sup>-6</sup></td></tr><tr><td>VI</td><td>No sedimentary geometry</td><td>Deterministic</td><td>213</td><td>1.1</td><td>1.1</td><td>1.1</td><td>1.1</td><td>4.82 × 10<sup>-6</sup></td></tr></tbody></table>

### 3.2.2. Comparison between scenarios

For scenarios produced at the facies and facies association scale, deterministic scenarios delayed particles' residence time results compared to their stochastic counterparts (Figs. 10A and B; C, and D). This is attributed to the fact that, as previously discussed in the horizontal flow section, when K is represented as a median or mean (deterministic) value, mesh elements with high K are not represented, hindering preferential pathways for the particles to flow. For this reason, the most significant differences between the stochastic and deterministic scenarios occur in the minimum particle residence time (Fig. 10; Table 7).

The cumulative curves of particles' residence time of stochastically distributed K (scenarios I, III and V) showed a smooth shape, with no sharp steps as observed in the deterministic models (scenarios II, IV, and IV) (Fig. 11). The sharp steps result from the absence of alternative pathways within the mesh elements of each facies in deterministic models, which leads particles to spend similar amounts of time traversing the distinct facies.

When sedimentary geometries were gradually removed, the particles' residence time displayed an opposite behaviour to that observed for horizontal flow, particularly regarding the minimum residence time. In horizontal flow, the first particles arrive later when the model is simplified; conversely, in vertical flow, they arrive sooner. This occurs because low-permeability layers, particularly HSS, are typically horizontal, forcing vertically migrating particles to move through them without access to higher-permeability pathways. Although weighting the hydraulic conductivity reduces the overall mean K, this effect is minor compared to eliminating the necessity for particles to cross low-permeability facies.

### 3.3. Upscaled K anisotropy

The upscaled K for horizontal flow showed no significant variation

<figure><img src="image_19.png" /><figcaption>**Fig. 11.** Cumulative curves of vertical flow particles' residence time across the distinct scenarios.</figcaption></figure>

across scenarios (Table 5). In contrast, for vertical flow, the upscaled K increased by approximately three to four times when the HSS facies layer was not discretised (Scenarios V and VI; Table 7).

The results indicate that facies within the fluvial channel facies association have minimal impact on upscaled K in both horizontal and vertical flow directions. In contrast, vertical flow is notably affected by the wind- and water-laid facies association (WWL) due to its tabular geometry and broad lateral extent. Most particles must cross these low-permeability intervals, which lowers the overall upscaled K.

Regarding upscaled anisotropy, Scenarios I to IV exhibit similar upscaled anisotropy values of approximately 0.3, whereas Scenarios V and VI, where no sedimentary architecture is represented, yield isotropic upscaled K ($K_v/K_h = 1$) (Table 8). Therefore, representing or not the individual facies within the FC resulted in minor differences in the upscaled anisotropy when compared to the discretisation of the WWL intervals, regardless of distributed K data across the mesh elements stochastically or as a single deterministic median value.

## 4. Discussions

### 4.1. Heterogeneities in the Lower Grès Vosgien formation and their impact on fluid-flow and transport modelling

The importance of discretising small-scale heterogeneity in fluid flow, heat-transfer, and solute transport modelling has been previously recognised (e.g. Possemiers et al., 2015; Comolli et al., 2019; Hakoun et al., 2019). However, when these heterogeneities are not associated with notable changes in grain sizes, they are often overlooked. This is particularly relevant for the Buntsandstein Gp. in the Upper Rhine Graben, where reservoir quality is largely governed by diagenesis (Busch et al., 2022; Busch, 2025; Quandt et al., 2022).

In the Lower Grès Vosgien formation (LGV), diagenetic processes affected the facies differently (Bofill et al., 2025). For example, within the medium-grained fluvial channel facies association (FC), the facies trough cross-bedded sandstone (St) was generally less affected by compaction than the facies low-angle cross-bedded sandstone (Sl), thereby preserving higher permeability (Bofill et al., 2025). As demonstrated in Scenarios III and IV, treating these facies as a single unit resulted in underestimation of preferential flow zones and,

<table><caption>**Table 8**
Upscaled anisotropy calculated for each scenario of heterogeneity representation.</caption><thead><tr><th>Scenario</th><th>Vertical upscaled K (Kv) (ms<sup>-1</sup>)</th><th>Horizontal upscaled K (Kh) (ms<sup>-1</sup>)</th><th>Upscaled anisotropy (Kv/Kh)</th></tr></thead><tbody><tr><td>I</td><td>1.64 × 10<sup>-6</sup></td><td>5.44 × 10<sup>-6</sup></td><td>0.30</td></tr><tr><td>II</td><td>1.26 × 10<sup>-6</sup></td><td>4.79 × 10<sup>-6</sup></td><td>0.26</td></tr><tr><td>III</td><td>1.84 × 10<sup>-6</sup></td><td>5.52 × 10<sup>-6</sup></td><td>0.33</td></tr><tr><td>IV</td><td>1.39 × 10<sup>-6</sup></td><td>4.65 × 10<sup>-6</sup></td><td>0.30</td></tr><tr><td>V</td><td>5.33 × 10<sup>-6</sup></td><td>5.26 × 10<sup>-6</sup></td><td>1.01</td></tr><tr><td>VI</td><td>4.80 × 10<sup>-6</sup></td><td>4.82 × 10<sup>-6</sup></td><td>1.00</td></tr></tbody></table>
consequently, delayed first-arrival times at the model outlet (Fig. 8).

Similarly, the medium-grained sandstones of the FC were less affected by compaction than the fine-grained sandstones of the wind- and water-laid facies association (WWL), represented here as composed only of the facies hybrid sand sheet (HSS). The combination of smaller grain size and higher compaction degree in the WWL produced permeability contrasts of up to four orders of magnitude between the WWL and the FC (Bofill et al., 2025). As a result, the WWL facies acts as significant flow baffles and solute-retention zones. When WWL deposits were not explicitly discretised, as in Scenarios V and VI, late-time particle arrivals were not reproduced. In all these cases, neglecting these small-scale features directly influenced the breakthrough times at the model outlet.

## 4.2. Stratigraphic architecture and macrodispersion implications

Core data from the EPS1 well show multiple WWL intervals within the upper 70 m of the LGV, with their frequency increasing upward (Bofill et al., 2024) (Fig. 1C). On average, FC deposits preserve 4 to 5 m of thickness between WWL intervals (Bourquin et al., 2009; Bofill et al., 2024). In contrast, the outcrop used for flow simulations includes only two WWL intervals separated by approximately 13 m of FC (Figs. 4, 5, and 9). Therefore, the higher frequency and cumulative thickness of WWL intervals in the subsurface are likely to enhance their influence on macrodispersion (Zhang et al., 2007; Zhang et al., 2013).

Thicker retention zones, represented here by the HSS facies, extend the tails of breakthrough curves (Zhang et al., 2007; Zhang et al., 2013); therefore, they prolong the non-Fickian behaviour and prevent dispersion from reaching asymptotic conditions. The stronger permeability contrasts observed between facies in core samples compared with outcrops (Fig. 3) indicate that even greater macrodispersion is expected in the subsurface (Dong et al., 2018). The effects demonstrated in this study therefore represent a conservative estimate when extrapolated to deep reservoirs.

## 4.3. Anisotropy and scaling effects

Each HSS interval exposed in the studied outcrop has an average thickness of ~30 cm (e.g., Fig. 5), yielding a cumulative ~0.6 m of HSS within the 20.6 m succession, and an upscaled anisotropy of ~0.3 (Table 8). These outcrop-derived data suggest that similar anisotropic effects are likely in the upper 70 m of the LGV, where twelve FC-WWL interactions are identified between 1161 m and 1231 m depth in the EPS1 well (Fig. 1C).

Assuming an average HSS thickness of 25 cm for each interval (Fig. 3) gives a total of 3 m of HSS within this 70 m section. Using the median hydraulic conductivities (K) of FC and WWL (Scenario IV; Fig. 5D), the upscaled horizontal (Kh) and vertical (Kv) conductivities were calculated as thickness-weighted arithmetic and harmonic means (Woessner and Poeter, 2024). The resulting Kv/Kh ratio is approximately 0.25.

These anisotropy values can directly inform the parameterisation of flow models at different scales: for example, a value of ~0.3 may be representative for cells on the order of the outcrop thickness (~20 m), whereas a value near 0.25 is appropriate for coarser-scale models with 70 m vertical resolution. This estimate likely represents a lower bound, as contrasts between WWL and FC are expected to be greater in the subsurface (Fig. 3) (Bofill et al., 2025). Such observations complement laboratory-scale anisotropy measurements previously reported for the Buntsandstein Gp. (Goupil et al., 2022), and highlight the potential influence of anisotropy on large-scale fluid flow and heat transfer within the Buntsandstein Gp. matrix reservoir.

## 4.4. Implications for geothermal and lithium resources

The Buntsandstein Gp. is recognised as a major lithium reservoir in Europe (Sanjuan et al., 2022), where lithium concentrations may exhibit

significant spatial variability (Kölbel et al., 2023). Within this context, the LGV represents the thickest formation of the Buntsandstein Gp., and the studied section, characterised by fluvio-aeolian interactions, corresponds to one of the most favourable zones for fluid circulation (Haffen et al., 2013). Prior to Bofill et al. (2025), this interval was generally regarded as a single “fluvio-aeolian” unit, and the WWL intervals were not recognised as potential flow baffles.

Given the strong coupling between advective flow and heat transfer (Sauty et al., 1982; Casillas-Trasvina et al., 2022), as well as the role of advection in solute transport (de Lange, 2020), the results presented here demonstrate that such heterogeneities, particularly the WWL intervals, exert a strong control on lithium-rich geothermal brine resources circulation within the matrix porosity and, therefore, the fluvio-aeolian deposits should not be represented as a single, homogeneous hydrostratigraphic unit in reservoir models.

The architecture of preferential flow zones, such as facies St within the FC, may also have thermal and hydraulic implications. Enhanced connectivity between these high-permeability pathways may increase thermal dispersivity (Aghaei et al., 2024) and reduce thermal breakthrough times in geothermal doublets (e.g. Sommer et al., 2013; Di Dato et al., 2022; Molnár et al., 2025).

## 4.5. Broader hydrogeological implications

Beyond implications for deep geothermal brines, unrecognised heterogeneities such as those tested here may also affect the performance of groundwater contaminant remediation (Zheng and Gorelick, 2003) in shallow aquifers west of the Vosges Mountains (BRGM, 1972; 2002; Vaute et al., 2007; 2013) (Fig. 1B). This underscores the importance of exchanging insights obtained in shallow and deep aquifers where similar heterogeneities architectures occur.

## 4.6. Methodological significance

Accurately assessing the subsurface architecture used in this study is rarely feasible, and even the most sophisticated models must simplify natural systems (Hunt and Welter, 2010). Typically, models incorporating fluvial sedimentary architecture consider the deposits of channels, floodplains, and levees as the main building blocks (e.g. Fogg et al., 1998; Maxwell et al., 2008; Zhang et al., 2013; Bianchi and Pedretti, 2017). However, in formations such as the LGV, where floodplain or levee deposits are absent and most deposits formed within channels (Bourquin et al., 2009; Bofill et al., 2025), this differentiation is not applicable.

The approach adopted here, which integrates flow simulations with quantitative sedimentological characterisation (Bofill et al., 2024) and petrophysical analysis (Bofill et al., 2025), proved effective for identifying “unknown-unknowns” (sensu Hunt and Welter, 2010). These are uncertainties that cannot be anticipated, for instance, within purely granulometric conceptual models. As Bredehoeft (2005) noted, modellers cannot correct for conceptual errors they do not recognise. For instance, after identifying the HSS intervals in outcrops and their lower permeability compared to the other facies in the LGV, it was easier to recognise them in core samples and to predict that they were not isolated low-permeability lenses, but stratigraphic intervals that can strongly influence the fluid circulation.

To adequately represent the heterogeneities tested here, it may be appropriate for modellers working on the matrix of the fluvio-aeolian deposits of the LGV to consider using a decimetre-scale cell (~30 cm) in the vertical (z) direction, or to conduct further analyses on the effectiveness of flow-based upscaling. When represented as a homogeneous block in larger-scale models, the model should account for the anisotropy yielded by the WWL intervals, as discussed in 4.3.
## 4.7. Uncertainties

Groundwater modelling is inherently uncertain due to the complexity and heterogeneity of subsurface environments, and this study is no exception. The main sources of uncertainties are outlined below.

### 4.7.1. Unaccounted small-scale heterogeneities

Fine-scale features, such as grain size variations within cross-bedding foresets, were not explicitly represented. While these features can influence fluid flow within individual sets of cross-bedded sandstone (e.g. Ringrose et al., 1993), the resulting *K* variations are minor compared with contrasts between facies. They are therefore unlikely to affect the late-time tail of the breakthrough curve (Zhang et al., 2013) and can, in general, be treated as part of a stochastic random field (Dassargues, 2019).

### 4.7.2. Flow direction assumptions

A constant hydraulic head at the inlet and outlet created a horizontal flow parallel to the palaeo-flow direction (dip). Whether similar solute transport patterns would be observed in the strike direction remains uncertain. However, given the complex erosional and depositional dynamics typical of braided fluvial systems, which result in amalgamated lenticular and tabular deposits (Collinson, 1996), and the lack of preserved floodplain deposits in the LGV, major strike-direction differences are not expected.

### 4.7.3. Reduced *K* variance in stochastic simulations

Monte Carlo simulations produced a lower *K* variance than the original dataset (Table 3), potentially underestimating advective-induced longitudinal dispersivity (Sommer et al., 2013; Di Dato et al., 2022; Ma et al., 2025). Nonetheless, longitudinal macrodispersivity is mainly controlled by permeability contrasts between facies (Zhang et al., 2013; Dong et al., 2018) and the scale and architectural organisation of the heterogeneity (Soltanian et al., 2017; Ma et al., 2025; Molnár et al., 2025), which were preserved.

### 4.7.4. Uniform porosity and only advective-induced longitudinal dispersion accounted

A single porosity value was applied throughout the model, and processes that affect heat and solute transport were omitted, such as: pore-scale longitudinal dispersion, transverse dispersion, molecular diffusion, sorption, and thermal-property heterogeneity. Incorporating facies-specific porosity and transport parameters would enhance model realism and improve representation of geothermal brine behaviour.

### 4.7.5. Near-surface conditions applied to deep systems

Flow simulations were conducted under shallow freshwater conditions, whereas deep geothermal brines are less viscous, slightly denser, and subjected to higher pressures. Natural hydraulic gradients in deep systems are also generally smaller, which may influence the balance between diffusion and convection (Possemiers et al., 2015). However, small-scale heterogeneities are expected to affect outcomes regardless of gradient magnitude, and the study focuses on relative contrasts between sedimentary features. Thus, despite physical differences between shallow and deep systems, the results remain relevant for transport-related uncertainties.

### 4.7.6. Two-dimensional model representation of a three-dimensional system

A 2D vertical model was used to represent the heterogeneities observed in the outcrop realistically; extending the section into 3D would require geometric information that was unavailable and would introduce further uncertainty. Although 2D models can yield reliable predictions (e.g., Prasad and Simmons, 2005; Xie et al., 2012; Ramasomanana et al., 2013; Knorr et al., 2016), they neglect cross-sectional flow, and differences from 3D models can be substantial (e.g., Burnett

and Frind, 1987; Kerrou and Renard, 2010; Pau et al., 2010; Knight et al., 2021; Gholami et al., 2025). Nevertheless, these differences can also be comparable to those arising from alternative representations of geological heterogeneity (e.g. Pau et al., 2010). Thus, the approach reflects a trade-off between achieving a realistic 2D representation of heterogeneity and incorporating additional uncertainty in the third dimension.

## 5. Conclusion

A hydrogeological conceptual model based on sedimentological and petrophysical outcrop data was used to simulate groundwater flow and transport within a realistic representation of aquifer heterogeneity. In a 2D vertical model, decimetre- to metre-scale fluvio-aeolian sandstone facies, often overlooked in groundwater flow simulations, were discretised.

In a reference model, hydraulic conductivity (*K*) values for each facies were assigned and distributed stochastically according to mapped geometries. Particle tracking in FEFLOW® was used to evaluate how the distinct facies influence advective transport across the model domain. Five simplified scenarios were also tested: (1) Facies-scale with *K* assigned deterministically; (2) Facies-association-scale with *K* assigned stochastically; (3) Facies-association-scale with deterministic *K*; (4) Structure-free with stochastic *K*; and (5) Homogeneous *K* field.

Results show that the high-permeability, well-connected fluvial trough cross-bedded sandstone (St) forms preferential flow pathways that control early particle arrivals, whereas the thinner, less permeable, and laterally extensive hybrid sand sheet facies (HSS) act as a retention zone producing the long-tailed late arrivals. When facies within the fluvial channel facies association (FC), St, sigmoidal cross-bedded sandstone (Ss), and low-angle cross-bedded sandstone (Sl), were homogenised, preferential pathways were not represented, altering early-arrival predictions. Stochastic *K* distribution is recommended due to the high variability within facies, which affects early arrival times in the model outlet.

Long-tailed breakthrough curves produced by HSS indicate non-Fickian, pre-asymptotic transport at the model scale, and it is uncertain whether asymptotic behaviour would emerge over longer distances.

Two decimetre-scale, low-permeability, and laterally extensive wind- and water-laid facies association (WWL) deposits yield an upscaled anisotropy of ~0.3 for the model domain (100 m length × 20.6 m height). Homogenising the facies from the FC (St, Ss and Sl) did not change this, indicating that anisotropy is controlled by contrast between WWL and FC rather than by heterogeneity within the FC.

Overall, heterogeneities within the fluvio-aeolian deposits of the Lower Grès Vosgien formation (LGV), often treated as a single hydrostratigraphic unit in the Upper Rhine Graben, strongly influence predictions of solute migration (e.g. lithium in geothermal brines or contaminants in drinking-water aquifers) and thermal distribution within the matrix reservoir. Capturing these features leads to more reliable models for managing lithium-rich geothermal brines and protecting drinking water resources.

## CRediT authorship contribution statement

**Lucas Bofill:** Writing – review & editing, Writing – original draft, Visualization, Validation, Project administration, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Gerhard Schäfer:** Writing – review & editing, Validation, Supervision, Software, Methodology, Funding acquisition, Conceptualization. **Guilherme Bozetti:** Writing – review & editing, Supervision, Resources, Funding acquisition, Conceptualization. **Jean-François Ghienne:** Writing – review & editing, Conceptualization. **Mathieu Schuster:** Writing – review & editing, Conceptualization.
Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This study is part of Lucas M. Bofill's PhD project, funded by the
École Doctorale des Sciences de la Terre et de l'Environnement (ED413).
This research was also supported by: (1) the Interdisciplinary Thematic
Institute (ITI) GeoT, part of the 2021-2028 ITI program of the University
of Strasbourg, the French National Research Centre (CNRS), and the
French National Institute of Health and Medical Research (Inserm). The
ITI GeoT was supported by the SFRI-STRAT'US project (ANR-20-SFRI-
0012) and by IdEx Unistra (ANR-10-IDEX-0002) under the Investments
for the Future (PIA) framework; and (2) The French-Azerbaijani Uni-
versity (UFAZ). We acknowledge GEIE EMC and ES-Géothermie, more
specifically Albert Genter, and Carole Glaas, for providing the structure
and access to the core EPS1. We also thank Dr Luca Colombera and the
anonymous reviewer for their time and effort invested, which greatly
contributed to this manuscript. Last but not least, special thanks to Laura
Menezes for her personal support and for enriching the content and
discussions of this study.

Supplementary materials

Supplementary material associated with this article can be found, in
the online version, at doi:10.1016/j.geothermics.2025.103581.

Data availability

Data will be made available on request.

References

Åberg, S.C., Aberg, A.K., Korkka-Niemi, K., 2021. Three-dimensional hydrostratigraphy and groundwater flow models in complex quaternary deposits and weathered/ fractured bedrock: evaluating increasing model complexity. *Hydrogeol. J.* 29, 1043–1074. https://doi.org/10.1007/s10040-020-02299-4.

Aghaei, H., Colombera, L., Yan, N., Mountney, N.P., Brenagozzi, G., Di Giulio, A., 2024. Impact of fluvial meander-belt sedimentary heterogeneity on the efficiency of low-enthalpy geothermal boullets: heat-transport simulations of forward stratigraphic models. *Geoenergy* 2, 16. https://doi.org/10.1144/geoenergy2024-024.

Aichholzer, C., 2019. Le log complet de la stratigraphie de la zone rhénane ainsi que les modalités stratigraphiques, sédimentaires et structurales de la transition « socle-couverture. Aplicacion a la géothermie profonde. PhD thesis. Université de Strasbourg. École doctorale des Sciences de la Terre et de l'Environnement (ED413), Intitut de Physique du Globe de Strasbourg (UMR7516). 546 p (In French).

Aigner, T., Aspiron, U., Hornung, J., Junghans, W.D., Kostrewa, R., 1996. Integrated outcrop analogue studies for triassic alluvial reservoirs: examples from southern Germany. J. Petrol. Geol. 19, 393–406.

Alexander, 1993. A discussion on the use of analogues for reservoir geology, 1993. In: Ashton, M. (Ed.), Advances in Rervoir Geology, Advances in Rervoir Geology, 69. Geological Society Special Publications, pp. 175-194.

Anderson, M.P., Aiken, J.S., Webb, E.K., Mickelson, D.M., 1999. Sedimentology and hydrogeology of two braided stream deposits. Sediment. Geol. 129, 187-199. https://doi.org/10.1016/S0037-0738(99)00015-9.

Aquilina, L., Pauwels, Genter A., Fouillac, C., 1997. Water-rock interaction processes in the triassic sandstone and the granitic basement of the Rhine Graben: geochemical investigation of a geothermal reservoir. *Geochim. Cosmochim. Acta* 61 (20), 4281–4295. N.

Aris, R., 1958. On the dispersion of linear kinematic waves. Proc. R. Soc. A and B 245 (1241), 268. https://doi.org/10.1098/rspa.1958.0082.

Behrens, H., Ghergut, J., Sauter, M., Wagner, B., Wiegand, B., 2022. Solute Co-
production from small geothermal reservoirs – How little is too little. In:
*Proceedings, 47th Workshop on Geothermal Reservoir Engineering*. Standford
University, California, p. 11. February 7-9, 2022. SGP-TR-223.

Bianchi, M., Pedretti, D., 2017. Geological entropy and solute transport in heterogeneous porous media. Water Resour. Res. 53, 4691–4708. https://doi.org/10.1002/2016WR020195.

Bianchi, M., Zheng, C., 2016. A lithofacies approach for modelling non-Fickian solute transport in a heterogeneous alluvial aquifer. Water Resour. Res. 52, 552–565. https://doi.org/10.1002/2015WR018186.

Bianchi, M., Zheng, C., Wilson, C., Tick, G.R., Liu, G., Gorileck, S.M., 2011. Spatial connectivity in a highly heterogeneous aquifer: from cores to preferential flow paths. Water Resour. Res. 47, 18. https://doi.org/10.1029/2009WR008966 p.

Bjørlykke, K., 2014. Relationships between depositional environments, burial history, and rock properties. Some principal aspects of diagenetic process in sedimentary basins. Sediment. Geol. 301, 1–14. https://doi.org/10.1016/j.sedgeo.2013.12.002.

Bofill, L.M., Bozetti, G., Schäfer, G., Ghienne, J.F., Schuster, M., Scherer, C., Souza, E., 2024. Quantitative facies analysis of a fluvio-aeolian system: lower Triassic Buntsandstein Group, eastern France. Sediment. Geol. 465, 22. https://doi.org/10.1016/j.sedgeo.2024.106634, 106634.

Bofill, L.M., Bozetti, G., Schäfer, G., Ghienne, J.F., Schuster, M., Heap, J.M., Knobelock, G., Scherer, C., Armelenti, G., Souza, E., 2025. Sedimentary control on permeability heterogeneity: the middle Buntsandstein continental sandstones (Lower Triassic, eastern France). Mar. Pet. Geol. 173, 22. https://doi.org/10.1016/j.marpetgeo.2024.107261.

Bouly, S., Demassieux, L., 1996. Structure, caractéristiques et fonctionnement du réservoir aquifère des Grès inférieurs du Trias en Lorraine. La Houille blanche (76). In French.

Bourquin, S., Guillocheau, F., Peron, S., 2009. Braided rivers within an arid alluvial plain (example from the Lower Triassic, western German Basin): recognition criteria and expression of stratigraphic cycles. Sedimentology 56, 2235–2264. https://doi.org/10.1111/j.1365-3091.2009.01078.x.

Bredehoeft, J., 2005. The conceptualization model problem-surprise. *Hydrogeol. J.* 13, 37–46. https://doi.org/10.1007/s10040-004-0430-5.

BRGM, 1972. Etude Hydrogeologique de La Napperre Aquifère des Gres Infratriassiques Dans Le Nord-Est de La France. Service de La Carte Géologique d'Alsace Et De Lorraine. Université de Strasbourg.

BRGM, 2002. Grès du Trias inférieur (GTI) - Code 210. Atlas hydrogéologique du bassin Rhin-Meuse (Octobre 2002). https://sigesrm.brgm.fr/La-BD-RHF-R-V1-dans-le-bassin-in-Rhin-Meuse accessed on 10th of May 2025.

Burkardt, J., 2023. The Truncated Normal Distribution. Department of Scientific Computing, Florida State University. Accessed through. https://people.sc.fsu.edu/~jburkardt/presentations/truncatednormal.pdf.

Burnett, R.D., Frind, E.O., 1987. Simulation of contaminant transport in three dimensions: 2. Dimensionality effects. Water Resour. Res. 23, 695–705. https://doi.org/10.1029/WR0231004p00695.

Busch, B., Spitzner, A.D., Adelmann, D., Hilgers, C., 2022. The significance of outcrop analog data for reservoir quality assessment: a comparative case study of Lower Triassic Buntsandstein sandstones in the Upper Rhine Graben. Mar. Pet. Geol. 141, 105701. https://doi.org/10.1016/j.marpetgeo.2022.105701.

Busch, B., Böke, L.L., Schmid, C., Warnecke, M., Hilgers, C., 2025. Reservoir quality controls in the Lower Buntsandstein and their relation to detrital grain size, URG, SW-Germany. Geothermics 132, 18. https://doi.org/10.1016/j.geothermics.2025.103453.

Casillas-Trasvina, A., Rogiers, B., Beerten, K., Wouters, L., Walraevens, K., 2022. Characterizing groundwater heat transport in a complex lowland aquifer using paleo-temperature recognition, satelitte data, temperature-depth profiles, and numerical models. Hydrol. Earth Syst. Sci. 26, 5577–5604. https://doi.org/10.5194/hess-26-5577-2022.

Castro, A.P., Guedes, C.C.F., Vesely, F.F., Silveira, D.M., Neves, L.F., 2023. Provenance, diagenesis and reservoir quality of sandstones facies of the Maracangalha Formation, Reconcavo Basin - Northeastern Brazil. Mar. Pet. Geol. 158, 26. https://doi.org/10.1016/j.marpetgeo.2023.106548.

Chandler, M.A., Kocurek, G., Goggin, D.J., Lake, L.W., 1989. Effects of stratigraphic heterogeneity on permeability in Eolian Sandstone Sequence, Page Sandstone, Northern Arizona. Am. Assoc. Petrol. Geol. Bull. 73 (5), 658–668.

Collinson, J.D., 1969. The sedimentology of the Grindslow shales and the Kinderscout Grit: a deltaic complex in the Namurian of northern England. J. Sediment. Res. 39, 194–221.

Collinson, J.D., 1996. Alluvial sediments, 1996. In: Reading, H.G. (Ed.), Sedimentary Environments: Processes, Facies and Stratigraphy, 3rd Edition. University of Oxford.

Comolli, A., Hakoun, V., Dentz, M., 2019. Mechanisms, upscaling, and prediction of anomalous dispersion in heterogeneous porous Media. Water Resour. Res. 55, 8197-8222. https://doi.org/10.1029/2019WR024919.

Cooijmans, R.A., Willems, C.J.L., Nick, H.M., Bruhn, D.F., 2016. The influence of facies heterogeneity on the doublet performance in low-enthalpy geothermal sedimentary reservoirs. Geothermics 64, 209–219. https://doi.org/10.1016/j.geothermics.2016.06.004.

Dachroft, W., 1985. Fluvial sedimentary styles and associated depositional environments in the Buntsandstein Wes of River Rhine in Saar Area and Pfalz (F.R. Germany) and Vosges (France). In: mader, D., aspects of fluvial sedimentation in the Lower Triassic Buntsandstein of Europe. Lect. Note Earth Sci. 4, 197–248. Berlin/Heidelberg /New York / Tokio (Springer).

Dassargues, A., 2019. Hydrogeology: Groundwater science and Engineering. CRC Press, Taylor & Francis Group. ISBN: 13: 978-1-4987-4400-3.

Davis, J.M., Lohmann, R.C., Phillips, F.M., Wilson, J.L., Love, D.W., 1993. Architecture of the Sierra Ladrones Formation, Central New Mexico: Depositional controls On the Permeability Correlation Structure, 105. Geological Society of America Bulletin, pp. 998-1007.

de Lange, WJ., 2020. Adveptive transport phenomena to better understand dispersion in field and modeling practice. Ground Water 58 (1), 46–55. https://doi.org/10.1111/gwat.12883.

Di Dato, M., D'Angelo, C., Casasso, A., Zarlenga, A., 2022. The impact of porous medium heterogeneity on the thermal feedback of open-loop shallow geothermal systems. J. Hydrol. 64, 14. https://doi.org/10.1016/j.jhydrol.2021.127205.
Diersch, H.J.G., 2014. FEFLOW: Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media. Springer, Heidelberg, Germany. https://doi.org/10.1007/978-3-642-38739-5-1.

Dong, S., Dai, Z., Li, J., Zhou, W., 2018. The scale dependence of dispersivity in multi-facies heterogeneous formations. Carbon. Evap. 33, 161–165. https://doi.org/10.1007/s13146-018-0421-6.

Dugamain, E.J.M., Richard, A., Cathelineau, M., Boiron, M.C., Despinois, F., Brisset, A., 2021. Groundwater in sedimentary basin as potential lithium resource: a global prospective study. Scientific Reports, Nature Portfolio. https://doi.org/10.1038/s41598-021-99912-7.

Duringer, P., Aichholzer, C., Orciani, S., Genter, A., 2019. The complete lithostratigraphic section of the geothermal wells in Rittershoffen (Upper Rhine Graben, eastern France): a key for future geothermal wells. BSGF Earth Sci. Bull. 2019, 190008. Available at www.bsgf.fr.

Farquharson, J.I., Heap, M.J., Lavallée, Y., Varley, N.R., Baud, P., 2016. Evidence for the development of permeability anisotropy in lava domes and volcanic conduits. J. Volcanol. Geotherm. Res. 323, 163–185.

Ferguson, G., 2007. Heterogeneity and thermal modeling of ground water. Ground Water 45, 485–490. https://doi.org/10.1111/j.1745-6584.2007.00323.x.

Fischer, C., Gaupp, R., Dimke, M., Sill, O., 2007. A 3D high resolution model of bounding surfaces in aeolian-fluvial deposits: an outcrop analogue study from the Permian Rotliefend, Northern Germany. J. Petrol. Geol. 30, 257–274.

Fogg, G.E., Noyes, C.D., Carle, S.F., 1998. Geologically based model of heterogeneous hydraulic conductivity in an alluvial setting. Hydrogeol. J. 6, 131–143. https://doi.org/10.1007/s100400050139.

Fries, D., Leboul, S., Maurer, V., Martin, C., Baujard, C., Ravier, G., Bobuais, R., Amari, S., 2022. Lithium Extraction Through Pilot Scale Test Under Real Geothermal Conditions of the Upper Rhine Graben. European Geothermal Congress, Berlin, Germany, p. 8, 202217–21 October 2022.

Gelhar, L.W., Gutjahr, A.L., Naff, R.L., 1979. Stochastic analysis of macrodispersion in stratified aquifer. Water Resour. Res. 15, 11.

Gelhar, L.W., Welty, C., Rehfeldt, K.R., 1992. A critical review of data on field-scale dispersion in aquifers. Water Resour. Res. 28, 1955–1974.

Ghergut, J., Wiegand, B., Behrens, H., Sauter, M., 2023. Model-independent and Model-dependent aspects of “Geothermal Solute” Co-production forecast for hydrothermal vs. In: Petrothermal Reservoirs. Proceedings, 48th Workshop on Geothermal Reservoir Engineering. Stanford University, California, p. 9. February 6-8. SGP-TR-224.

Ghienne, J.F., Deynoux, M., Jeannetter, D., 1994. Structures sédimentaires, discontinuités et caractères pétrophysiques des grès du Buntsandstein des Vosges du Nord. T rias de l'Est de la France. C.R. Acad. Sci. Paris, t. 319, série II, pp. 1351–1358 (in French, abstract in English).

Gholami, A., Jazayeri, A., Werner, A.D., 2025. Cross-sectional models of groundwater flow: review and correlation for transverse flow. Groundwater 63, 752–763. https://doi.org/10.1111/gwat.70017.

Glaas, C., Salmon, M., Baron, F., Bozetti, G., Knobelock, G., Rumbach, G., Patrier, P., Baczkowski, C., Genter, A., 2025. Are Triassic rocks of the Upper Rhine Graben a geothermal reservoir? Insights from Rittershoffen wells with Geophysical Logs and SWIR Spectroscopy. Presented at the European Geothermal Congress, Zürich, Switzerland.

Goupil, M., Heap, M., Baud, P., 2022. Permeability anisotropy in sandstones from the Soultz-sous-Forêts geothermal reservoir (France): implications for large-scale fluid flow modelling. Geotherm. Energ. 10 (32), 23. https://doi.org/10.1186/s40517-022-00243-1.

Guillocheau, F., Robin, C., Allemand, P., Bourquin, S., Brault, N., Dromart, G., Friedenberg, R., Garcia, J.P., Gaulier, J.M., Gaumet, F., Grosdoy, B., Hanot, F., Le Strat, P., Mettraux, M., Nalpas, T., Prijac, C., Rigollet, C., Serrano, O., Grandjean, G., 2000. Meso-cenozoic geodynamic evolution of the Paris Basin: 3D stratigraphic constraints. Geodin. Acta. 13, 189–246.

Haffen, S., Geraud, Y., Diraison, M., Dezayes, C., 2013. Determination of fluid-flow zones in a geothermal sandstone reservoir using thermal conductivity and temperature logs. Geothermics 46, 32–41. https://doi.org/10.1016/j.geothermics.2012.11.001.

Hakoun, V., Comolli, A., Dents, M., 2019. Upscaling and prediction of lagrangian velocity dynamics in heterogeneous porous media. Water Resour. Res. 55, 3976–3996. https://doi.org/10.1029/2018WR023810.

Heap, M.J., Kushnir, A.R., Gilg, H.A., Wadsworth, F.B., Reuschlé, T., Baud, P., 2017. Microstructural and petrophysical properties of the Permo-triassic sandstones (Buntsandstein) from the Soultz-sous-Forêts geothermal site (France). Geotherm. Energ. 5, 1–37.

Heidsiek, M., Butscher, C., Blum, P., Fischer, C., 2020. Small-scale diagenetic facies heterogeneity controls porosity and permeability pattern in reservoir sandstones. Environ. Earth Sci. 79 (425), 14. https://doi.org/10.1007/s12665-020-09168-z.

Herivaux, C., Marechal, J., 2019. Prise en compte des services dépendants des aquifères dans les démarches d'évaluation des services écosystémiques. Rapport final, BRGM/RP-68929-FR. 107 p (In French).

Herr, M., Schäfer, G., Spitz, K., 1989. Experimental studies of mass transport in porous media with local heterogeneities. J. Contam. Hydrol. 4, 127–137.

Høimyr, O., Kleppe, A., Nystuen, J.P., 1993. Effects of heterogeneities in a braided stream channel sandbody on the simulation of oil recovery: a case study from the Lower Jurassic Statfjord Formation, Snorre Field, North Sea. From: Ashton, M. (ed.), 1993. Advances in Reservoir Geology. Geol. Soc. Publ. 69, 105–134.

Howell, J.A., Martinius, A.W., Good, T.R., 2014. The application of outcrop analogues in geological modelling: a review, present status and future outlook. Geol. Soc. Lond. Spec. Pub. 387 (1), 1–25. https://doi.org/10.1144/sp387.12.

Hunt, R.J., Welter, D.E., 2010. Taking account of “unknown unknowns”. Groundwater 48, 477. https://doi.org/10.1111/j.1745-6584.2010.00681.x.

Huysmans, M., Peeters, L., Moermans, G., Dassargues, A., 2008. Relating small-scale sedimentary structures and permeability in a cross-bedded aquifer. J. Hydrol. 361, 41–51. https://doi.org/10.1016/j.jhydrol.2008.07.047.

Iversen, B., Moldrup, V., Schjonning, P., Jacobsen, O.H., 2003. Field application of a portable air permeameter to characterize spatial variability in air and water permeability. Vadose Zone J. 2, 618–626. https://doi.org/10.2136/vzj2003.6180.

Jung, R., 2013. EGS – Goodbye or Back to the future, Effective and Sustainable Hydraulic Fracturing. Intech, pp. 95–121. https://doi.org/10.5772/56458Jussel. Chapter 5P., 1989. Stochastic description of typical inhomogeneities of hydraulic conductivity in fluvial gravel deposits. In: Kobus, H.E., Kinzelbach, W., 1989. Proceedings of the International Symposium on Contaminant Transport in Groundwater, Stuttgart, 4-6 April 1989.

Kerrou, J., Renard, P., 2010. A numerical analysis of dimensionality and heterogeneity effects on advecteddispersive seawater intrusion processes. Hydrogeology 18, 55–72. https://doi.org/10.1007/s10040-009-0533-0.

Klise, A.K., Weissmann, G.S., McKenna, S.A., Nichols, E.M., Frechette, J.D., Wawrzyniec, T.F., Tidwell, V.C., 2009. Exploring solute transport and streamline connectivity using lidar-based outcrop images and geostatistical representations of heterogeneity. Water Resour. Res. 45, 11. https://doi.org/10.1029/2008WR007500.

Knight, A.C., Irvine, D.J., Werner, A.D., 2021. Alongshore freshwater circulation in offshore aquifers. J. Hydrol. 593, 125915. https://doi.org/10.1016/j.jhydrol.2020.125915.

Knorr, B., Xie, Y., Stumpp, C., Maloszewski, P., Simmons, C.T., 2016. Representativeness of 2D models to simulate 3D unstable variable density flow in porous media. J. Hydrol. 542, 541–551. https://doi.org/10.1016/j.jhydrol.2016.09.026.

Kölbel, L., Kölbel, T., Sauter, M., Schäfer, T., Siefert, D., 2020. Identification of fracture zones in geothermal reservoirs in sedimentary basins: a radionuclide-based approach. Geothermics 85, 10.

Kölbel, L., Kölbel, T., Herrmann, L., Kaymakci, E., Ghergut, I., Poirrel, A., Schneider, J., 2023. Lithium extraction from geothermal brines in the Upper Rhine Graben: a case study of potential and current state of the art. Hydrometallurgy 221, 106131. https://doi.org/10.1016/j.hydromet.2023.106131.

Koltermann, C.E., Gorelick, S.M., 1996. Heterogeneity in sedimentary deposits: a review of structure-imitating, process-imitating, and descriptive approaches. Water Resour. Res. 32, 2617–2658.

Kreft, A., Zuber, A., 1978. On the physical meaning of the dispersion equation and its solutions for different initial and boundary conditions. Chem. Eng. Sci. 33, 1471–1480.

Loll, P., Moldrup, P., Schjonning, P., Riley, H., 1999. Predicting saturated hydraulic conductivity from air permeability: application in stochastic water infiltration modeling. Water Resour. Res. 35, 2387–2400.

Ma, Z., Qi, L., Dai, Z., Yang, Z., Ma, Y., Wang, D., Carroll, K.C., Soltanian, M.R., 2025. Impact of multiscale heterogeneous sediments and boundary conditions on dispersivity spatial variations. Water Resour. Res. 61. https://doi.org/10.1029/2024WR039151 e2024WR039151.

Makeen, Y.M., Shan, X., Lawal, M., Ayinla, H.A., Su, S., Yelwa, N.A., Liang, Y., Ayuk, N. E., Du, X., 2021. Reservoir quality and its controlling diagenetic factors in the Bentiu Formation, Northeastern Muglad Basin, Sudan. Sci. Rep. 11, 18442. https://doi.org/10.1038/s41598-021-97994-x.

Maxwell, R.M., Carle, S.F., Tompson, A.F.B., 2008. Contamination, risk, and heterogeneity: on the effectiveness of aquifer remediation. Environ. Geol. 54, 1771–1786. https://doi.org/10.1007/s00254-007-0955-8.

Miall, A.D., 1988. Architectural elements and bounding surfaces in fluvial deposits: anatomy of the Kayenta formation (Lower Jurassic), Southwest Colorado. Sediment. Geol. 55, 233–262.

Molnár, B., Lenkey, L., Pedretti, D., Mádl-Szönyi, J., Galsa, A., 2025. Elucidating the relationship between aquifer heterogeneity and thermal dispersivity using stochastic heat transport simulations. Hydrogeol. J. 33, 493–513. https://doi.org/10.1007/s10040-025-02875-6.

Newell, A.J., Pourmalek, A., Butcher, A.S., Shariatipour, S.M., 2019. The importance of lithofacies control on fluid migration in heterogeneous aeolian formations for geological CO2 storage: lessons from observational evidence and modelling of bleached palaeoreservoirs at Salt Wash Graben. 91, 17 p https://doi.org/10.1016/j.ijggc.2019.102841.

O'Sullivan, J., Araya, N., Popineau, J., Renaud, T., Riffault, J., 2023. A natural State and production forecast model of the Santon Sea geothermal field for lithium extraction. GRC Trans. 47, 20.

Pau, G.S.H., Bell, J.B., Pruess, K., Almgren, A.S., Lijexski, M.J., Zhang, K., 2010. High-resolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers. Adv. Water Resour. 33, 443–455. https://doi.org/10.1016/j.adwwatres.2010.01.009.

Posamantier, H.W., Jervey, M.T., Vail, P.R., et al., 1988. Eustatic controls on clastic deposition I – conceptual framework. eds. In: Wilgus, C.K., et al. (Eds.), Sea-level changes: an Integrated approach: Society of Economic Paleontologists and Mineralogists. Special Publications, pp. 109–124, 42.

Possemiers, M., Huysmans, M., Peeters, L., Batelaan, O., Dassargues, A., 2012. Relationship between sedimentary features and permeability at different scales in the Brussels Sands. Geologica Belgica 15/3, 156–164.

Possemiers, M., Huysmans, M., Batelaan, O., 2015. Application of multiple-point geostatistics to simulate the effect of small-scale heterogeneity on the efficiency of aquifer thermal energy storage. Hydrogeol. J. 23, 971–991. https://doi.org/10.1007/s10040-015-1244-3.

Prasad, A., Simmons, C.T., 2005. Using quantitative indicators to evaluate results from variable-density groundwater flow models. Hydrogeol. J. 13, 905–914. https://doi.org/10.1007/s10040-004-0338-0.
Ramanathan, R., Ritzi Jr., R.W., Allen-King, R.M., 2010. Linking hierarchical stratal architecture to plume spreading in a Lagrangian-based transport model: 2. Evaluation using new data from the Borden Site. Water Resour. Res. 46, 1. https://doi.org/10.1029/2007WR006282.

Ramasomanana, F., Younes, A., Ackerer, P., 2013. Estimation of macrodispersion in 2-D highly heterogeneous porous media using the Eulerian-Langrangian localized adjoint method. Water Resour. Res. 49, 43–53. https://doi.org/10.1029/2012WR012228.

Ren, W., Ershadnia, R., Wallace, D.W., LaBolle, E.M., Dai, Z., Barros, F.P.L., Soltanian, M. R., 2022. Evaluating the effects of multiscale heterogeneous sediments on solute transport mixing and effective dispersion. Water Resour. Res. 58, 9. https://doi.org/10.1029/2021WR031886.

Ringrose, P.S., Sorbie, K.S., Corbett, P.W.M., Jensen, J.L., 1993. Immiscible flow behaviour in laminated and cross-bedded sandstones. J. Petrol. Sci. Eng. 9, 103–124. https://doi.org/10.1016/0920-4105(93)00701-L.

Rogiers, B., Beerten, K., Smeekens, T., Mallants, D., Gedeon, M., Huymans, M., Batelaan, O., Dassargues, A., 2013. The usefulness of outcrop-analogue air-permeameter measurements for analysing aquifer heterogeneity: testing outcrop hydrogeological parameters with independent borehole data. Hydrol. Earth Syst. Sci. 17, 5155–5166. https://doi.org/10.5194/hess-17-5155-2013.

Sanjuan, B., Gourcerol, B., Millot, R., Rettenmaier, D., Jeandel, E., Rombaut, A., 2022. Lithium-rich brines in Europe: an up-date about geochemical characteristics and implications for potential Li resources. Geothermics 101, 18.

Sauty, J.P., Gringarten, A.C., Frabris, H., Thiery, D., Menjoz, A., Delnel, P.A., 1982. Sensible energy storage in aquifers 2. Field experiments and comparison with theoretical results. Water Resour. Res. 18, 253–265. https://doi.org/10.1029/WR018i002p00253.

Schäfer, G., 1991. Einfluß von schichtenstrukturen und lokalen einlagerungen auf die längsdispersion in Poregrundwasserleitern. M itteilungen des Instituts für Wasserbau, Heft 75, Universität Stuttgart, 154p, ISBN 3-921694-75-2.

Scorgie, J.C., Worden, R.H., Utley, J.E.P., Roche, I.P., 2021. Reservoir quality and diagenesis of triassic sandstones and siltstones from arid fluvial and playa margin environments: a study of one of the UK's earliest producing oilfields. Mar. Pet. Geol. 131, 28. https://doi.org/10.1016/j.marpetgeo.2021.105154.

Snider, S., Griffiths, C.M., Owen, A., Hartley, A.J., Howell, J.A., 2021. Stratigraphic forward modelling of distributive fluvial systems based on the Huesca system, Ebro Basin, northern Spain. Basin Res. 33, 3137–3158. https://doi.org/10.1111/bre.12597.

Soltanian, M.R., Sun, A., Dai, Z., 2017. Reactive transport in the complex heterogeneous alluvial aquifer of Fortymile Wash, Nevada, Chemosphere. doi: 10.1016/j.chemosphere.2017.03.136.

Sommer, W., Valstar, J., Gaans, P.V., Grotenhuis, T., Rijnaarts, H., 2013. The impact of aquifer heterogeneity on the performance of aquifer thermal energy storage. Water Resour. Res. 49, 8128–8138. https://doi.org/10.1002/2013WR013677.

Talon, L., Ollivier-Triquet, E., Dentz, M., Bauer, D., 2023. Transient dispersion regimes in heterogeneous porous media: on the impact of spatial heterogeneity in permeability and exchange kinetics in mobile-immobile transport. Adv. Water Resour. 174, 104425. https://doi.org/10.1016/j.advwatres.2023.104425.

Theel, M., Huggenberger, P., Zosseder, K., 2020. Assessment of the heterogeneity of hydraulic properties in gravelly outwash plains: a regionally scaled sedimentological analysis in the Munich gravel plain, Germany. Hydrogeol. J. 28, 2657–2674. https://doi.org/10.1007/s10040-020-02205-y.

Théobald, N., Perriaux, J., Langenfeld, F., Both, J., 1967. Carte géologique 1/50000 bitche- Walschbronn XXXVIII3. (In French).

Thompson, J., Parkes, D., Hough, E., Wakefield, O., 2019. Using core and outcrop analogues to predict flow pathways in the subsurface: examples from the triassic sandstones of north Cheshire, UK. Adv. Geosci. 49, 121–127.

Valocchi, A.J., 1985. Validity of the local equilibrium assumption for modeling sorbing solute transport through homogeneous soils. Water Resour. Res. 21 (6), 808–820.

Vaute, L., Gigleux, S., Ngyen-The, 2007. Eaux souterraines du département des Vosges: caractérisation des principals ressources exploitables et revision du modèle de gestion de la nappe des grès du Trias inférieur. Rapport final, BRGM/RP-55653-FR. 396 p (In French).

Vaute, L., Innocent, C., Fourniguet, 2013. Actualisation du modèle hydrogeologique de la nappes des grès du Trias en Lorraine. Rapport final. BRGM/RP-62405-FR. 60p (In French).

Vidal, J., Genter, A., 2018. Overview of naturally permeable fractured reservoirs in the central and southern Upper Rhine Graben: insights from geothermal wells. Geothermics 74, 53–73. https://doi.org/10.1016/j.geothermics.2018.02.003.

Wadsworth, F.B., Vasseur, J., Heap, M.J., Carbillet, L., Dingwell, D.B., Reuschlé, T., Baud, P., 2023. A universal model for the permeability of sintered materials. Acta Mater. 250, 118859.

Walker, R.G., 1992. Facies, facies models and modern stratigraphic concepts. In: Walker, R.G., James, N.P. (Eds.), Facies Models - in Response to Sea Level Change, 1–14. ISBN 0-0919216-49-8.

Webb, E.K., Anderson, M.P., 1996. Simulation of preferential flow in three-dimensional heterogeneous conductivity fields with realistic internal architecture. Water Resour. Res. 32, 533–545. https://doi.org/10.1029/95WR03399.

Weissmann, G.S., Pickel, A., McNamara, K.C., Frechette, J.D., Kalinovich, I., Allen-King, R.M., Jankovic, I., 2015. Characterization and quantification of aquifer heterogeneity using outcrop analogs at the Canadian Forces Base Borden, Ontario, Canada. Geol. Soc. Am. Bull. 127, 1021–1035. https://doi.org/10.1130/B31193.1.

Woessner, W.W., Poeter, E.P., 2024. Hydrogeologic properties of earth materials and principles of groundwater flow. The Groundwater Project. 7th ed. https://doi.org/10.21083/978-1-7770541-2-0.

Xie, Y., Simmons, C.T., Werner, A.D., Diersch, H.J.G., 2012. Prediction and uncertainty of free convection phenomena in porous media. Water Resour. Res. 48. https://doi.org/10.1029/2011WR011346.

Yan, N., Colombera, L., Cosgrove, G.I.E., Mountney, N.P., 2024. A 3D forward stratigraphic model of aeolian dune evolution for prediction of lithofacies heterogeneity. Comp. Geosci. 187, 15. https://doi.org/10.1016/j.cageo.2024.105594.

Yin, M., Ma, R., Zhang, Y., Lin, J., Guo, Z., Zheng, C., 2023. Competitive control of multiscale aquifer heterogeneity on solute transport in an alluvial aquifer. J. Hydrol. 616, 13. https://doi.org/10.1016/j.jhydrol.2022.128819.

Zhang, Y., Benson, D.A., Baeumer, B., 2007. Predicting the tails of breakthrough curves in regional-scale alluvial systems. Ground Water 45, 473–484. https://doi.org/10.1111/j.1745-6584.2007.00320.x.

Zhang, Y., Green, C.T., Fogg, G.E., 2013. The impact of medium architecture of alluvial settings on non-fickian transport. Adv. Water Resour. 54, 78–99. https://doi.org/10.1016/j.advwatres.2013.01.004.

Zheng, C., Gorelick, S.M., 2003. Analysis of solute transport in flow fields influenced by preferential flowpaths at the decimeter scale. Groundwater 41, 142–155. https://doi.org/10.1111/j.1745-6584.2003.tb02578.x.