Stratigraphic forward modelling of distributive fluvial systems based on the Huesca System, Ebro Basin, northern Spain

To better understand the stratigraphic development of sedimentary systems, it is necessary to link the controls on sedimentary processes to the resulting deposits, which in turn allows predictions of stratigraphic architectures at a range of scales. We use a stratigraphic forward model to link the governing parameters to the distribution of deposits within a distributive fluvial system (DFS). The numerical model has been validated against outcrop observations to establish how the depositional processes needed to form the specific sedimentary system have been reproduced. We chose the previously studied Oligocene to Miocene Huesca DFS in northern Spain to investigate and calibrate the model. Additionally, downstream profiles from modern DFS in northern India, and hydrological measurements from the High Island Creek, Minnesota, USA, were used as input parameters for the model in addition to the outcrop data from the Huesca DFS. The resulting model adequately reproduced the real‐world system. Once validated, the analysis of the modelled DFS led to key findings, which expand our understanding of DFS stratigraphic architecture. Reservoir characteristics in radial DFS are dependent on the angle away from the meridian (straight line from the source through the apex to the distal zone of the DFS). The greater the angle is, the coarser the average grain size in the proximal zone is but the finer the average grain size in the medial and distal zones. Lateral variability of net to gross, sandbody thickness and number, and amalgamation ratio is greatest at the transition between the proximal and medial zone and is still significant in the distal part of the DFS. Stratigraphic forward modelling enhanced our understanding of DFS, which leads to reducing risk associated with exploration, production and storage of fluids in subsurface DFS.


| INTRODUCTION
To better understand the stratigraphic evolution of fluvial sedimentary systems, their controlling parameters need to be linked to their deposits. Numerical forward modelling can be used to make this connection. Hydrodynamic stratigraphic forward models (SFM) require no prior knowledge of the depositional geometry such as channel planform, as physical formulas for fluid flow and sediment erosion, transport, and deposition compute the distribution of sediments over natural temporal and spatial scales. This echoes Best and Fielding (2019) who said that the next step in describing fluvial systems is to link processes to deposits and stratigraphy. They also proposed the use of facies models which should be based on physical process distinctions to achieve this next step. Ultimately, the numerical or facies models will lead the way to a predictive model of sediment body architecture, which is crucial for derisking subsurface exploitation and storage of fluids (e.g., hydrocarbons, ground water, geothermal energy, hydrogen and carbon dioxide) in these sedimentary systems.
Distributive fluvial systems were chosen as the fluvial system to analyse as they constitute a significant component of fluvial systems in modern continental basins  and occur in every climate, tectonic, and basin setting Weissmann et al., 2010). As a result, they may be considered prevalent in the rock record. Commonly, they are characterized by downstream decreases in channel dimensions and grain size, with channel confinement decreasing towards the toe of the systems . The term DFS was coined by Hartley et al. (2010) and defined as 'the deposit of a fluvial system which in planform displays a radial, distributive channel pattern'. Specific DFS character had been described previously in both modern and ancient settings, with these deposits termed 'alluvial fans', 'fluvial fans', 'humid alluvial fans' and 'megafans' (e.g. Beaty, 1963;Cain & Mountney, 2009;Chakraborty et al., 2010;Dade & Friend, 1998;Gumbricht et al., 2004;Heward, 1978;Mack & Leeder, 1999;Mohindra et al., 1992;Nichols & Fisher, 2007;Whipple et al., 1998).
The DFS model was derived from interpreting remote sensing data of 700 modern continental sedimentary basins  and 415 DFS (Davidson et al., 2013;Hartley et al., 2010). Even though these studies provided an understanding of modern DFS in planform, they mostly covered a short time period (Late Pleistocene to present day). The analysis of ancient fluvial systems in outcrop has provided a temporal scale to the DFS model, particularly quantified outcrop studies such as from the Brownstones, western UK (Allen, 1983); the Huesca and Luna DFS, northern Spain (Burnham, 2016;Hirst, 1991;Hirst & Nichols, 1986;Kulikova, 2013;Martin et al., 2021;Nichols, 1987); the Salt Wash DFS, southwestern USA (Owen et al., , 2017; the Blackhawk Formation, southwestern USA (Rittersbacher et al., 2014); the Marília Formation, southeastern Brazil (Dal' Bó et al., 2019;Soares, Basilici, Lorenzoni, et al., 2020;Soares, Basilici, Silva Marinho, et al., 2020) and the Tordillo Formation, central Argentina (Coronel et al., 2020). However, outcrops are commonly limited by their lateral extent, and studies of multiple outcrops within the same system are frequently beset with correlation problems that make it difficult to understand the spatiotemporal evolution of the system in detail (Nichols, 1987;Turner et al., 1984). Flume tank experiments have also been used to better understand fluvial sedimentology at a system scale (Connell et al., 2012;Terwisscha van Scheltinga et al., 2020) but have associated scaling issues. Remote sensing, outcrop and flume tank studies can be complemented by numerical models, which can reproduce the behaviour and resultant stratigraphic architecture of fluvial systems at temporal and spatial scale. This has been achieved for deltas (Huang et al., 2015;van der Vegt et al., 2016), alluvial fans (Clevis K E Y W O R D S distributive fluvial system, Ebro Basin, fluvial fan, Huesca system, process-based forward model, stratigraphic forward model | 3139 EAGE et al., 2003;Welivitiya et al., 2020), single-channel belts (Karssenberg & Bridge, 2008) and parts of sedimentary basin fills (Bridge & Leeder, 1979;Yan et al., 2017).
We suggest that linking governing parameters (initial surface, hydrology and sediment characteristics) to the resulting topology of architectural elements leads to a predictive stratigraphic model of DFS. We show that this forward loop can be closed by using stratigraphic forward modelling through the following steps: Find a suitable DFS in outcrop to model, find modern systems for the governing parameters that cannot be derived from outcrop, build a stratigraphic forward model based on these measurable input parameters, use an iterative modelling workflow to adapt model input parameters not based on measurements until the modelled output closely aligns to the observed sedimentological parameters of the DFS in outcrop and finally use the resulting stratigraphic forward model of the DFS to improve our understanding of this sedimentary system.
A hydrodynamic stratigraphic forward modelling approach was chosen as it can best approximate the sedimentary processes that are crucial for the development of DFS.
The model was quality checked against the Huesca DFS in northern Spain, which is well studied with quantified sedimentological data encompassing a representative proportion of the system (Burnham, 2016;Hirst, 1991;Kulikova, 2013;Martin et al., 2021). A successful model is deemed to be the one where output data closely align with sedimentological observations such as grain size distribution and channel thickness at the respective Huesca DFS outcrop locations as well as reproducing the sedimentary features described in the Huesca DFS. Once validated, the model allows the grain size distribution and bed thicknesses to be sampled at the true temporal and spatial scales. SFM allows the prediction of stratigraphic architecture in three dimensions enabling the quantification of sandbody connectivity, the result of longitudinal, lateral and temporal channel and splay contacts. Sandbody connectivity constitutes the primary control of reservoir performance with implications for subsurface exploitation and storage of fluids in DFS deposits.

| Description of the outcrop system
The Huesca DFS is located in the Ebro Basin, which is part of the southern Pyrenean foreland basin, northern Spain ( Figure 1). The Pyrenees formed in the Late Cretaceous due to collision between the Eurasian plate to the north and the Iberian sub-plate to the south (Muñoz, 1992;Puigdefàbregas et al., 1992;Teixel, 1996). The Ebro Basin was initiated in the Late Eocene (Costa et al., 2010) but did not become a unique entity until the Oligocene when F I G U R E 1 Overview map of the central Ebro Basin and the Huesca distributive fluvial system. Arrows indicate the mean drainage direction derived from paleocurrent measurements (Hirst, 1991;Nichols, 1987). The outcrop locations from (Martin et al., 2021) are in blue. The map below shows the location of the Ebro Basin and the upper map as a red square within the Iberian Peninsula (modified from Hirst & Nichols, 1986) the southward migration of the Pyrenean thrust front resulted in uplift of the older Jaca Basin and emergence of the External Sierras (Cámara & Klimowitz, 1985;Muñoz, 1992;Vergés, 1993). The Ebro Basin has a triangular shape and is endorheic, bound to the north by the External Sierras, to the south by the Iberian Range and to the east by the Catalan Coastal Range (Figure 1). It remained internally drained until the Late Miocene (Garcés et al., 2020) aggrading kilometres of sediment (Gaspar-Escribano et al., 2001). At that time, headward erosion of the Proto-Ebro River breached the Catalan Coastal Range and connected the Ebro Basin to the Mediterranean Sea. This caused the switch from aggradation to erosion, and the present-day landscape started to form (Garcés et al., 2020).
The Huesca DFS comprises part of the Sariñena Formation (Soler-Sampere & Puigdefábregas, 1970), which has been dated using microvertebrate fossils as Chattian to Burdigalian (Luzón, 1997(Luzón, , 2005 of which the initial 250 kyr has been modelled here. It encompasses fluvial and overbank deposits, which terminated in a playa lake as a series of unconfined splay deposits (Hirst, 1991). Recent work by Martin et al. (2021) combined with mapping conducted by Hirst (1991) and Vincent and Elliott (1996) positioned the apex of the DFS close to the village of Naval (Figure 1). Outcrops of the Huesca DFS are limited to the sector west and southwest of the Barbastro Anticline ( Figure 1). The proximal zone has been tectonically altered, most notably through the rise of the syndepositional Barbastro anticline (Martínez Peña & Pocoví Juan, 1988), but the remainder of the DFS is largely undeformed (Hirst, 1991). The Huesca DFS has been divided into proximal, medial and distal zones (Hirst, 1991), and this division was also used for the results of the SFM.
The proximal zone occupies a radius of 20 km downstream of the apex (Hirst, 1991;Martin et al., 2021) and is characterized by amalgamated pebbly and sandy channel deposits of braided rivers without preservation of overbank facies. Cross-bedding and preservation of bar structures points towards highly mobile fluvial channels, migrating within large channel belts. Sandbodies are fully connected in the proximal zone, due to the fully amalgamated multistorey channel architecture and the absence of laterally extensive overbank fines (Nichols & Fisher, 2007).
The medial zone occupies a radius of approximately 20-60 km downstream of the apex (Hirst, 1991;Martin et al., 2021) and is dominated by channel belt deposits with both ribbon and sheet geometries, enclosed within overbank mudstones. The channel deposits are predominantly composed of sandstone with sparse mudstone (Hirst, 1991). Sandbodies may be connected laterally, but vertical connectivity is less common due to the presence of intervening overbank fines, which are only locally removed by erosion (Nichols & Fisher, 2007).
The distal zone extends in a radius beyond 60 km from the apex (Hirst, 1991;Martin et al., 2021) and is dominated by overbank fines with occasional thin sheet sandstones formed by unconfined flow and sparse ribbon (isolated) channels. These fines interfinger with lacustrine muds, marls and limestones (Hirst, 1991). Sandbody connectivity in the distal zone is observed in outcrop to be lower than in the medial deposits as the smaller volume sandbodies connect neither laterally nor vertically (Nichols & Fisher, 2007).
Downstream trends were quantified in the Huesca DFS by Hirst (1991) and Martin et al. (2021): average grain size decreases from very coarse sand to fine sand and channel body percent decreases from 66% to 4%. Although these trends vary, most notably in the proximal and medial zone, no discernible vertical changes in the fluvial architecture were observed (Hirst, 1991;Martin et al., 2021). This indicates that the outcrops only encompass the aggradational phase of the Huesca DFS.
Surrounding the toe of the DFS are lacustrine deposits composed of fine clastic, carbonate and evaporite deposits (Hirst & Nichols, 1986). A sediment accumulation rate of approximately 85 mm/kyr was estimated from magnetostratigraphy (Garcés et al., 2020;Pérez-Rivarés et al., 2002). This rate was derived from deposits approximately 80 km southwest of the Huesca DFS apex.
The quantified outcrop study by Martin et al. (2021) was used to validate the model output. This study is based on the outcrops Pertusa, Monzón, Montearagón, Piraćes, Castleflorite, Torrollón, Bolea and Sigena ( Figure 1). They are in the medial to distal zone of the system (distance from the apex of 31,32,43,43,44,45,59 and 62 km, respectively).

| Stratigraphic forward modelling software
The stratigraphic forward modelling software Sedsim was used to link the governing parameters (model inputs) to the resulting stratigraphic architecture (model outputs). Sedsim is a deterministic hydrodynamic, threedimensional, process-based program. It is fully scalable in time (seconds to millions of years) and space (centimetres to kilometres; Tetzlaff & Harbaugh, 1989). Sedsim was selected over other software for its hydrodynamic approach, which can model all physical processes needed for the formation of a DFS. It was also chosen as measured data (e.g. hydrological parameters from river gaging stations) can be directly used as model input parameters. In contrast, other software use input parameters unrelated to real-world data which must be fine-tuned to match the desired model output. A description of all process-based forward models tested for this project and why Sedsim was selected is detailed in the Appendix S1.
Sedsim uses a computational approximation to the Navier-Stokes equation to compute the hydrodynamics (Griffiths et al., 2012;Tetzlaff & Harbaugh, 1989). Open channel flow is modelled in two horizontal dimensions with a depth-averaged flow velocity. To approximate an actual depth-velocity profile with velocity decreasing towards the water-channel interface, a frictional parameter (Manning Coefficient) is used at this interface. Continuous flow is represented by fluid elements, fixed volumes of water, which are released at the source at each time step. The fluid elements move between nodes (cells have nodes at each corner) following the steepest downhill gradient of the underlying surface at that time step. Sediment erosion, transport and deposition are modelled in three dimensions and computed using a mixed Lagrangian/Eulerian computational scheme for every node at each time step. These processes are modelled according to the principle of mass preservation. The boundary between erosion and transportation is determined by the critical shear stress, calculated as a function of particle diameter (Griffiths et al., 2012). Up to four clastic sediment types can be included as well as carbonates and organics. The model output, grain size distribution and sediment thickness are stored at each node for each display interval (Tetzlaff & Harbaugh, 1989).

| Model input parameters
The simulation parameters are summarized in Table 1 and described in more detail in the Appendix S1. The workstation used for modelling was running Windows and had an Intel Xeon Gold 5122 processor with four cores, 3.6 GHz and 64 GB of installed memory.
Time: The time parameters (Table 1) were chosen as a compromise between temporal resolution and computer run time. An interval of 250 kyr within Aquitanian to Burdigalian times (23-16 Ma) was modelled within the lacustrine sediment accumulation rates (Garcés et al., 2020;Pérez-Rivarés et al., 2002). The computational time step (Table 1) was kept as short as possible given the available computing power. The temporal resolution of the model is determined by the display interval, the time interval at the end of which the software saves the model output. The model output file sizes were limited by the post-processing software to 1,000 display intervals. Thus, a display interval of 250 years was chosen (250 kyr model duration/1,000 display intervals). Grid: The model size (Table 1)   10 km. The grid resolution (Table 1) was a compromise between temporal resolution and computer run time. The initial surface ( Figure 2a) is a fan building out from the apex. The fan has a concave profile with a 0.002 gradient at the apex, decreasing to an 0.0002 gradient at the toe ( Figure 2b). Downstream of the fan is a lacustrine slope with a 0.0002 gradient. Depositional gradients could not be determined in outcrop as the tectonic reconstruction and measurement errors are too great. The gradients used were, thus, derived from average gradients (measured from digital elevation models) of modern DFS . These DFS are located in northern India and are active in a similar climate zone and basin type and have a similar downstream length to the Huesca DFS. The flattening of the gradient downstream of the apex represents the graded profile of fluvial systems (Davis, 1806;Gilbert, 1877). The modelled fluvial input is at the top of a valley, which cuts through the high ground to the north and widens towards the apex ( Figure 2a). All model boundaries besides the boundary upstream of the source (to the north) were open to fluid and sediment transport out of the model. Lake level: During the modelled period, the lake level rose to keep pace with the lacustrine sediment accumulation rate (Garcés et al., 2020;Pérez-Rivarés et al., 2002) so that the lake water depth remained constant. This is used as a proxy for subsidence.
Sources: Of the source parameters (Table 1), only the source location could be derived from outcrop data and was set 10 km upstream of the DFS apex location (Figure 2a). For the source activity periods and the hydrodynamic parameters (flow velocity, discharge and sediment concentration) measurements from the modern High Island Creek gaging station in Henderson, Minnesota, USA (Groten et al., 2016) were used. This creek was chosen as it encompasses a similar grain size range (Groten et al., 2016; F I G U R E 2 (a) Initial surface for the stratigraphic forward model of the Huesca distributive fluvial system. The distributive fluvial system is shown in yellow and the lake in blue. North is up and the vertical exaggeration is 100 times. The cross-section in Figure 4 is along the red dashed line named 'Southern section'. The cross-sections in Figure 5 are along the black lines named 'Proximal section', 'Medial section' and 'Distal section'. (b) Cross-section along 'Southern section' in (a) with the slope gradients in red Hirst, 1991;Martin et al., 2021) and its watershed is in the same climate zone as the Huesca DFS (Groten et al., 2016;Hamer et al., 2007;Kottek et al., 2006). The High Island Creek sediment load was matched to the uncompacted sediment volume deposited in the Huesca DFS during the model duration (see Appendix S1). To create realistic flow regimes, only flood events ranging from annual to decennial floods were modelled. The occurrence of each flood event and the magnitude of the hydrological parameters were randomized (see Appendix S1).
Sediments: Four representative grain sizes (Table 1) were selected to cover the grain size range seen in the Huesca DFS outcrops (Hirst, 1991;Martin et al., 2021): pebbles (P), coarse to medium sands (CMS), fine sands and silts (FSS) and clays (C). Outcrop data were not used for the grain size distribution at the source as they represent a biased sample of the entire system. The grain size distribution from the High Island Creek data (Groten et al., 2016) was used instead. A petrographic study for outcrops near Pertusa (Turner et al., 1984) on sediments ranging in grain size from medium sand to coarse silt showed that the main framework minerals are quartz, clay and oxide lithoclasts and limestone lithoclasts. In this petrographic analysis, carbonate lithoclasts and cements (calcite) were not distinguished. As the samples were heavily cemented, it was assumed that most of the carbonate content was introduced after deposition. Diagenesis was not included in the numerical model. Only quartz and carbonate (lithoclasts and cement) content was measured by Turner et al. (1984) and quartz volume ranged from 30% to 60%. Given that clastic quartz predominates, a density of 2,650 kg/m 3 (Anthony et al., 2015) was used for the three coarsest grain sizes (pebbles, CMS, and fine sands to silts) in the model. The clay minerals found in the Huesca system are mainly chlorite and illite (Turner et al., 1984). As no known quantitative petrography of clays has been performed on the sediments in the Huesca DFS, a clay grain size density of 2,650 kg/m 3 was assumed, falling within the range of chlorite and illite densities measured in labs (Anthony et al., 2015). The transport type was taken from the High Island Creek as well: the two coarser grain sizes were transported as bed load and the two finer grain sizes were transported as suspended load (Groten et al., 2016).
The modelling process is iterative starting with a simple model using only the mandatory Sedsim input modules (time, grid, sea level, sources and sediments). Additional modules were added incrementally (slope angles, source heights, parametric sampling interval, sediment transport parameters and Manning coefficients) and, as understanding of the depositional system advanced, the input parameters for the source module were made more complex, varying the number of sources, their duration of activity and their parameters. The optional Sedsim modules are detailed in the Appendix S1.
During the modelling process, the most sensitive input parameters proved to be the shape of the initial surface and the sediment source parameters. Adding source variability as a proxy for climate oscillations created a fluvial pattern that matched observations of a DFS more closely. The variability was fine-tuned to match the outcrop data while keeping the total sediment volume constant.

| Model output processing and visualization
The model outputs were processed to extract grain size distribution, sediment thickness, net to gross (NTG), sandbody thickness, number of sandbodies and amalgamation ratio. NTG and sandbody thickness were also compared with outcrop data from Martin et al. (2021). The detailed description of the processing and visualization can be found in Appendix S1.
NTG is a key consideration in hydrocarbon reservoir studies, and there are numerous definitions (Ringrose & Bentley, 2015). In simple terms, it is the proportion of the rock volume that has the potential to be reservoir. For this analysis, we consider net to encompass all grain sizes besides clay.
Amalgamation ratio is an important characterization besides NTG for reservoir performance as it quantifies vertical connectivity. A ratio of 1 indicates a fully connected reservoir without baffles or barriers, whereas a low ratio indicated a highly compartmentalized reservoir in the vertical direction. The amalgamation ratio is defined as the fraction of sand-to-sand contacts (amalgamation surfaces) relative to all bed/cell contacts (sand-to-sand, sand-to-mud and mud-to-mud) in a vertical pseudo-well (Zhang et al., 2017). It was calculated as (sum of cell contacts within each sandbody)/(total number of cells -1).
To visualize the three-dimensional data in a twodimensional section, the Euclidian distance between each node and the apex was calculated and all nodes with the same distance from the apex were combined to 1 point along a proximal to distal radius. The equidistant nodes were grouped in 1 by 1 km bins as this equals the spatial resolution of the model. Only the nodes south of the apex were sampled (Figure 2).

| Model limitations
All models are a simplification of the real world and as such have limitations. The main limitations of this model stem from temporal and spatial resolution, using multiple data sets as input for the model and inherent limitations of the stratigraphic forward modelling software Sedsim.
The key limitation is in the model resolution, which is limited by the available computing power. The hardware used for this project allowed a 500 × 500 m grid size and a 5-year time step. Following the Nyquist-Shannon theorem (Nyquist, 1928;Shannon, 1949) only spatial features with a minimum length/width of 1 km can be resolved. That grid size is larger than the width of many features such as the fluvial channels observed in the Huesca DFS, and consequently the modelled channel bodies are too wide. The temporal resolution is limited by the same theorem to a minimum of ten years. The hydrological events modelled in this project are floods, and their duration is measured in days. Using 5-year time steps upscales these floods. Models with shorter time steps show a greater degree of detail in the depositional pattern but more important than the time step is the display interval as the model output is computed for every time step but averaged and then saved at every display interval. The maximum number of display intervals allowed by the post-processing software is 1,000, which for a model duration of 250 kyr results in display intervals no shorter than 250 years. A shorter time step would ultimately not result in a higher resolution if it is shorter than half the duration of the display interval (Nyquist, 1928;Shannon, 1949).
The model input parameters are a combination of model dimensions and sediment characteristics from the Huesca DFS, hydrodynamic parameters from the High Island Creek and the initial surface topography from modern DFS in northern India. As the goal of the project was to link controls on sedimentary processes (model input parameters) to the resulting deposits (model output), it was paramount to use as many measured parameters for the model input as possible.
Sedsim only models fluid flow in two dimensions with a depth-averaged velocity profile to approximate the vertical component (Griffiths et al., 2012;Tetzlaff & Harbaugh, 1989). Because of this simplification, sedimentary features that are thinner than the flow depth such as dunes or bars within the channels are not modelled. Inter-grain cohesion is also not modelled in Sedsim (Griffiths et al., 2012;Tetzlaff & Harbaugh, 1989), which means that clay-rich deposits in the model are more easily erodible and do not maintain steep depositional angles. This might increase channel mobility and avulsion rates in the medial zone where floodplainchannel interactions are common.  Figure 2a. The location of the apex (A) and distances from it in kilometres are marked in red. The lake level is shown in blue. The grain size distribution between two nodes of the same time step is interpolated and given the corresponding colour out of the quadrangle of colours with the end members shown in the legend. The black square in 250 kyr outlines the location of Figure 12 Subsidence was not modelled as, to the authors' knowledge, no subsidence rate was measured for the Huesca DFS. The model was deposited by allowing the system including the source point to aggrade. In Sedsim as in reality, rivers progressively downcut to the base level, in this case the lake level. As aggradation and not downcutting was observed in the rock record (Hirst, 1991;Martin et al., 2021) the lake level was raised during the model duration to prevent downcutting and allow aggradation. This only represents one of many possible subsidence-lakelevel interactions, which could have taken place during the deposition of the Huesca DFS.
The uncertainty introduced by model limitations was not quantified but for each model iteration, the goal was to reduce the effect of these limitations on the model output. This was done by quality checking each model iteration against the Huesca DFS outcrop data with the goal of achieving an ever-closer match between model and outcrop data.
Quantifying the uncertainty linked to temporal and spatial resolution was done in preliminary runs on cropped parts of the model (only part of the entire system and model duration) as the final model uses the highest temporal and spatial resolution possible for the given computer capacity. SFM need to model the entire system to be able to reproduce all the sedimentary processes observed. If only a small part of the system is modelled, not all sedimentary processes can be reproduced. Therefore, these cropped models to test the effects of different temporal and spatial resolutions cannot be directly compared with the final model and the uncertainty introduced by resolution not quantified.
To reduce the uncertainty introduced by combining three data sets as model input parameters, the data sets were taken from areas with the same climate and tectonic setting to reduce mismatch and during the modelling process. Optional modules in the stratigraphic forward modelling software were used to match the three data sets. Unavoidably, inherent uncertainty remains when combining multiple data sets into a single model. A sensitivity analysis on the impact of the different input data sets on the model output were beyond the scope of this project.
It was impossible to conduct a direct and quantified comparison between the tested process-based forward modelling software as each software takes different input parameters and uses different physical equations. That is why the uncertainty introduced due to the inherent limitations of Sedsim could not be quantified.

| RESULTS
During the model duration, 550 km 3 of sediment was introduced at the source of the model and 450 km 3 of those sediments were ultimately deposited within the bounds of the model. The remaining 100 km 3 were moved out of the model through the open-model boundaries to the east, south and west. All the Ps and CMS, 98% of the FSS and 75% of the C were retained in the model. The deposited volume equals 1.64 times the estimated minimum uncompacted sediment volume of the Huesca DFS over the same time period.

| Stabilisation and progradational phase of the model
The initial surface (Figure 2a) was devoid of sediments, and initially CMS were deposited within channels from the source to the lake shoreline. FSS and C were mainly deposited not only at the lake shore in terminal fans F I G U R E 4 Cross-section of the modelled grain size distribution from north to south through the source and apex. The depositional surfaces after 20, 125 and 250 kyr are marked as red lines. The colour scheme is the same as in Figure 3. The location of the cross-section can be seen in Figure 2a but also adjacent to the apex as sheet flows (Figure 3 -0.25 kyr). For the first 20 kyr, this resulted in a finingupward grain size distribution due to retrogradation of terminal fans over the initial sand bypass zone all the way to the proximal part of the DFS. Coarsening-upward successions are limited to the proximal zone due to the progradation of sediments adjacent to the apex (Figure 3 -0.25, 20 kyr, and Figure 4). After approximately 20 kyr the retrograding and prograding deposits coalesce and blanket the entire area later occupied by the DFS. This marks the stabilisation of the model and the beginning of DFS deposition (Figure 3 -20 kyr).
From approximately 20-125 kyr, the stable DFS progrades (Figure 4), as shown by the coarsening and thickening upward trends ( Figure 5). In plan-view the progradations relate to active fluvial channels and travel  Figure 3. The location of the cross-section can be seen in Figure 2a further downstream than the sediments not confined to the active channel (Figure 3 -20 and 125 kyr). In a downstream cross-section, this is shown as multiple progradation cycles with each cycle occurring when the active fluvial channel flows along the selected downstream section (Figure 4).

| Aggradational phase of the model
After 125 kyr of modelled time, the proximal to medial transition has reached its final distance of 20 km from the apex and thereafter remains stable in space (Figure 3 -125 and 250 kyr) even though oscillations around the average distance still occur, which is expected within an avulsion driven system (Figure 4). The radius of the modelled DFS is continuously extending in the aggradational phase but mostly through deposition of C. The coarser grain sizes are not deposited any further downstream after 125 kyr, which can be seen in the aggradational stacking pattern of the proximal, medial and distal deposits ( Figure 5).
The proportion of P in the proximal zone (0-20 km downstream of the apex) decreases from 98% at the apex to 2% at the boundary between the proximal and medial zones ( Figure 6). In plan-view, grain size P is mainly deposited in the feeder channel connecting the source to the apex and two lobes on either side of the apex. Only at the sides and downstream end of the feeder-channel, finer-grained sediments especially CMS and C were deposited. These deposits lead to a decrease in in NTG, sandbody thickness and number and amalgamation ratio. They also greatly increase the variability in these parameters, especially sandbody thickness and number (Figure 7). Further downstream, the depositional pattern changes to dendritic and CMS is the dominant grainsize and makes up 2% at the apex and 80% at the boundary between proximal and medial zones (Figures 6 and 8). The stratigraphic architecture is characterized by amalgamated channelized deposits (centre of Figures 4 and 5a,b) and sheet flow deposits (Figure 9). The sandbodies are massive in planform ( Figure 5 -250 kyr). The average total sediment thickness is stable throughout the zone between 32 and 33 m with little spread of measured thicknesses. This results in a sediment accumulation rate of 260 mm/kyr (Figure 10b). NTG decreases from 100% at the apex to 88% on average at the transition to the medial zone, whereas the variability in NTG increases strongly downstream of 16 km from the apex. It reaches its maximum at the transition between proximal and medial zone of 51% (range: 49%-100%; Figure 7a). For the first 6 km from the apex, the sediments are mostly part of a single sandbody with a thickness between 31 and 36 m. From 7 to 20 km, the number of sandbodies increases to 20 (at 19 km from the apex). The maximum sandbody thickness remains constant around 35 m, but the minimum thickness drops to around 0 m while the average thickness is 12 m at 20 km (Figure 7b,c). When the number of sandbodies increases dramatically starting at 14 km from the apex, the amalgamation ratio starts to drop from around 1 to 0.87 at 20 km from the apex. The range of amalgamation ratios increases at a faster rate to span 0.34 to 1 at 20 km from the apex (Figure 7c,d).
In the medial zone (20-60 km downstream of the apex), the proportion of CMS is reduced from 80% at the proximal to medial boundary to 5% at the transition to the distal zone. P is only present in a small percentage and disappears by 30 km from the apex. The proportion of FSS increases from 4% at the proximal to medial boundary to 15% at the boundary between the medial and distal zones. Most of the medial zone is comprised of C ( Figure 6). In plan-view, the grain size distribution shows a similar trend but highlights the variability of grain size distribution with a distance from the apex (Figures 3 and  8). Although the fan is radial, the CMS propagate approximately 10 km further downstream towards the south than towards the west. (Figure 8b-d). The averaged sediment thickness increases from 32 m at the proximal to medial boundary to 35 m at 48 km from the apex and then reduces to 32 m at the medial to distal boundary. This results in a sediment accumulation rate of 266 mm/ kyr. The variability of thickness measurements with similar distances from the apex is much higher in the medial zone and reaches up to 35 m (Figure 10b). In plan-view, the greatest thicknesses are located along channelized deposits (Figure 10a). The cross-section through the medial zone (Figure 6c) shows flood plain deposits of FSS and C separating the CMS and FSS filled channelized deposits. The channel deposits are incised into the underlying F I G U R E 6 Plot of averaged and stacked grain size distribution vs. distance from the apex. The modelled data comes from the aggradational time interval (125-250 kyr). All nodes with the same distance from the apex are averaged to a single grain size distribution and plotted at the respective distance from the apex sediments, and there are also crevasse splays of CMS to C interbedded with the FSS and C (Figure 9). Some channels are plugged by FSS and C, with CMS only preserved on accretion surfaces created by lateral migration of the channel (Figure 5c). The sandbodies are channelshaped in planform (Figure 3 -250 kyr). In the medial zone, the average NTG drops from 88% at the boundary with the proximal zone to 19% at the boundary with the distal zone. The variability of NTG values reaches its maximum in the first few kilometres of the medial zone but decreases downstream to around 40% (Figure 7a). The maximum number of sandbodies are encountered at 25 km from the apex and decrease steadily until the medial to distal transition (Figure 7c). Around 25 km from the apex, the sandbody thickness also decreases at a much lower rate. The average thickness only reduces by about 2.4 m from 25 to 60 km from the apex, whereas in the transition zone (15-25 km from the apex) the average thickness decreases by 25.4 m. The variability in thickness decreases downstream as well (Figure 7b). The amalgamation ratio also decreases steeply until 25 km from the apex and then more gradually until 60 km from the apex to reach 0.04. As with the number of sandbodies, the amalgamation ratio range is largest around 25 km from the apex and decreases downstream (Figure 7d).
In the distal zone (60-70 km downstream of the apex), C dominates and FSS decreases continuously until it disappears by 85 km from the apex. CMS is only present as a few percent and disappears by 80 km downstream (Figure 3). In plan-view, the grain size distribution shows a clustered texture (Figure 8b,c). The sediment clusters occur in areas with reduced sediment thickness (Figure 10a). These sediment accumulations are situated at the end of channels and represent terminal fans (Figure 3 -125, 250 kyr and Figure 9). This creates patchy and isolated sandbodies in plan-view (Figure 3 -250 kyr). The stratigraphic architecture reveals only isolated channelized deposits mostly filled with FSS and to a lesser extent CMS and terminal fans mainly depositing C (Figures 5d and 9). The average total sediment thickness in the distal zone drops from 32 to 29 m (Figure 10b). This results in a sediment accumulation rate of 241 mm/kyr. NTG decreases from 19% to 13% in the distal zone (Figure 7a). Sandbody thickness F I G U R E 7 Plots of (A) net to gross, (B) sandbody thickness, (C) number of sandbodies, and (D) amalgamation ratio vs. distance from the apex. The modelled data comes from the aggradational time interval (125 to 250 kyr). The parameters at every node with the same distance from the apex are averaged to a single value and plotted at the respective distance from the apex. The bars represent the minimum (bottom) and maximum values (top) to show lateral variability and number also decrease from 0.6 to 0.4 m and 2.5 to 2, respectively (Figure 7b,c). The amalgamation ratio decreases from 0.04 to 0.02 (Figure 7d). Variability of all parameters decreases further (Figure 7).
In the lacustrine zone (more than 70 km downstream of the apex), only FSS and C are present, but deposition does not occur over the entire model area covered by the lake. The average total sediment thickness drops from around 29 m at the lake shore to below 1 m beyond 101 km away from the apex. This results in a sediment accumulation rate of 155 mm/kyr. NTG drops from 13% to 0% on average at 85 km from the apex. The sandbody thicknesses and numbers reach zero at the same distance, but the amalgamation ratio already reached 0 on average at 79 km from the apex (Figure 7).  (Figures 5 and 6) and channel confinement decreases distally (Figure 9).
At a system scale, the Huesca DFS has been divided into proximal, medial and distal zones (Hirst, 1991): The proximal zone is characterized by amalgamated pebbly and sandy channel deposits without preservation of overbank fines (Hirst, 1991). In the model, the grain size distribution is also dominated by pebbly and sandy deposits but downstream of 15 km from the apex, overbank fines are deposited as well as increasing to 12% at the proximal to medial transition ( Figure 6). No grain size distribution F I G U R E 8 Maps for the modelled grain size distribution averaged at every node over the aggradational time interval (125-250 kyr) with the distribution of (a) pebbles, (b) coarse to medium sands, (c) fine sands to silts and (d) clays. The unit of these maps is the fraction of each grain size to the total grain size distribution is available from outcrops in the proximal zone of the Huesca DFS, so overbank fines adjacent to the medial zone could be present but not described from the outcrop. Strong channelization was reproduced by the model (Figure 5a), but sheet flows were also modelled (Figure 9), which are not described in outcrop. This results from a minor mismatch between the depositional gradient and the hydrodynamic input parameters, which could not be eliminated during the modelling workflow. The proximal zones of the Huesca DFS was described as a single sandbody (Nichols & Fisher, 2007). Figure 7c shows that in the model most of the proximal zone is encompassed by a single amalgamated sandbody that starts to become separated by floodplain deposits towards the proximal to medial transition. As with the grain size distribution, the sandbody connectivity in the transition zone is not specifically described in the Huesca DFS.
The medial zone in the Huesca DFS is described as sandy channel belt deposits enclosed within overbank mudstones (Hirst, 1991), which corresponds closely to the modelled grain size distribution ( Figure 6) and channel pattern (Figure 3 -125, 250 kyr and Figure 9). Subordinate channels with mud plugs were observed in outcrop (Hirst, 1991) and also reproduced in the model (Figure 5c). Channel-fill, ribbon and sheet sandstones are described in the outcrop (Hirst, 1991). In the model, channel-fill sandstones can be seen in the east-west crosssection (Figure 5c), and unconfined sheet sandstones are highlighted in Figure 9. Ribbon sandstones were not observed in the model that might be due to the coarse grid resolution. The sandbody connectivity observed in outcrop is mostly lateral with limited vertical connectivity (Nichols & Fisher, 2007). In the medial zone of the model, the rapidly decreasing amalgamation ratio (Figure 7d) clearly shows a low vertical connectivity of sandbodies. In cross-section, horizontal connectivity in the model is, however, apparent (Figure 5c).
Overbank fines with occasional thin-sheet sandstones characterize the distal zone of the Huesca DFS (Hirst, 1991). The modelled grain size distribution matches this description (Figure 6), and most channelized deposits terminate in fans within the distal zones ( Figure 9). Strong flood events transport CMS as thalweg deposits or proximal parts of terminal fans and create isolated channel deposits. These channel deposits are mainly filled by C and FSS (Figure 5d). In outcrop, sandbodies are not connected laterally nor vertically and are of small volume. The model has a very low amalgamation ratio and sandbody thickness (Figure 7b,d), and in cross-section, the sandbodies are mostly isolated (Figure 5d).
Deposition in the distal lacustrine zone is oversimplified in the model as it was not a focus of the study. The deposits are sourced from only one fluvial system, whereas the lacustrine deposits of the Ebro Basin were a combination of clastic influx from multiple rivers, carbonates and evaporites (Garcés et al., 2020). Even though the Sedsim stratigraphic forward modelling software can recreate all these additional sedimentary systems, they were not included in this study as it focussed on a single siliciclastic fluvial system. This means that the interfingering of distal DFS and F I G U R E 9 Surface-view of the southwestern section of the model at 250 kyr. The outline is shown in Figure 3 -250 kyr. The vertical exaggeration is 100 times. The model is shown in perspective view. For scale, the location of the apex (A) and distances from the apex in kilometres are marked in black. Red indicates the depositional surface of the previous time step. Green indicates the pebble and coarse to medium sand grain sizes deposited during the time step. Blue indicates the fine sand to silt and clay grain sizes deposited during the time step. (1) highlights the sheet flow of pebbles and coarse to medium sands in the proximal zone. (2) indicates multiple crevasse splays of coarse to medium sands (green) and clays (blue). (3) is a terminal splay mainly depositing clays and some fine sands and silts (blue) as well as coarse to medium sands (green) closer to the feeder channel. The channel in green is incised from the proximal to medial transition up to the beginning of the terminal fan lacustrine deposits seen in outcrop (Hirst, 1991) is not modelled. Not modelling evaporite and carbonate sediments in the lacustrine zone also led to a sediment thickness rapidly decreasing away from the lake shore (Figure 10b) unlike the shallow depositional gradient of lacustrine deposits observed in the Ebro Basin (Valero et al., 2014). The lake level in the Ebro Basin varied due to Milanković-driven climatic cycles at a million-year scale (Valero et al., 2014), which was not included in this project as the model duration is much shorter than the cycle length.
To compare the model to outcrops, the model output parameters (Figure 11a,b) were extracted at the equivalent outcrop locations and compared with outcrop measurements ( Figure 1) from Martin et al. (2021). To quantify the comparison, a ratio between the model parameters and the outcrop measurement was calculated (Table 2). For the NTG/in-channel component, the match between the model and outcrop for Monzón, Montearagón, Piraćes and Torrollón is good (model-outcrop ratio: ±15%), for Pertusa and Castleflorite is intermediate (model-outcrop ratio: between ±15% and ±50%) and for Bolea and Sigena is poor (model-outcrop ratio: more than ±15% and ±50%). On average the NTG is 1.5 times higher in the model than in the outcrops. For the sandbody/channel thickness, the F I G U R E 1 0 (a) Map of the modelled sediment thickness deposited during the aggradational time interval (125-250 kyr) for every node. (b) Plot of total sediment thickness deposited during the aggradational time interval (125-250 kyr) vs. distance from the apex. The total sediment thickness at every node is shown as grey crosses and the averaged thickness (all nodes with the same distance from the apex are averaged to a single value) is shown as a black line match between the model and outcrop is generally poor as only Monzon has a model-outcrop ratio below ±50%. On average, the modelled sandbody thickness is one third the channel thickness measured in outcrop. In summary, the match between the model output and the outcrop measurements at the outcrop locations is poor.
Although the model cannot recreate the Huesca DFS exactly, it aims to create a modelled DFS, which matches the outcrop measurement at their respective distances from the apex. Therefore, in Figure 11a,b, the minimum and maximum values encountered at the equivalent distance from the apex to the outcrops was also plotted for each outcrop location. If the outcrop measurements fall within the minimum to maximum range of the model, the model is a good approximation of the system-scale distribution of grain size and sediment thickness. Figure 11a shows that the NTG/in-channel component measurements from the outcrop always plot within the minimum to maximum range of the model at equidistant locations. The mismatch between the average NTG of the model at equidistant locations and the NTG/in-channel component from the outcrop as well as the NTG from the model at the outcrop locations shows that there is a strong lateral variability within a DFS, especially at the Bolea and Sigena outcrops where NTG measured at the outcrop is much lower than extracted from the model. In the model at the outcrop equivalent location more channel and splay deposits resulted in a higher net deposition than observed in the outcrop.
As sandbody/channel thickness can have multiple measurement at a single location, the minimum, average and maximum values from the outcrop (red), the model at the outcrop locations (green) and the model at equidistant locations (black) were plotted (Figure 11b). The average channel thickness from the outcrop lies within the minimum to maximum range of the modelled sandbody thickness at equidistant locations. The average channel thickness measured in outcrop is larger than the sandbody thickness modelled at equidistant locations. This could be due to the lateral variability and/or indicate that the modelled sandbodies are too thin (by a factor of 3). Another reason for this mismatch is that the comparison is between sandbodies and channels. Sandbodies encompass thinner sheet sandstones and channel sandstones which can lead to a thinner average sandbody thickness. Only in Pertusa and Piraces is the maximum thickness of channel body measured in outcrop thicker than the maximum thickness modelled in equidistant locations (0.2 and 2 m, respectively).

| How does the model increase our understanding of the DFS in the rock record?
DFS tend to form a radial pattern in planform if they are not laterally confined ). This DFS model uses a perfectly semi-circular initial surface and is not laterally confined downstream of the apex (Figure 2), F I G U R E 1 1 Comparison of model output (in green) and outcrop measurements from (Martin et al., 2021; in red) at eight outcrops. (a) Comparison of modelled net to gross and net to gross (virtual outcrops) and in-channel component measured in the outcrop (other outcrops). To show the lateral variability in the model, for each equivalent outcrop location all equidistant nodes (same distance from the apex) were sampled to give a maximum (upper bar), average (circle) and minimum (lower bar) value (in black). (b) Comparison of modelled sandbody thickness and channel thickness measured in the outcrop. The maximum (upper bar), average (symbol) and minimum (lower bar) for the modelled equidistant nodes, for the modelled nodes at the equivalent outcrop locations and for the actual outcrop locations was plotted. The modelled data comes from the aggradational time interval (125 to 250 kyr). For outcrops Pertusa, Piraćes, Torrolón, Bolea and Sigena, the parameters were measured in virtual outcrops. For the other three the parameters were measured in the outcrop so a radial pattern in planform should result. When plotting the stacked grain size distribution as in Figure 6 but for transects from the apex to the eastern, southeastern, southern, southwestern and western corner (Figure 12), a similar grain size distribution with distance from the apex is observable but the greater the angle from the northsouth median, the greater the difference in grain size distribution from the radial average one shown in Figure 6.
The eastern and western transects have a 2-km-wide section close to the apex with a 60% and 69% reduction of P, respectively, and replacement by CMS and C. The map view (Figure 8) shows that this zone is situated at the edges of the canyon connecting the source to the basin and where the canyon opens to the basin (Figure 2). This is the result of a reduction in transport capacity of the modelled river. This phenomenon has not been described in the Huesca DFS as the proximal zone close to the apex has been tectonically altered (Martínez Peña & Pocoví Juan, 1988). It has important implications for the reservoir quality predictions of the proximal zone in DFS as it was previously described as a single fully interconnected sandbody (Leleu & Hartley, 2010;Leleu et al., 2009;Nichols & Fisher, 2007).
The eastern and western transect also show that more, P but less CMS is deposited along these transects in comparison with the other ones (Table 3). In total, they have around 10% less net (P, CMS and FSS) than the other transects (Table 3). The strong decline in P is also T A B L E 2 Ratio of model output to outcrop measurements (Martin et al., 2021) for net to gross/in-channel component and sandbody/channel thickness F I G U R E 1 2 Stacked grain size distribution plotted as in Figure 6 but along transects from the apex to the eastern, southeastern, southern, southwestern and western corner of the model. Note that the distance from apex axes is not uniform as the transects are of different length further downstream starting at 10 km from the apex for the eastern and western transect and right at the apex for the other ones. The eastern and western transects are not mirror images though as 1% more P and 2% CMS are deposited in the western transect (Table 3). In summary, the greater the angle from the meridian (straight line from the source, through the apex to the most distal part of the model, see southern transect in Figure 2a), the coarser the average grain sizes are in the proximal zone but the finer in the medial and distal zones ( Figure 12). This directly translates to reservoir quality as the flow properties improve with grain size, for the same degree of sorting. The advantage of a numerical model is that the output can be sampled at every node/cell unlike outcrops, which are limited by exposure, which in turn is controlled by differential rates of erosion. This allows a quantification of variability, both laterally and vertically. As DFS are characterized as radial fluvial systems, the lateral variability encompasses all the locations with the same downstream distance from the apex. Figure 7 shows this variability for NTG, sandbody thickness, number of sandbodies and amalgamation ratio. The lowest variability occurs in the proximal zone as it is a single homogenous sandbody (approximately 0-15 km from the apex) but at the transition zone between the proximal and medial zone (approximately 15-25 km from the apex) the variability increases exponentially to its highest level found throughout the DFS model. Delineating the single sandbody part of the proximal zone from the increasingly heterogenous proximal to medial transition zone is crucial for reservoir quality predictions. The vertical connectivity (inferred from the amalgamation ratio) decreases from fully connected to 49% connected on average at 25 km from the apex. This rapid increase in reservoir heterogeneity is also supported by the peak in the number of sandbodies at 25 km from the apex. In the model, this variability is due to the varying intensity in fluid and sediment input at the source. The annual flood events are the most common but weakest in discharge, whereas the centennial flood events are the rarest but strongest in discharge. The stronger the discharge at the source the further the sediments will be transported before they are deposited. This can be seen in Figure 3 -250 kyr as the channelized flow deposits of sands (yellow) cutting through the flood plain deposits (green and grey). As the display interval includes multiple flood events, at this temporal resolution, the variability cannot be directly linked to single flood events. The variability in the calculated parameters decreases distally but only reaches zero in the lacustrine zone (87 km from the apex). This means that throughout the DFS and even on the margins around the DFS, variability is an important factor to quantify grain size distribution and ultimately reservoir performance. This lateral variability has been described by Hirst (1991) and Martin et al. (2021) in the proximal and medial zone of the Huesca DFS and by Owen et al. (2015), Owen et al. (2017) for the entire Saltwash DFS, southwestern USA, but their limited coverage of outcrops did not allow them to fully quantify this variability at a system scale.
Continental sedimentary deposits especially within DFS rarely contain datable material so that temporal correlation and accumulation rates are frequently difficult to obtain from outcrops. Numerical models, on the other hand, have full temporal and spatial control in the outputs. This work has revealed that the accumulation rate is slightly higher on average in the medial zone than in the proximal zone even through the average grain size is smaller. In the model settings, the angle of repose increases with grain size. This means that other factors lead to increased deposition in the medial zone such as more sediment bypassing the proximal zone and being deposited in the medial zone as the flow rate decreases enough to deposit most of the sediments it was transporting.

| CONCLUSIONS
As Best and Fielding (2019) said: 'The DFS debate demonstrates the essential need to link a knowledge of the processes of deposition with evolving models for their recognition, and that the two should not be divorced from each other'. We used the hydrodynamic stratigraphic forward modelling software Sedsim to reproduce the main characteristics of a DFS as constrained by outcrop measurements from the Huesca system. Linking the processes with the resulting deposits has led to three key findings, which expand our understanding of DFS. Even though the modelled DFS is laterally unconfined and displays a radial T A B L E 3 Grain size distribution along the five transects shown in Figure 12 and of the radially averaged representation shown in Figure 6 Transects Pebbles (%) pattern in planform, grain size distribution varies depending on the angle from the meridian of the DFS (straight line from the source through the apex to the distal zone of the DFS). The larger the angle, the coarser the average grain size in the proximal zone and the finer the average grain size in the medial and distal zones. In the proximal zone where the canyon connecting the source to the basin opens to the basin, finer-grained sediments are deposited. The variability of net to gross, sandbody thickness, number of sandbodies and amalgamation ratio is greatest at the transition between the proximal and medial zone of the DFS. Finally, the accumulation rate is highest in the medial zone of a DFS. This means that other processes like bypass of the proximal zone and decreasing transportation capacity of the river downstream from the apex have a measurable influence on the location of sediment deposition. This shows that the DFS debate demonstrates the essential need to link a knowledge of the processes of deposition with evolving models for their recognition and that the two should not be divorced from each other.
Creating and analysing the model brings us closer to understanding the link between controls and resulting deposits in a DFS. Although the model built within this study can help to improve our understanding of the Huesca system and DFS as a whole, it also presents significant opportunity for further study. SFM allow the user to modify individual input parameters to determine their impact on the subsequent stratigraphic architecture. Having a robust and realistic model, which is benchmarked against a real-world system is the starting point for further sensitivity studies, such as spatial and temporarily varying subsidence rates or changing the sediment supply through time. Now established, this model will be used as a base-case for such studies.