Geobody 5 records stable and low SAR on the overbanks (compare to geobody 4 overlap within their standard deviation; **Figure 10** and **Table 1**). As reported by De Ruig and Hubbard (2006), the uppermost UPF and BHF show a decrease in grain size and an increase in mud content, which suggests a reduction of turbidity-current height, density and/or competence. The turbidity currents in the channel became less efficient in spilling their diluted tops onto the overbanks. This process also explains the two characteristics of the BHU in the seismic-reflection data

(**Figure 2**): Whereas on the northern overbanks the section above and below the BHU is characterized by continuous, high-amplitude reflectors, the channel around the BHU is highly erosive – evidenced by truncated reflectors and the biostratigraphic results from Well W (**Figures 2**, **6**, **11**, late appearance of H. ampliaperta, see section Formation Boundaries for discussion).

The SARs from the channel and overbank model in the time interval from 19.6–19.0 Ma are directly linked to the process which has caused the formation of the BHU. Zweigel (1998) pointed out that on the Bavarian Shelf the SAR became negative (erosion) in her Sequence 4a (interval below BHU, **Figure 3**). Sequence stratigraphy (Jin et al., 1995; Zweigel, 1998) as well as foraminifera analysis (Wenger, 1987) indicate a major sea-level decrease in the eastern NAFB. Erosion on the shelf (Zweigel, 1998) points to a redeposition of shelf sediments into the deeper basin in the period from 19.6 – 19.0 Ma also evident by the foraminiferal fauna in the channel sediments in Hochburg1 (**Figure 1**) (Grunert et al., 2013). If from 19.6 – 19.0 Ma the channel was mainly fed by redeposited material from the Bavarian Shelf, sediment-supply rate from the hinterland must have dropped accordingly within this period to maintain the relatively constant and low SARs. Consistently, Kuhlemann (2000) described on the Egerian/Eggenburgian boundary a ∼45% decrease in sediment-supply rates from the Eastern Alps and a decrease of ∼35% in rates from the Central Alps.

Zweigel (1998) explained her findings on the shelf (sediment starvation, an angular unconformity, uplift in the northern part of the NAFB, decrease in sea level, deepening of the basin at the axis) with the visco-elastic relaxation of the European Plate that has caused an uplift of the forebulge, propagation of the bulge toward the orogenic front and a decrease in elevation in the hinterland. This would explain the decrease in sea level, the erosion and redeposition of older foreland sediments into the deeper basin through uplift and subareal exposure of parts of the shelf, and the decrease in sediment-supply rates and grain size related to a decrease in elevation in the Alpine hinterland. Zweigel's (1998) model is further in line with the ceasing of northward progradation of the Eastern Alpine front during the UPF (Hinsch, 2013). The uplift and progradation of the forebulge does not require a northward progradation of the thrust front as the stress relaxation of the lithospheric plate is dependent on temperature, rheology, and initial stress (Beaumont, 1981). The doubling in lag times in apatite fission track ages from the proximal archives of the paleo-Inn from 2–4 Ma in Aquitanian to 6–8 Ma in Burdigalian sediments (Kuhlemann et al., 2006) further support the idea of a relief reduction in the hinterland. We therefore interpret the BHU as the result of a period of non-sedimentation on the overbanks and, likely, in the channel, followed by a period when the channel was fed by remobilized shelf sediments. Grunert et al. (2013) pointed out that the channel in Hochburg1 (**Figure 1**) cuts into the UPF, deposited middle Eggenburgian strata but lower Eggenburgian strata is missing. Therefore, we suggest that the period of non-sedimentation extends over the entire early Eggenburgian.

After the mfs at 19 Ma, the SAR indicates a 8-fold increase compared to the geobody 5 (**Figure 10**, 4.5-fold increase compared to geobody 3). During this period – the beginning of the Upper Marine Molasse – channel sedimentation terminated and hemipelagic sedimentation (**Figures 2–4**) was established (middle Eggenburgian), which was replaced by clinoform progradation from the southern margin into the basin in the late Eggenburgian. These clinoforms indicate a sediment transport direction roughly from S to N (**Figures 2**, **7** and **Supplementary Figure S1**) (Grunert et al., 2013). It is likely that the clinoform progradation caused the tremendous increase in the SAR from 19 to 18.1 Ma (**Figure 10**).

During the period from 19.0 to 18.1 Ma, the basin shallowed from 1,000 – 1,500 m water-depth (Grunert et al., 2013) to water-depths <200 m (Grunert et al., 2012). The basin was transformed from an underfilled foreland basin (characterized by wedge-shaped deposits) to a filled one (sedimentary deposits tabular) (**Figure 2**), an evolution that occurred ∼11 Ma earlier in the Swiss part of the NAFB (Sinclair, 1997). The shallowing trend cannot be explained by eustatic sea-level fall (Haq et al., 1987) or a drastic change in subsidence patters in the basin (Genser et al., 2007). Instead, the basin was filled-up by sediments deposited by prograding clinoforms (**Figures 2**, **10**). The migration direction of the clinoforms (from S to N, **Figure 2** and **Supplementary Figure S1**) indicate a southern source. Detrital zircon U-Pb ages further suggest a change in provenance, when compared to the underlying UPF (Sharman et al., 2018). The provenance change and the predominantly northward directed sediment supply is consistent with redeposition of the Augenstein Formation on top of today's Northern Calcareous Alps (Frisch et al., 2001).

Two of the four expected external forcings are mirrored in the SARs of the Oligocene/Miocene Upper Austrian NAFB; the visco-elastic relaxation of the European Plate at 19.6 – 19 Ma and the unroofing and redeposition of the Augenstein Formation from 19.0 Ma onward. However, our data do not record the initially hypothesized doubling in SAR in connection with the external signals of the exhumation of the Lepontine Dome and the Tauern Window. Both signals were described within the proximal sedimentary archives in the Swiss NAFB (Schlunegger and Norton, 2015) and the south-east German NAFB (**Figure 4**; Ganss, 1977). In case of the external forcing of the Lepontine Dome, signal shredding may have occurred within the fluvial sedimenttransfer zone. Sediments traveled a distance of ∼500–600 km from the proximal archives into the deep-marine sink. The intermittent basins, the Swiss and German NAFB, created transient accommodation space along the sediment-routing system in which the excessive sediment load was stored (Kuhlemann and Kempf, 2002) and, hence, the tectonic signal is not manifested in the SARs of the marine transfer zone of the study area.

The increase in SAR of ∼65–100% (Ganss, 1977) described from the deep-marine section (Wenger, 1987) of the Chiemgau Fan (**Figure 4**) and related to the exhumation of the Tauern Window in the Eastern Alps is also not recorded in our SAR results. At present, the Chiemgau Fan is located

merely ∼70 km south-west of our research area (**Figure 4**). This implies that the Puchkirchen channel system has been very efficient in bypassing sediment. The channel system is characterized by low sinuosity (1.05 – 1.5; De Ruig and Hubbard, 2006) and the absence of avulsions (Hubbard et al., 2009) throughout the >7.4 Myrs of channel activity, possibly due to the lack of available space in the narrow (<20 km) foreland (Hubbard et al., 2009). These channel characteristics may increase bypass efficiency through the reduction of areas of sediment deposition, such as channelavulsion strata and lateral accretion at inner meander bends. The bypass efficiency of the system may be an additional reason for the non-recording of the signal from Lepontine-Dome exhumation.

After the reorganization of the NAFB at 19 Ma, resulting in the termination of channel sedimentation and the decrease of bypass efficiency of the sediment-routing system, the basin in the study area evolved from the marine transfer zone to the marine, terminal sediment sink. Hence, the sedimentary archive in the Upper Austrian NAFB was able to record the redeposition of the Augenstein Formation.

#### Geodynamic Implications

With the new stratigraphic results, we are able to further constrain the evolution of the basin. The NSU formed between 28.1 (Soliman, 2012) and 26.9 Ma (**Figure 14**). Masalimova et al. (2015) suggested a tectonically induced oversteepening of the northern basin slope due to a deepening of the basin axis. A temporal correlation exists with a northward advance of the Alpine thrust front (Kuhlemann and Kempf, 2002) as well as the final nappe stacking and crustal thickening in the Eastern Alps (Frisch et al., 2000). Both processes lead to an extra load onto the European Plate and could account for the basin deepening.

During the deposition of the Puchkirchen Group, sedimentation rates in the basin were mainly controlled by autogenic processes (**Figure 14**). However, minor tectonic activity is evident in the basin on the southern slope where Penninic and Helvetic thrust sheets overthrusted and partly incorporated Puchkirchen strata (Covault et al., 2009; Hinsch, 2013). The northward progradation of the Alpine front ceased in the UPF and the deformation stepped back into the orogenic wedge (Hinsch, 2013).

The northern NAFB was uplifted (Zweigel, 1998) at 19.6– 19.0 Ma, reflected in the change in sediment-transport direction (**Figures 2**, **4**, **14**). The Puchkirchen channel consistently funneled sediments to the east (De Ruig and Hubbard, 2006) indicating an eastward-directed slope in the basin. The northward migration of the clinoforms in the HFM (**Figure 2** and **Supplementary Figure S1**) indicate that either the deep eastern basin part was cut off from the Puchkirchen Trough or the eastern part was uplifted. In the upper HFM (19.0–18.1 Ma), the Upper Austrian NAFB turned from an underfilled basin to a filled one further indicating an absence of subsidence (**Figure 14**). In contrast, the Southern Alpine Foreland Basin (Venetian and Friuli Basin) of the Eastern Alps from late Burdigalian (∼17 Ma) onward experienced higher subsidence rates and South-Alpine thrusts emerged (Mellere et al., 2000). These findings might indicate that the former pro-wedge foreland basin, the NAFB, was turned into a retrowedge position, and the former retro-wedge basin, the Southern Alpine Foreland Basin, into a pro-wedge one. Taking this interpretation further, our and Zweigel's (1998) findings could be interpreted as the surface expression of the postulated northward retreat and tearing event of European lithospheric mantle underneath the Eastern Alps at 25–20 Ma (Handy et al., 2015) as part of an early Miocene switch in subduction polarity to present northward subduction of the Adriatic Plate (Lippitsch et al., 2003). However, the subduction polarity switch and the exact timing and geometry of slab tearing are subject of extensive ongoing research (Handy et al., 2015, ongoing research of the Priority Program "SPP - Mountain building processes in 4D" of the German Research Foundation). At the current stage, the model of a slab polarity switch fails

to explain the uplift of the forebulge in the German NAFB (Zweigel, 1998) and the redeposition of shelf strata into the deeper basin. Therefore, at present we prefer the model of visco-elastic relaxation of the lithosphere (Zweigel, 1998) to explain our results.

#### CONCLUSION

The results from the analysis of foraminifera, calcareous nannoplankton, geochemistry (TOC, S, CaCO3, δ 13Ccarb) and spatially averaged SARs provide new constraints for the evolution of the Upper Austrian NAFB. Deep-marine sedimentation in the basin started with the Zupfing Formation (NP24), the Puchkirchen Group formed during NP25 to NN2 (26.9 – 19.6 Ma), and the Hall Formation accumulated during NN2 – NN3 (19.6 – 18.1 Ma). Deposits of the Puchkirchen Group are diachronic between the central basin axis and the northern overbanks. The differentiation between Upper and Lower Puchkirchen Formation first postulated by Papp and Steininger (1975) and perpetuated in most following studies proved to be difficult to identify and facies dependent.

Both large-scale unconformities in the basin were dated. The NSU was established between 28.1 – 26.9 Ma, possibly due to a tectonically induced oversteepening of the basin caused by an extra load onto the European Plate either related to a thrust front progradation or the final nappe stacking period in the Eastern Alps. The BHU (19.6 – 19.0 Ma) resulted from a period of non-sedimentation in the basin, followed by a period of redeposition of sediments from the shelf in Bavaria (Grunert et al., 2013).

The results of our SAR calculation show changes in SARs which are possibly related to two of the four expected tectonic forcings in the Oligocene/Miocene Upper Austrian NAFB. We did not detect the hypothesized doubling of SAR caused by the exhumation of the Lepontine Dome (Central Alps) and the Tauern Window (Eastern Alps), which is present in proximal (Bavarian and Swiss) sediment archives. These signals were either 'shredded' in the foreland during long (500–600 km) eastward sediment transport or not recorded in the deepmarine setting due to efficient sediment bypass within the submarine channel system.

A change in the foraminiferal assemblages suggests sediment reworking on the shelf and stable SARs despite a sea-level fall from 19.6-19.0 Ma. This period of sediment remobilization is interpreted as the result of the visco-elastic relaxation of the European Plate. This relaxation induced uplift in the shelf area of the basin and a decrease in relief in the hinterland, leading to stable SARs by the enhanced input of reworked shelf sediment.

The maximum SARs from 19.0 Ma onward are interpreted to result from unroofing and redeposition of the Augenstein Formation. This interpretation is furthermore supported by a shift in sediment-transport direction from east to north via northward prograding clinoforms and a major shift in sediment provenance (Sharman et al., 2018). During this time, the study area was transformed from a deep-marine transfer zone to a final sediment sink, efficiently preserving the tectonic signal in the SARs.

In the Hall Formation (19.6 – 18.1 Ma), the Upper Austrian NAFB experienced major changes in basin geometry (from underfilled to filled), subsidence, drainage network, and the ceasing of deformation along the Alpine thrust front (Hinsch, 2013). These changes may point to a large-scale orogen-wide reorganization triggered by a change in slab polarity (Lippitsch et al., 2003; Handy et al., 2015) which turned the former pro-wedge (>19.0 Ma) Upper Austrian NAFB into a retrowedge setting.

# DATA AVAILABILITY STATEMENT

All datasets generated for this study are included in the article/**Supplementary Material**.

## AUTHOR CONTRIBUTIONS

JH and GF wrote the manuscript with major contributions from PG, GA, and AB. JH and GF drew the figures. GF carried out the foraminiferal analyses under the supervision of PG. GA carried out the calcareous nannoplankton analyses. JH undertook the interpretation of the seismic-reflection data and sediment-volume calculation. AB developed the project idea and provided funding.

# FUNDING

RAG Austria AG provided funding for this study.

# ACKNOWLEDGMENTS

We would like to thank the RAG Austria AG for providing the 3D seismic-reflection data set and access to the samples. Claudia Puschenjak (University of Graz) is thanked for Leco analysis. Prof. Joachimski (GeoZentrum Nordbayern, University of Nürnberg-Erlangen) coordinated the isotopic analysis. Schlumberger is acknowledged for the donation of the Petrel (∗ ) license. JH would like to thank Schlumberger's Support for their help and the Deutsche Bahn for unexpected extension of working time. We thank Mark Handy for comments on an earlier version of the manuscript. We also thank R. Sachsenhofer and J. Lorenzo-Trueba, for their constructive criticism and Associate Editor JF for carefully handling the manuscript.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2019. 00302/full#supplementary-material

## REFERENCES

feart-07-00302 November 26, 2019 Time: 18:16 # 23

Zweigel, J. (1998). Eustatic versus tectonic control on foreland basin fill. Contrib. Sedimentary Geol. 20:140.

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer RS declared a past co-authorship with one of the authors PG to the handling Editor.

Copyright © 2019 Hülscher, Fischer, Grunert, Auer and Bernhardt. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Grain-Size and Discharge Controls on Submarine-Fan Depositional Patterns From Forward Stratigraphic Models

#### Nicolas Hawie1 \*, Jacob A. Covault 2 and Zoltán Sylvester 2

*1 Beicip Franlab, Rueil-Malmaison, France, 2 Bureau of Economic Geology, Jackson School of Geosciences, The University of Texas at Austin, Austin, TX, United States*

#### Edited by:

*William Helland-Hansen, University of Bergen, Norway*

#### Reviewed by:

*Michal Janocko, Equinor ASA, Norway Cornel Olariu, Jackson School of Geosciences, The University of Texas at Austin, United States Ian Antony Kane, University of Manchester, United Kingdom*

> \*Correspondence: *Nicolas Hawie<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Sedimentology, Stratigraphy and Diagenesis, a section of the journal Frontiers in Earth Science*

Received: *15 August 2019* Accepted: *29 November 2019* Published: *13 December 2019*

#### Citation:

*Hawie N, Covault JA and Sylvester Z (2019) Grain-Size and Discharge Controls on Submarine-Fan Depositional Patterns From Forward Stratigraphic Models. Front. Earth Sci. 7:334. doi: 10.3389/feart.2019.00334* Submarine fans are important components of continental margins; they contain a stratigraphic record of environmental changes and host large accumulations of oil and gas. The grain size and volume of sediment supply to fans is thought to control the heterogeneity of deep-water deposits; predicting spatial variability of sandy and muddy deposits is an important applied challenge in the characterization of fans. Here, we use DionisosFlow stratigraphic-forward models to evaluate the sensitivity of submarine-fan deposition to a range of grain sizes, with corresponding diffusion coefficients ranging from 10 to 100 km2 /kyr for coarse sand to silt/clay, and discharges. In general, finer grains are transported farther in our models because they have larger diffusion coefficients. Coarser grains typical of a sand-rich fan tend to pile up and compensationally stack at the mouth of a proximal feeder channel. Increasing sediment-gravity-flow discharge resulted in a thicker depositional system; however, relatively coarse sediment piled up at the mouth of the feeder channel, which created a slope that promoted basinward sediment transport. Our modeling results can be applied to predict the overall geometry, stacking, and grain-size distribution of submarine fans. Improved understanding of grain-size and discharge controls also informs interpretation of the stratigraphic record of submarine fans. For example, outcrop observations of heterogeneity and compensational stacking of depocenters can be related to changing boundary conditions, namely changes in the caliber and overall supply of sediment delivery to deep-water basin margins.

Keywords: sediment supply, grain size, compensational stacking, forward stratigraphic modeling, deep-water depositional settings, submarine fan, diffusion coefficient, sediment transport

# INTRODUCTION

Submarine fans are deep-water depositional systems that received terrigenous sediment from canyons and channels that extend across continental margins (Piper and Normark, 2001). The deposits of fans host a relatively complete and readily dated record of environmental changes, including tectonic deformation and climate, in their upstream sedimentary systems compared to onshore records (Clift et al., 2000; Clift and Gaedicke, 2002; Romans et al., 2016). The deposits of fans are also globally important oil and gas reservoirs (Weimer and Pettingill, 2007). The reservoir properties of sandy architectural elements of submarine fans and their lateral continuity and vertical connectivity are important issues for petroleum geology (Piper and Normark, 2001). Piper and Normark (2001) suggested that the distribution of sandy architectural elements is primarily controlled by grain size and sediment supply, and the overall geometry of submarine fans is influenced by basin setting.

Forward stratigraphic modeling can be applied to predict the location and heterogeneity of depositional systems and petroleum reservoirs (Miller et al., 2008), as well as the depositional response to controlling factors (Piper and Normark, 2001). For example, in an exploration scenario with lowresolution seismic-reflection data (generally two-dimensional, 2-D, profiles with frequencies of 5–20 Hz; Normark et al., 1993; Prather et al., 2012) and no lithologic control from well penetrations, seismic-stratigraphic interpretation and structural restoration can be applied to create a paleotopographic surface for modeling the location, size, shape, and sub-seismic heterogeneity and stacking of deposits (Groenenberg et al., 2010; Aas et al., 2014; Deville et al., 2015; Hawie et al., 2017; Barabasch et al., 2019). Commonly used geostatistical methods in reservoir modeling use semivariograms, geometric parameters, and/or training images to reproduce spatial statistics from available seismic-reflection and well data (Pyrcz and Deutsch, 2014; Pyrcz et al., 2015). Recently, stratigraphic-forward modeling has been used to incorporate quantitative, process-based geologic information to constrain reservoir modeling. For example, Sacchi et al. (2016) used a depth-averaged flow in the 2-D horizontal plane, basin-scale stratigraphic-forward model called SimClast (Dalman and Weltje, 2008, 2012) to simulate fluviodeltaic stratigraphic evolution and create a 3-D probability distribution of facies proportions. This probability cube was used as additional input for a geostatistical reservoir model. Similarly, Falivene et al. (2014) improved DionisosFlow basinscale stratigraphic-forward model predictions of stratigraphic trap and reservoir presence by calibrating the models to independent constraints, such as thicknesses from seismicreflection or well data.

Hawie et al. (2018) used DionisosFlow at finer temporal and spatial resolution (i.e., 103 -104 yr time steps within an area of 17 × 17 km with cell sizes of 200 × 200 m) to simulate the stratigraphic evolution and sub-seismic heterogeneity of a Pleistocene submarine fan on the continental slope of tectonically active eastern Trinidad (**Figure 1**). A regional seismic-stratigraphic horizon was used as an initial paleotopographic surface input to the forward stratigraphic model. Over a range of time steps, compensational-stacking patterns governed the lateral continuity and vertical connectivity of sandy and muddy architectural elements of the submarine fan (Hawie et al., 2018); similar compensational stacking patterns in fans are common in other settings (e.g., Deptuck et al., 2008; Prelat et al., 2010; Jobe et al., 2017; Sweet et al., 2019). However, in all models of Hawie et al. (2018), thicknesses in the proximal areas of the models exceeded thicknesses observed in the field (**Figure 2**).

Here, we revise the reference-case forward stratigraphic model of Hawie et al. (2018) to achieve a better thickness match with the field example offshore Trinidad. Then, we evaluate the sensitivity of this model to input variables, namely: diffusion coefficients related to a range of grain sizes and sedimentgravity-flow discharges. We use an automated multi-simulation workflow using a Latin Hypercube Experimental Design (McKay et al., 1979) with thickness calibration to quantify the variance of thickness and sand distribution, which can be applied to de-risk petroleum-reservoir presence. Our forward stratigraphic models show the influence of key controlling factors, namely grain size, with corresponding diffusion coefficients ranging from 10 to 100 km2 /kyr for coarse sand to silt/clay, and sediment supply (Piper and Normark, 2001), on fan deposition. Understanding controls can be applied to the interpretation of fan stratigraphy in outcrops and subsurface datasets (e.g., Burgess et al., 2019).

# PREVIOUS WORK

Hawie et al. (2018) used seismic-reflection horizons from the tectonically active continental slope east of the Columbus basin and along the southern margin of the Barbados accretionary wedge (Brami et al., 2000; Wood and Mize-Spansky, 2009; the "NW" depocenter between the Darien and Haydn ridges in block 25A offshore of Trinidad) to define the top and base of a submarine fan. They simulated the evolution of this fan (124.5– 24.5 ka) with a series of DionisosFlow forward stratigraphic models within a domain of 17 × 17 km (200 × 200 m cell sizes), and assessed the impact of duration of time step (20, 10, 5, and 1 kyr) on sub-seismic stacking and heterogeneity of leveedchannel and lobe architectural elements (**Figure 1**). They used the regional seismic-reflection horizon at the base of the fan, which overlies mass-transport deposits mapped by Moscardelli et al. (2006), as the initial topography of the model. A single 300–600 m-wide feeder channel was located in the south of the model domain; it delivered grain sizes ranging from silt/clay to coarse sand. No differential subsidence was used in the model. Hawie et al. (2018) simulated the transport of relatively coarse grains: 20% coarse sand, 30% medium sand, 40% fine sand, and 10% silt/clay. The transport parameters used for the referencecase model ranged from 10 to 100 km2 /kyr for water-driven diffusion (or, in the case of submarine fans, sediment-gravityflow-driven diffusion; Kw) and 0.001–0.1 km2 /kyr for slopedriven transport (Ks).

In all of the models of Hawie et al. (2018), varying the simulated time step (20, 10, 5, and 1 kyr) resulted in a similar thickness trend and compensational stacking of depocenters. All models showed three to four major phases of sediment diversion during the migration of a relatively coarse depocenter. Moreover, the overall proximal-to-distal trend from relatively coarse leveedchannel to finer lobe deposits was similar in all models. However, the proximal and distal areas of all the models exhibited thickness differences, locally greater than two-fold, compared to the field (**Figure 2**). In light of this previous work, we are motivated to explore two questions: (1) How can we achieve a better thickness match between model and field cases? (2) Furthermore, how do variable diffusion coefficients related to a range of grain sizes and discharges influence submarine-fan depositional patterns (Piper and Normark, 2001)?

FIGURE 1 | Reference-case model (10 kyr time step) modified from Hawie et al. (2018). (A) Isochore map of the entire model. (B) Cross sections of the model. Left (L) and right (R) orientations in cross sections are left and right in map in part (A). Cross-section locations are indicated in part (A). (C) Isochore maps of depositional sequences within the model showing major phases of sediment diversion.

# METHODS

#### Forward Stratigraphic Modeling

DionisosFlow software is a 4-D process-based deterministic multi-lithology forward stratigraphic model that simulates basin filling (Granjeon, 1997, 2014; Granjeon and Joseph, 1999). A range of sedimentary processes are modeled including diffusive sediment transport, delta autoretreat, incision, largescale avulsion, and slope failure in response to tectonic, climate, and sea-level fluctuations during millennia and longer time scales (e.g., Pinheiro-Moreira, 2000; Rabineau et al., 2005; Alzaga-Ruiz et al., 2009; Gvirtzman et al., 2014; Harris et al., 2016; Hawie et al., 2017). Detailed fluid dynamics are not considered in this model; the goal is to simulate the large-scale (102 -103 m cell size) and long-term (103 -105 yr time steps) evolution of basin fill.

Sediment transport equations are used to simulate the transport of various classes of grain size (e.g., clay to sand) across a basin. This stratigraphic model combines (1) linear slope-driven diffusion (transport proportional to slope), referred to as hillslope creep, and (2) non-linear water- and slopedriven diffusion, referred to as water-discharge-driven transport (Willgoose et al., 1991; Tucker and Slingerland, 1994; Granjeon, 1997; Granjeon and Joseph, 1999; Deville et al., 2015):

$$Q\_s = -\left(K\_s / \overrightarrow{\nabla} \, h + \, K\_{\text{w}} Q\_{\text{w}}^m \text{S}^n\right) \tag{1}$$

where Qs is sediment discharge (km3 /Myr), Ks and Kw are the slope- and water-driven diffusion coefficients, respectively (km2 /kyr), Qw is water discharge (m3 /s), n and m are exponents

that affect sediment transport capacity with values between 1 and 2 (Tucker and Slingerland, 1994), S is the dimensionless local gradient of the basin, and h (m) is topographic elevation (Granjeon, 2014). Sedimentation and erosion rates are quantified by a mass balance equation in 3-D for each class of grain size (Euzen et al., 2004). In general, if slope is not sufficiently large to promote sediment transport downstream, deposition occurs. Large Kw and Qw values can promote sediment transport across lower slopes, as indicated in Equation (1). Sediment-gravity flows, commonly turbidity currents, are the primary agents of sediment transport, erosion, and deposition in submarine fans (Bouma et al., 1985). We liken the water-driven diffusion coefficient and the water discharge to a sediment-gravityflow-driven diffusion coefficient and gravity-flow discharge, respectively, which govern the rate of sediment transport though the system. That said, the diffusion equation (Equation 1) does not account for stratification of flow; both the velocity and concentration maxima of turbidity currents occur near the base of the flow (Sequeiros et al., 2010; Eggenhuisen and McCaffrey, 2012). Choosing Kw values for different grain sizes is a challenge because published values span orders of magnitude and they depend on many factors in addition to grain size, such as lithology, roundness, and discharge (Rivenaes, 1992; Falivene et al., 2014; Gvirtzman et al., 2014; Harris et al., 2016). Following Hawie et al. (2018), we used a range of Kw from 10 to 100 km2 /kyr for coarse sand to silt/clay. In a recent publication using DionisosFlow to simulate sedimentation across a larger region offshore Trinidad, Deville et al. (2015) used a range of Kw from 100 to 1,000 km2 /kyr. They used much larger Kw values in order to transport sediment across a much larger area (1,200 × 1,200 km); their individual cell sizes are nearly the size of our entire model domain (Deville et al., 2015). Ideally, we would tune diffusion coefficients to produce geometries and grain-size distribution similar to deposits of a field case (Harris et al., 2016). We created a reference-case model during an initial phase of model calibration. In lieu of grain-size information confirmed by well penetrations, we tuned the variables of the diffusion equation, Kw and Qw, to achieve a thickness trend that is similar to that published in Hawie et al. (2018) (**Figure 1**).

In addition to the assumptions presented above, there are limitations to the use of diffusion-based models to simulate sediment transport and deposition. For example, diffusion does not model particle-particle nor particle-fluid interactions (Equation 1). Moreover, cell size within the model domain is hundreds of meters by hundreds of meters and time steps are thousands of years; therefore, diffusion only provides an approximation of particle transport and depocenter location across the model domain during a time step. Stratigraphic forward model results should be compared to stratigraphic interpretations based on real geological and geophysical data; we compare the thickness of our model to that published in Hawie et al. (2018).

#### Experimental Designs

In order to achieve a better thickness match between model and the Trinidad field case of Hawie et al. (2018), we manually generated simulations of various grain sizes (i.e., sand- vs. silt/clay-rich) during 10 kyr time steps. We slightly modified some of the other input variables of Hawie et al. (2018), such as the discharges and diffusion coefficients (**Table 1**). Similar to Hawie et al. (2018), we used the regional seismic-reflection horizon at the base of the fan, which overlies mass-transport deposits mapped by Moscardelli et al. (2006), as the initial topography of the model. A single 300–600 m-wide feeder channel was located in the south of the model domain; it delivered grain sizes ranging from silt/clay to coarse sand. We extended the model domain 5 km to the south (17 km wide × 22 km long) to better match the feeder-channel location observed in the seismic-reflection data (**Figure 3**). Sediment enters the model domain from the south (**Figure 3A**). Once we improved the thickness match with the field case of Hawie et al. (2018), we tested two experimental designs: (1) 80–20% silt/clay proportion (22 simulations) and (2) ±20% gravity-flow discharge (22 simulations). We performed multiple automated simulations using a Latin Hypercube Experimental Design to quantify the standard deviations of thickness and sand-depositional patterns. Latin Hypercube Experimental Design samples variables from uniform distributions and ensures that the ensemble of samples is representative of the natural variability of the system, contrary to simple random sampling in Monte Carlo studies (McKay et al., 1979). In statistical sampling, a square grid containing sample positions is a Latin square if there is only one sample in each

#### TABLE 1 | Reference-case model inputs.

row and each column. A Latin Hypercube is the generalization of this concept to an arbitrary number of variables, whereby each sample is the only one in each axis-aligned hyperplane containing it. When sampling a function of N variables, the range of each variable is divided into M equally probable intervals. M sample points are then placed to satisfy the Latin Hypercube requirements; this forces the number of intervals, M, to be equal for each variable. The maximum number of combinations for a Latin Hypercube of N variables and M intervals can be computed with the following equation (Audze and Eglais, 1977; McKay et al., 1979; Iman et al., 1980, 1981):

$$\left(\prod\_{n=0}^{M-1} (M-n)\right)^{N-1} = \left(M!\right)^{N-1} \tag{2}$$

The following are features of Latin Hypercube Experimental Design: (1) it treats every variable as equally important and ensures uniformly distributed sampling; (2) it samples the full range of a variable; (3) it allows for multiple, multidimensional variables (e.g., 1-D sediment load and/or 3-D topographic variables); and (4) we, the designers, determine the number of simulations, which include output grids of properties for comparison to a reference case (Hawie et al., 2017).

Some caveats to our experimental designs include: (1) we restricted our model domain to the submarine fan, excluding the upstream canyon-channel system; and (2) we opened the northern boundary of the model and closed the eastern and western boundaries. Focusing only on the submarine fan required us to strongly vary diffusion coefficients in order to achieve a fan shape (e.g., Hawie et al., 2018). For example, more diffusion is required to produce the relatively flat, sheet-like geometries of the distal fan.

#### RESULTS

We improved the thickness match between our manually calibrated stratigraphic-forward model and the field case of Hawie et al. (2018) by modifying the grain-size inputs to be 5% coarse sand (10 km2 /kyr), 10% medium sand (15 km2 /kyr), 30% fine sand (30 km2 /kyr), and 55% silt/clay (100 km2 /kyr) (**Figure 3B**). This grain-size distribution is more similar to mudrich continental margins (e.g., Reading and Richards, 1994), such as offshore Trinidad, compared to the relatively sand-rich models of Hawie et al. (2018). Approximately 1/3 of the total sediment load was transported beyond the northern boundary of the model, and ∼80% of this sediment load was silt/clay. Although we improved the thickness match, especially in the proximal and central parts of the model domain, with model-thickness values >75% of the field case, there is still a mismatch in the western and thinner, peripheral parts of the model (∼30–60% of the field case). The overall architecture of the model is similar to the models of Hawie et al. (2018): compensationally stacked proximal leveed-channel depositional elements transition downstream to lobes (**Figures 3C**, **4**; **Supplementary Animations 1, 2**). The model output shows at least four phases of depocenter migration (**Figures 3C**, **4**; **Supplementary Animations 1, 2**). Initially, the depocenter was oriented southwest-to-northeast (0.124–0.104 Ma), then it shifted to the west and was oriented more south-to-north (0.104–0.094 Ma), followed by a more gradual shift to the east (0.094–0.064 Ma), where it split into three channels covering the model domain, and it was predominantly oriented southwest-to-northeast at the end of the simulation (0.054–0.024 Ma) (**Figure 3C**; **Supplementary Animations 1, 2**).

Following this initial phase of manual calibration, we implemented the two experimental designs. In the first experimental design, we varied silt/clay load, hereafter simply called mud load, 80–20%. **Figures 5**, **6** show the thicknesses and the distributions of coarse and medium sand in the maximum and minimum mud-load simulations (**Supplementary Animations 3**–**6**). Lower mud load results in relatively thick accumulations in the northeast distal part of the model domain, immediately downstream from the feeder channel; higher mud load results in a relatively thick band across the northern, central, and eastern regions of the model (**Figure 7A**). Both high- and low-mud-load cases resulted in similar patterns of shifting depocenters (**Figures 5**, **6**; **Supplementary Animations 3**–**6**). The map of thickness standard deviation of the 22 simulations shows larger variance in thickness (± nearly 20 m) near the proximal feeder channel (**Figure 7C**). The thickness standard deviation is also greater along channel forms extending to the northeast away from the feeder channel (**Figure 7C**). At first glance, the coarse and medium sand proportion maps look similar in both the high and low mud-load simulations (**Figures 5**, **6**). However, the map of sand proportion standard deviation shows large variance (± >10%) along the boundaries of the model domain, especially the western and southeastern boundaries (**Figure 7E**).

In the second experimental design, we varied the sedimentgravity-flow discharge ±20% relative to the manually calibrated model. **Figures 8**, **9** show the thicknesses and the

FIGURE 3 | New extended reference-case model. (A) Isochore map of the new model. (B) Difference in thickness of the model and the slope fan in block 25A, offshore Trinidad (Hawie et al., 2018). (C) Isochore maps of depositional sequences within the model showing major phases of sediment diversion.

distributions of coarse and medium sand in the maximum (+20%) and minimum (−20%) discharge simulations (**Supplementary Animations 7**–**10**). Discharge appears to play an important role in thickness difference, with nearly uniformly greater thickness corresponding with higher discharge (**Figure 10A**). In detail, higher discharge results in fewer major depocenters (at least two to three dominant depocenter- and channel-orientation directions) compared to lower discharge (similar to the manually calibrated model, at least four or five major depocenters) (**Figures 8**, **9**; **Supplementary Animations 7**–**10**). Moreover, the largest difference in thickness, with excess thickness in the case of higher discharge, is observed in the most proximal and distal parts of the model (**Figure 10A**). The map of thickness standard deviation of the 22 simulations shows larger variance (± nearly

20 m) in the central region of the model domain, extending from the proximal feeder channel basinward to the northern distal part (**Figure 10C**). As in the first experimental design, the coarse and medium sand proportion maps look similar in both the maximum and minimum discharge simulations (**Figures 8**, **9**). The map of sand proportion standard deviation of the 22 simulations shows larger variance (± nearly 10%) along the northwestern region of the model domain (**Figure 10E**).

#### Interpretations

#### Improved Model-Field Match

We improved the thickness match between our manually calibrated stratigraphic-forward model and the field case of Hawie et al. (2018) by modifying the grain-size inputs to be finer (**Figure 3**). In our modeling, sediment discharge is a function

of the diffusion coefficient Kw (Equation 1), which is larger for smaller grain sizes (**Table 1**). So, larger grain sizes assigned smaller Kw values will tend to pile up at the mouth of the feeder channel, like in the models of Hawie et al. (2018). By modifying the grain sizes in our new manually calibrated model to be finer, with correspondingly larger Kw values, the sedimenttransport equations of the model diffused sediment farther across the model domain and achieved a better match with the field. Similar to Hawie et al. (2018), the model shows repeated cycles of channel avulsion, compensational stacking, and unconfined deposition at the mouths of channels (e.g., Sun et al., 2010). In particular, compensational stacking is a key characteristic of submarine-fan deposits (e.g., Deptuck et al., 2008), and our results show that relatively simple diffusion-based models can produce realistic compensation patterns. In our model, topographic build up by deposition promotes compensation and depocenter migration around the model domain. However, the thickness of compensationally stacked depocenters in the proximal region of the model, near the feeder channel, is less pronounced than in the models of Hawie et al. (2018). This is because our new model comprises relatively fine grains with larger diffusion coefficients, which are more easily transported over topography and generate an overall smoother, more elongate fan geometry.

#### Experimental Design 1: Grain Size

The first experimental design highlights the influence of changes in diffusion coefficients related to a range of grain sizes on

fan deposition (**Figures 5**–**7**; **Supplementary Animations 3**–**6**). In the higher-mud-load case, a relatively thick band of deposits is spread across the northern, central, and eastern regions of the model as a result of the larger diffusion coefficient Kw (Equation 1) of finer sediment compared to the lower-mud-load case (**Figure 7A**). The larger diffusion coefficient promotes the basinward transport of finer sediment across the model. In the lower-mud-load case, coarser sediment piles up at the mouth of the feeder channel, which creates a steeper proximal slope that bypasses sediment to the distal northeast part of the model (**Figure 8A**). Moreover, the lower-mud-load case is relatively thick near the western and southeastern boundaries of the model domain, reflecting more compensational stacking (**Figure 6**). As diffusion is predominantly driven by slope and the diffusion coefficient, higher slopes and lower diffusion coefficients will promote more topographically directed, compensationally stacked deposits. The difference in depocenters in the high- vs. low-mud-load cases is also reflected in the maps of thickness and sand proportion standard deviations; the largest variance in thicknesses of the 22 simulations is at the mouth of the feeder channel in the proximal, central region of the model domain, and the largest variance in sand is at the western and southeastern boundaries of the model domain, where thicker, coarser deposits accumulated in the lower-mud-load case (**Figure 7**).

#### Experimental Design 2: Discharge

The second experimental design highlights the influence of sediment-gravity-flow discharge on fan deposition (**Figures 8**–**10**; **Supplementary Animations 7**–**10**). Intuitively, higher discharge (+20%) results in a thicker depositional system (**Figure 10A**). However, higher discharge results in fewer major depocenters because, initially, voluminous relatively coarse sediment accumulates in front of the proximal feeder channel and builds up a slope that dominates the direction of sediment transport during the simulation. The buildup of a slope in the proximal region of the model promotes bypass of finer sediment to the distal region, thereby causing some of the largest thickness differences between the high and low discharge cases along a

region extending from the proximal feeder channel basinward to the northern distal part of the model (**Figure 10A**). Moreover, the largest variance in thicknesses of the 22 simulations is along this north-to-south region, where thickness probably depends on whether high enough discharge promotes the development of a single major north-to-south depocenter, as in the higher discharge case (**Figure 10C**). The large variance in sand proportions of the 22 simulations is located in the northwest because sandy depocenters did not visit that region of the model domain in every simulation (**Figure 10E**). In contrast, the central region is immediately down slope of the proximal feeder channel and consistently received the bulk of the relatively coarse sediment load. Mud transport is less sensitive to topography and, as a result, can more easily spread uniformly across the model.

# DISCUSSION

The primary goals of our experiments were to understand the submarine-fan depositional response to changes in diffusion coefficients related to a range of grain sizes and discharges. In the models of Hawie et al. (2018), thicknesses in proximal areas exceeded thicknesses observed in the field. We increased the proportion of mud (i.e., silt/clay), with correspondingly larger Kw values, to the model by a factor of 5.5 in order to improve the thickness match between our reference-case model and the field case of Hawie et al. (2018). Rather than sand piling up in the proximal region of the model, at the mouth of the feeder channel, a mixed sediment load was more uniformly distributed across the model. Indeed, this is consistent with general models of submarine-fan run-out based on their mud- vs. sand-rich

sediment supply (e.g., Reading and Richards, 1994; Richards et al., 1998). Some of the largest, longest run-out submarine fans in the world are characterized as mud rich with highly "efficient" sediment transport (Richards et al., 1998; and references therein). Finer grain sizes are transported farther in our model; this has also been demonstrated for deposition on natural submarinefan systems, such as the Amazon fan (Pirmez and Imran, 2003), and in physical experiments (Baas et al., 2004). Simple models of sediment transport and deposition using an advection-settling scheme also predict this behavior (Straub et al., 2008; and references therein). For example:

$$\varkappa = U \frac{z\_i}{w\_s} \tag{3}$$

where the distance the particle travels (x) is equal to the velocity of the particle (U) multiplied by the of the initial height of the particle above the bed (zi) divided by the particle settling velocity (ws). Although this advection-settling scheme is different than how Equation (1) diffuses sediment across our model domain, it supports our experimental results that show, in general, keeping all other variables constant, a finer, slower settling grain will be transported a longer distance compared to a coarser, faster settling grain.

We also explored an experimental design of variable mud load (80–20%). Decreasing mud load (i.e., increasing sand load) resulted in more widely distributed, coarser depocenters across the model, from the western to southeastern boundary of the model domain (**Figure 7A**). This is consistent with conceptual models of sandy submarine fans based on outcrops

and high-resolution, shallow-subsurface data in tectonically active settings; they tend to have limited channel extension and more frequent avulsions associated with steep basinmargin slopes (Hoyal et al., 2014). These concepts have been further developed in physical experiments (e.g., Spinewine et al., 2009; Hamilton et al., 2015; Postma et al., 2016) and reduced-complexity modeling of submarine-fan evolution over an evolving topography (e.g., Burgess et al., 2019). Although our models do not capture the high-resolution detail of morphodynamic interaction between flow and topography as in some physical experiments, including a potentially important hydraulic jump flow transformation and resultant depositional architecture (Mutti and Normark, 1987), our models show sandier sediment supply associated with the general trend of more frequent avulsions resulting in a larger number of depocenters.

Increasing sediment-gravity-flow discharge had the intuitive result of overall thicker deposits. However, this also resulted in fewer depocenters because relatively coarse sediment piled up at the mouth of the feeder channel, which created a slope that promoted basinward sediment transport along a dominant south-to-north fairway. Piper and Normark (2001), in their analysis of the distribution of sandy and muddy architectural elements of submarine fans, interpreted that steepening of proximal channels promotes basinward sand bypass and the development of a pattern of sand distribution similar to detached lobes of Mutti (1979), with lobe deposition displaced from the limit of sandy channel deposition. Although this interpretation was related to depositional patterns of channels and lobes in response to avulsions on the Amazon fan (Pirmez and Flood, 1995), steeper channels appear to promote basinward sediment bypass in our models as well (**Figures 8**–**10**).

Our modeling results can be applied to predict the overall geometry, stacking, and grain-size distribution of submarinefan oil and gas reservoirs. In tectonically active, stepped-slope profiles, like offshore eastern Trinidad, increasing sand delivery to deep water (e.g., Moscardelli et al., 2012), can result in more compensational stacking of relatively sand-rich, proximal depocenters (Piper and Normark, 2001). Increasing sediment supply, but maintaining grain-size distribution, can result in basinward sediment transport along a dominant fairway, thereby potentially producing isolated, "detached lobe" depositional architectures downstream of primary feeder channels. Our modeling also informs the interpretation of the depositional record of submarine fans. Observing some of the aforementioned depositional patterns, such as more/less frequent avulsions and compensational stacking of depocenters (e.g., Straub and Pyles, 2012), can be interpreted in the context of changing boundary conditions, namely changes in the caliber and overall supply of sediment to deep water. In order to better understand the results of our modeling experiments, we maintain our input variables (e.g., silt/clay proportion and discharge) as constant for the entire simulation time of 100 kyr. Of course, in nature, boundary conditions change during periods as long as 100 kyr, and future work will follow some of the more recent forward modeling work of Sylvester et al. (2015) and Burgess et al. (2019), which vary inputs, including substrate mobility and discharge. We also aim to pursue more numerous simulations, of the order of thousands, in order to achieve better calibration with the subsurface; new developments in high-performance parallelized computing processes can accelerate computations of high-resolution DionisosFlow models by as many as five times (Granjeon et al., 2018).

#### CONCLUSION

We used DionisosFlow forward stratigraphic models to understand the submarine-fan depositional response to changes in diffusion coefficients related to a range of grain sizes and discharges. Finer-grained loads were transported farther across the model domain; coarser loads resulted in the active migration of the depocenter across the model domain. This is consistent with general models of submarine-fan run-out and compensational stacking based on their mud- vs. sandrich sediment supply. Higher discharge resulted in a thicker depositional system, but fewer channel avulsions because an initial relatively coarse-grained sediment buildup at the mouth of the proximal feeder channel promoted bypass of finer sediment to the distal region of the model domain. Our modeling results can be applied to predict the depositional architecture of submarine fans; changes in diffusion coefficients related to a range of grain sizes and discharges have an effect on the overall geometry, stacking, and heterogeneity of our models. Our results can also be applied to the interpretation of fan stratigraphy in outcrops and subsurface datasets. For example, observations of compensational stacking in outcrops can be related to changes in sediment supply to the depositional system. Our future work will evaluate the depositional response to temporally varying the inputs of more numerous simulations. We envision thousands of

# REFERENCES

Aas, T. E., Basani, R., Howell, J., and Hansen, E. (2014). Forward modelling as a method for predicting the distribution of deep-marine sands: an example simulations to generate facies probability maps to be integrated in reservoir models.

#### DATA AVAILABILITY STATEMENT

All datasets generated for this study are included in the article/**Supplementary Material**.

## AUTHOR CONTRIBUTIONS

NH, JC, and ZS designed the experiments and wrote the manuscript. NH performed the experiments. NH and JC created figures.

#### ACKNOWLEDGMENTS

We thank the sponsors of the Quantitative Clastics Laboratory (http://www.beg.utexas.edu/qcl) and Beicip-Franlab for access to DionisosFlow and CougarFlow forward stratigraphic modeling and multi-simulation software. We are grateful for thoughtprovoking comments from Tim Demko and Mauricio Perillo, as well as recommendations from reviewers Michal Janocko, Cornel Olariu, and Ian Kane and editors David Hodgson and William Helland-Hansen.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart. 2019.00334/full#supplementary-material

Supplementary Animation 1 | Reference case model of the coarse sand distribution (overhead view).

Supplementary Animation 2 | Reference case model of the coarse sand distribution (perspective view).

Supplementary Animation 3 | High mud load model of the coarse sand distribution (overhead view).

Supplementary Animation 4 | High mud load model of the coarse sand distribution (perspective view).

Supplementary Animation 5 | Low mud load model of the coarse sand distribution (overhead view).

Supplementary Animation 6 | Low mud load model of the coarse sand distribution (perspective view).

Supplementary Animation 7 | High discharge model of the coarse sand distribution (overhead view).

Supplementary Animation 8 | High discharge model of the coarse sand distribution (perspective view).

Supplementary Animation 9 | Low discharge model of the coarse sand distribution (overhead view).

Supplementary Animation 10 | Low discharge model of the coarse sand distribution (perspective view).

from the Peïra Cava Sub-basin. Geol. Soc. Lond. Spec. Publ. 387, 247–269. doi: 10.1144/SP387.9

Alzaga-Ruiz, H., Granjeon, D., Lopez, M., Seranne, M., and Roure, F. (2009). Gravitational collapse and Neogene sediment transfer across the western margin of the Gulf of Mexico: insights from numerical models. Tectonophysics 470, 21–41. doi: 10.1016/j.tecto.2008.06.017

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer CO declared a shared affiliation, with no collaboration, with one of the authors, JC and ZS, to the handling editor at time of review.

Copyright © 2019 Hawie, Covault and Sylvester. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Paleohydrology and Machine-Assisted Estimation of Paleogeomorphology of Fluvial Channels of the Lower Middle Pennsylvanian Allegheny Formation, Birch River, WV

#### Oluwasegun Abatan\* and Amy Weislogel\*

Department of Geology and Geography, West Virginia University, Morgantown, WV, United States

#### Edited by:

Brian W. Romans, Virginia Tech, United States

#### Reviewed by:

Adrian John Hartley, University of Aberdeen, United Kingdom Tara Nicole Jonell, The University of Queensland, Australia

#### \*Correspondence:

Oluwasegun Abatan<EMAIL_ADDRESS>Amy Weislogel<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Sedimentology, Stratigraphy and Diagenesis, a section of the journal Frontiers in Earth Science

Received: 15 August 2019 Accepted: 30 December 2019 Published: 22 January 2020

#### Citation:

Abatan O and Weislogel A (2020) Paleohydrology and Machine-Assisted Estimation of Paleogeomorphology of Fluvial Channels of the Lower Middle Pennsylvanian Allegheny Formation, Birch River, WV. Front. Earth Sci. 7:361. doi: 10.3389/feart.2019.00361 Rivers transport sediments in a source to sink system while responding to allogenic controls of the depositional system. Stacked fluvial sandstones of the Middle Pennsylvanian (Desmoinesian Stage, ∼310–306 Ma) Allegheny Formation (MPAF) exposed at Birch River, West Virginia exhibit change in sedimentary structure and depositional style, reflecting changes in allogenic behavior. Paleohydrologic and numerical analysis were used to quantify geomorphological and paleohydrologic variations reflected by MPAF fluvial deposits with the goal of understanding the controls on resulting fluvial sandstone architecture in these different systems. Channel body geometry, sedimentary structures, and sandstone grain size distribution were used to reconstruct the paleoslope and flow velocity of the MPAF fluvial systems. In order to enhance paleohydrological estimates, machine learning methods including multiple regression and support vector regression (SVR) algorithms were used to improve the dune height, and channel depth estimated from cross-set thickness. Results show that the channel depths of the lower MPAF beneath the Lower Kittanning coal beds tend to decrease upsection; this decrease is interpreted to reflect a transition from fluvial systems formed in a humid ever-wet climate to fluvial systems formed in less humid, seasonally wet, semi-arid climate. Paleohydrologic estimations enabled the evaluation of hydraulic changes in the fluvial depositional systems of the Appalachian Basin during the Desmoinesian stage. Paleoslope estimates indicated that the slope was low, which indicated that the fluvial gradient response was not driven by the effect of tectonic subsidence or uplift and sea-level change.

Keywords: fluvial, paleohydrology, sedimentology, allegheny, source-to-sink, machine-learning

# INTRODUCTION

Fluvial systems are the main terrestrial conduits for transporting the sediment load of a source to sink system. The source to sink system involves source rock erosion in the initial catchment area and sediment transportation through fluvial environments and to ultimate deposition in a basinal sink (e.g., Bhattacharya et al., 2016; Lin and Bhattacharya, 2017). The fluvial system

responds to external factors, such as climate, tectonics, and eustasy and is driven to maintain equilibrium while efficiently routing sediments. Changes in fluvial hydrology lead to changes in sediment transport and deposition, which alter channel aggradation, channel incision, and channel morphology (Leeder, 1993, 2009; Holbrook and Wanas, 2014). For example, fluvial channels in seasonal semi-arid climates have different geomorphology and hydrologic processes from fluvial channels of ever-wet humid climates (Fielding et al., 2009; Allen et al., 2014; Plink-Björklund, 2015). Modern fluvial depositional system analogs indicate that fluvial systems in seasonal semi-arid climatic regions typically have a greater channel width to depth ratio than fluvial channels of ever-wet humid climatic region (Fielding et al., 2009; Gibling et al., 2014).

Channel depth and width data combined with sedimentologic data from outcrop can be used to estimate paleohydrology for ancient fluvial systems (Rubin and McCulloch, 1980; Bhattacharya et al., 2016), which can then be tied to climate controls. Improved knowledge of the relationship between paleohydrology and depositional products of fluvial systems can also be used to improve reservoir characterization and reservoir quality prediction. In particular, the continuity and quality of fluvial sandstone reservoirs are dependent on the channel style of the fluvial depositional system (Miall, 1996; Bridge, 2009). Fluvial systems with high net-to-gross sandstone ratios form reservoirs with higher quality compared to fluvial systems with abundant overbank fine-grained sediments. Braided fluvial systems produce laterally continuous sandstone bodies with sheet geometries, while sinuous fluvial systems (meandering or anastomosing) produce laterally restricted sandstone bodies with ribbon and lens geometry (Miall, 1996). Braided channels have higher flow velocities because they are formed in areas with high slope, whereas sinuous fluvial system has relatively lower velocities because they are formed in areas with relatively lower slope (Schumm, 1981; Miall, 1996).

This paper proposes an enhanced methodology with which to estimate the paleohydrology and paleo-geomorphology of fluvial channels, using the fluvial sandstone deposits of the lower part of the Middle Pennsylvanian Allegheny Formation (MPAF) of Central West Virginia as a case example (**Figures 1**, **2**). The MPAF is characterized by repetitive cycles of clastic and chemical sediments known as cyclothems (Cecil, 1990). The MPAF at the Birch River area central West Virginia lacks marine zones where it is well exposed along a continuous road cut ˜110 m high and 500 m long along US 19 as it crosses Powell Mountain near Birch River in central West Virginia (**Figure 1**). The lower part of the MPAF (from here on referred to as MPAF) includes sandstones overlying the Lower Kittanning coal (LKC) beds, the Upper No. 5 Bock coal beds, and the No. 5 Block coal beds (**Figures 2**, **3**). Facies analysis determined channel style and geometry of the lower MPAF sandstones and revealed a range of channel forms, including high sinuosity, low sinuosity, and braided. This interval was selected for paleohydrological analysis because previous coal paleobotany studies indicate fluctuation between a humid and a seasonally wet-dry climate during MPAF deposition (Cecil, 1990; Eble, 2002; Cecil et al., 2003; Falcon-Lang, 2004; Greb et al., 2008; Falcon-Lang and Dimichele, 2010), and, thus, provides important independent constraints on paleoclimate variability with to investigate fluvial system response to paleohydrological controls.

Fluvial paleohydrology can be modeled from numerical equations based on grain size along with channel depth and width measurements and augmented by flow depth estimates from estimated dune bedform height (Ethridge and Schumm, 1977; Bridge and Tye, 2000; Leclair and Bridge, 2001; Leclair, 2002; Bhattacharya et al., 2016). These empirical equations relate sandstone grain size and channel geometry to estimates of paleohydrology. To build upon previous attempts at reconstructing paleohydrology of ancient fluvial systems, machine-assisted algorithms were developed to improve the accuracy of the estimated dune height from cross-set thickness using data of cross-set thickness and dune height from flume experiments reported by Leclair (2002) and Leclair and Bridge (2001). Multi-variate regression analysis was performed on the original data set to highlight the statistical significance (p-value) of the relationship between the variables in the data set. Support vector regression algorithm (herein and after referred to as SVR) was selected to better assess the relationship between variables with acceptable statistical significance (i.e., p-value < 0.05) because it can be used where bivariate relationships are established between geological properties with multivariate relationships (Ethridge and Schumm, 1977; Davis and Sampson, 1986; Bridge, 2009). Through this approach, the paleohydrological controls on MPAF fluvial architecture can be assessed to provide insights into the evolution of fluvial style and fluvial basin-fill record of the Alleghany foreland basin.

#### GEOLOGICAL SETTING

#### Geologic History

The MPAF (Desmoinesian Stage, ∼310–306 Ma) is part of an Upper Paleozoic cratonward prograding clastic wedge shed from the adjacent orogenic highlands of the Allegheny orogeny during the late Middle Pennsylvanian (Arkle et al., 1979; Donaldson and Shumaker, 1981; Ettensohn, 2008). The collision of Laurasia and Gondwanaland (∼325 Ma) initiated the Alleghenian orogeny, which was characterized by collision and compressional deformation structures that formed the Allegheny fold-thrust belt (Donaldson and Shumaker, 1981; Ettensohn, 2005, 2008; Sak et al., 2012). The Alleghenian orogeny resulted in the formation of a broad shallower foreland basin than the Acadian and Taconic orogeny (Ettensohn, 2005, 2008). Paleoclimate models developed using coal beds, paleosol, soil carbonate-based, and fossil leaf-based proxies indicate that paleoclimate shifted from ever-wet humid to seasonally arid conditions during the Middle Pennsylvanian (Cecil et al., 2003, 2004; Tabor and Poulsen, 2008; DiMichele et al., 2010; Falcon-Lang and Dimichele, 2010; Montañez et al., 2016). In particular, palynomorph studies of the MPAF showed tree ferns, which are common in less humid environments, increased and became more common in No. 5 Block and Upper No. 5 Block coal beds sections of the MPAF, whereas lycopsids, which are common in very humid environments, dominated the LKC bed (Kosanke and Cecil, 1996; Eble, 2002; Falcon-Lang and Dimichele, 2010).

The major driver of the paleoclimate change was attributed to the low paleo-latitudinal position of the Appalachian Basin during Middle Pennsylvanian; and the effect of glacial volumes at the poles on Hadley Cell circulation patterns along the Intertropical Convergence Zone (ITCZ) (Cecil and Dulong, 2003; Cecil et al., 2004). Changes to the Haley Cell circulation patterns along the ITCZ resulted in the seasonality of rainfall in low latitudes during glacial minimum and high rainfall during glacial maximum. The development of a rain shadow on the Alleghenian foreland basin, which is located on the downwind side of the orogenic highlands, may have also contributed to the drier climate (Tabor and Montanez, 2002; Tabor and Poulsen, 2008). Paleobotanical and sedimentologic studies indicate that earlier MPAF depositional systems formed in a humid climate, while the MPAF above the LKC beds were deposited in a semi-arid climate (Cecil, 1990; Cecil et al., 2003; Greb et al., 2008; DiMichele et al., 2010; DiMichele, 2013; Montañez et al., 2016). Paleogeographic reconstructions of the North American

craton suggest that the Appalachian basin was near the paleoequator with Appalachian highlands to the northeast and coastal lowlands located to the west (Archer and Greb, 1995; Cecil et al., 2004). The resulting paleo-gradient resulted in south and western drainage directions (Donaldson and Shumaker, 1981; Cecil et al., 2003, 2004). Paleodrainage models based on sedimentary analysis indicate the MPAF clastic wedge is composed of swamp, lacustrine, fluvial, and deltaic deposits (Donaldson and Shumaker, 1981; Cecil, 1990). Marine fossils observed in MPAF sandstones suggest the downdip extent of the fluvial segments of the MPAF prograding clastic wedge is located in southeast Ohio (Stubbs, 2018).

#### MPAF Channel Belts, Birch River, WV

The MPAF clastic units are subdivided based on coal beds, which stratigraphically oldest to youngest include: No. 5 Block, Upper No. 5 Block, Lower Kittanning (No. 6 Block Coal), and Middle Kittanning coal beds (Arkle et al., 1979; Blake et al., 2002; Eble, 2002). Palynomorph studies found that more lysosomes (fungi) spores, which are common in humid ever-wet environment are more abundant in early Middle Pennsylvanian deposits below the MPAF, whereas herbaceous fern plants, which are common in less humid environments, were more abundant in late Middle to early Upper Pennsylvanian deposits (Cecil et al., 1985; Kosanke and Cecil, 1996; Peppers, 1996; Eble, 2002). Sedimentologic models used lithologic climate indicators such as presence of caliche, calcareous pedogenic concretions, and siderite to assess climatic fluctuations during the Pennsylvanian, including parts of the upper MPAF which includes the Middle Kittanning, Upper Kittanning, Lower Freeport, and Upper Freeport coal beds and their associated clastic deposits (Donaldson and Shumaker, 1981; Cecil et al., 1985; Cecil, 1990; Cecil and Dulong, 2003).

# MATERIALS AND METHODS

#### Facies Architecture of the MPAF Channel Belts

Facies associations and architecture were used to interpret the fluvial styles (Miall, 1996; Bridge, 2009). Facies and facies association were identified from a 45-m thick and 495-m wide outcrop (Miall, 1996; Bridge, 2009). Data were measured using

a Wentworth calibrated grain-size card, measuring staff and ruler. Paleocurrent data were acquired from the left and right limbs of trough cross-strata using the best-fit circle method to determine paleocurrent direction on a stereographic plot (DeCelles et al., 1983).

#### Paleochannel Geometry Measurements and Estimation

Sedimentological data for the study was acquired from road cut (outcrop) along Route 19, Central West Virginia (**Figures 1**, **3**). Units of MPAF present at Birch River outcrop include the No. 5 Block coal bed, the LKC bed, shale and sandstone units above and below No. 5 Block, the Upper No. 5 Block, and LKC coal beds (Blake et al., 2002; Eble, 2002; Cecil et al., 2004; **Figures 1**, **2**). Sedimentologic data acquired from the outcrop include grain size, cross-bedding height, and barform height. These data were used to determine channel geometry (width and depth) by the methods outlined below.

#### Channel Depth

Channels were measured from preserved paleochannel boundaries corrected for compression during burial. Channel depths that were estimated from bar height involved the measurement of fully preserved channel bars in outcrop (Allen, 1970; Lin and Bhattacharya, 2017). The thickness of lateral or downstream accretion bars from outcrop, adjusted for 10% compaction factor, is representative of the bankfull channel depth (Ethridge and Schumm, 1977; Davidson and Hartley, 2010). The thickness of lateral or downstream accretion bars was determined using the fining upward sequence concept, where the lower and mid-section of the bar is characterized by planar, trough cross, planar cross, and inclined bedded sandstone that is relatively coarser than the upper section of the bar, which is characterized by massive and ripple bedded sandstone with plant debris and/or rooting structure (Bridge and Tye, 2000). The paleochannel flow depth was estimated from the thickness of lateral or downstream accretion macroforms using the equation by Ethridge and Schumm (1977):

$$d = D^\*/0.9\tag{1}$$

where D ∗ is the maximum channel depth, which is represented by the thickness of the sandstone macroform, 0.9 compensates for the compaction factor. Errors associated with this method can be up to 100% if it is used for muddy sections (Holbrook and Wanas, 2014).

Channel depths were also estimated from dune-scale crossset thickness, using empirical equations and machine-assisted algorithms. Previous work developed empirical equations which have been applied in the estimation of paleochannel dimension and morphology for ancient fluvial channel deposits. These equations determined relationships between the mean value of the exponential tail of the probability density function (PDF) for cross-set thicknesses and dune heights (Leclair and Bridge, 2001). The work by Leclair and Bridge (2001) has shown that the dune height (hm) can be estimated from mean cross-set thickness (Sm) using a regression equation:

$$h\_{\rm m} = 5.3\beta + 0.001\beta^2\tag{2}$$

$$
\beta = \text{ S}\_{\text{m}} / 1.8 \tag{3}
$$

which can be simplified as:

$$h\_{\rm m} = \text{ 2.9 S}\_{\rm m} \tag{4}$$

where hm is the mean dune height, Sm is the mean cross-set thickness, β is the mean value of the exponential tail of the PDF for topographic height relative datum. The range of error in this empirical equation is ∼20% (Leclair and Bridge, 2001). Hence, the authors suggest the equation be used on data set with similar standard deviation. Based on the observation that the ratio of bankfull depth to dune height is common between 6 and 10 (Bridge and Mackey, 1993; Bridge and Tye, 2000), dune height (hm) can be used to estimate channel bankfull flow depth (d):

$$6 \times d/h\_{\rm m} < 10\tag{5}$$

#### Machine-Assisted Estimation of Channel Belt

A support vector machine regression algorithm was developed to generate a new empirical equation that relates preserved crossset thickness to dune height to improve channel depth estimates from cross-set thickness. These relationships were established from measurements of dunes, and corresponding cross-set geometry produced under known hydrological conditions of a flume. First, multiple regression analysis using least squares elimination method was applied to the data set of Leclair (2002), which includes measurements of flow conditions and resulting bedform and cross-set heights, in order to determine the statistical relationship between dune heights and crossset thicknesses.

#### **Multiple regression analysis**

This study employs multiple regression to highlight the relationship between all the variables measured in flume studies that were used to explain the relationship between cross-set thickness and dune height (Leclair, 2002). The equation used for multiple regression is given as:

$$Y = \beta\_0 + \beta\_1 X\_1 + \beta\_2 X\_2 + \dots + \beta\_n X\_n + \text{ (6)}$$

where Y is the dependent variable represented as cross-set thickness, β0 is the intercept, βn are the coefficients, Xn are the independent variables describing flow conditions and depositional products, and ε is the random error. The accuracy of the multivariate regression was scored using R 2 . Then, SVR was applied to create and improve empirical relationships between the variables with the highest level of statistical significance as determined from the multiple regression analysis.

#### **SVR analysis**

Support vector regression is a type of supervised machine learning algorithm that fits as many instances in the model by taking into consideration the outliers in the dataset while developing an empirical relationship. The SVR machine

learning model was selected because it performs linear or nonlinear regression in a higher-dimensional space using linear, polynomial, or Gaussian kernels. The kernels transform the data into a higher-dimensional space by creating a vector from the evaluation of the test positions of all the data and establishes a linear, polynomial, or Gaussian relationship among the variables in the data. The Gaussian kernel uses normal distribution curves around data points to try to establish a relationship with the variables being considered. The advantage of SVR over linear regression is that SVR allows the model to be less fitted to the training data but more flexible for predicting new data (Zhang et al., 2014). SVR can also be used for multivariate regression; hence, new variables can be added to try to improve predictions. The simplified equation for predicting the dependent variable (Y) using the SVR model (Bao and Liu, 2006; Awad and Khanna, 2015) is given by:

$$Y = \; \boldsymbol{w}^T \boldsymbol{\phi}(\boldsymbol{x}) + \boldsymbol{b} \tag{7}$$

where Y is the dependent variable, w = (w0, w1, w2, . . .)T is the fitting coefficient in the higher dimensional space, φ is the kernel function transforming the independent variable x (crossset thickness in this paper) to a higher dimensional feature space, and b is the intercept. The model's performances compared to the previously used empirical equation (Eq. 2) were evaluated by testing the accuracy of the model's predictions using mean square errors (MSEs). Grid search algorithm was used to determine the best penalty parameter (C), fitting error (ξ ), and the kernel line of best fit for the data. Details of the algorithm and selected parameters to develop the SVR model are in the **Supplementary Material**. Algorithms for the SVR was written using Python and scikit-learn libraries (Pedregosa et al., 2011; Python, 2019). The steps taken to derive the SVR model for predicting dune height from cross set thickness include data preprocessing, kernel and parameter selection, and model fitting (**Figure 4**).

Data preprocessing. This includes sorting of the independent variable from lowest to highest value and normalizing the data. The cross-set thickness was set as the independent variable, while the dune height was set as the dependent variable for the SVR model. The cross-set thickness data were sorted and both data set were normalized. Normalization removes any disparity in the model that may be due to different units of measurements and large variance between values in the data that might skew the regression model in favor of the data set with larger values. The normalization of the independent and dependent variables involved adjusting both data set to a common scale. The normalization method used was MinMaxScaler, which has the ability to scale the data set between any range of values stipulated. The data were scaled into values between 0 and 1 using methods described in Pedregosa et al. (2011).

Kernel and parameter selection. Grid search was used to crosscheck all kernels and parameters until there was convergence, i.e., the ideal kernel and parameters that will give the best solution are determined. Kernel is a weighing factor between two sequences of linear and/or non-linear data, which enables the correlation of the data set in higher dimension space (Pedregosa et al., 2011). Three types of kernel considered are linear, polynomial, and

Gaussian kernels. The parameters considered include gamma, C, and epsilon (ξ ). The gamma parameter defines how far the influence of a single training example reaches and can be seen as the inverse of the radius of influence of samples selected by the model as support vectors (Pedregosa et al., 2011). The gamma range considered was from 0.5 to 0.8. The C parameter trades off correct classification of training examples against maximization of the decision function's margin; hence, the C parameter behaves as a regularization parameter in the SVM (Pedregosa et al., 2011). The range of C parameters considered was from 0.1 to 100. The epsilon defines a margin of tolerance where no penalty is given to errors (Pedregosa et al., 2011). The larger epsilon is, the larger errors you admit in your solution. The epsilon range considered was from 0.01 to 0.5. The gamma, C, and epsilon values, 0.8, 10, and 0.01, respectively, were selected because they produced the best SVR model. The accuracy of the SVR model, when using the Gaussian kernel and selected parameters 0.8, 10, and 0.01 for gamma, C, and epsilon values, respectively, is 84%.

SVR model fitting. Support vector regression model is fitted to the data set using the kernel and parameters from the grid search analysis. The model can be used to predict the dune height from cross-set thickness data inputted into the model. An inverse normalization is used to revert the normalized data and normalized model prediction.

#### Channel Width

Full channel widths were determined from Channel belt 1. Full channel belt widths could not be determined from the other channel belts because of the erosive nature of the channel boundaries. True channel width was derived from apparent channel widths measured from Channel belt 1 by correcting for

the orientation of the MPAF outcrop and paleoflow direction. Channel width was also estimated using published scaling relationships for channel geometry that takes into consideration the channel style as well as the tectonic and climatic setting of the fluvial systems (Gibling, 2006; Blum et al., 2013). The common range of channel width to depth scaling ratios selected from Gibling (2006) includes 5–50 for fixed river systems, which were used in Channel belt 1, 50–1000 for braided and lowsinuosity rivers used for Channel belts 2 and 3, and 30–250 for Channel belt 4.

#### Paleoslope

Paleoslope was estimated using grain size and density of sediment grains following the empirical equation of Holbrook and Wanas (2014):

$$\mathfrak{r}\_{\text{bf50}}^{\*} = \ (d\_{\text{m}}\text{S})/(\text{RD}\_{50}) = \text{constant} \tag{8}$$

where S is the slope, τ ∗ bf50 is the bankfull Shields number for dimensionless shear stress, dm is the mean bank full flow depth, R is the submerged dimensionless density of sand–gravel sediment (qs–qw), and D50 is the median grain size. τ ∗ bf50 is assumed to be 1.86 after Holbrook and Wanas (2014).

Grain size for this study was quantified from thin-section petrography, as well as estimated from observations of rocks in outcrop using a grain size card with graphical representation of Wentworth grain size classes. The error in grain size made from grain size cards has an error of about 1/2 phi (Lin and Bhattacharya, 2017). Four thin sections were selected that were representative of average flow in the four fluvial channel types interpreted in the lower MPAF. The thin-sections were acquired from above the scour deposits, which should be representative of deposits of moderate flow conditions. For each thin section, the maximum axis of at least 100 grains was measured and used to calculate median grain size (D50) for use in the empirical equation to estimate the paleoslope for MPAF paleochannels (Holbrook and Wanas, 2014).

#### Paleohydrology

Channel dimensions and paleoslope combined with flow velocity permit paleohydrologic reconstruction of MPAF channels using channel width derived from scaling factors so as to account for variabilities in channel cross-sectional area due to the depositional environment. Paleodischarge was estimated using the continuity equation (Eq. 9):

$$Q = \text{VA} \tag{9}$$

where Q is the instantaneous discharge and A is the crosssection area, which is the product of channel width and depth. Flow velocity was estimated by using sedimentary structures to infer the bedform for comparison with the bedform phase diagram of Rubin and McCulloch (1980). The dominant bedform observed in the channel belts was used to estimate flow velocity under the assumption that the dominant bedform reflects dominant bedload transport conditions during flooding events (Bhattacharya et al., 2016; Lin and Bhattacharya, 2017). The cross-sectional area was derived from estimated depth using the SVR machine-assisted model and width from width to depth scaling the relationship of modern and ancient fluvial channels (Gibling, 2006). We elected to not use empirical equations to estimate channel width from channel depth estimates, as this approach does not consider channel style as a variable in constraining channel width.

# RESULTS

#### MPAF Channel Belts

Facies and facies architectural analysis revealed nine lithofacies that represented fluvial channel deposits. The channel lithofacies include horizontally stratified sandstone, ripple-stratified sandstone, poorly sorted sandstone, planar cross-stratified sandstone, trough cross-stratified sandstone, massive sandstone, and low angle cross-bedded sandstone or convex upward sandstone, laminated mudrock, and massive mudrock facies (**Figure 5**). The coal beds overlie mudrocks and non-channel sandstones. We categorize the fluvial channel deposits into channel belts based on the channel planform.

#### Channel Belt 1: Low Sinuosity Fluvial System

Channel belt 1 is made up of multiple stories up to 10-m thick of tabular and lenticular, fine to medium-grained sandstones with sharp, sub-horizontal to horizontal, undulating erosional basal contact and sharp, curved erosional bounding surface above (**Figure 5**). Channel belt 1 overlies the No. 5 Block coal bed (**Figure 2**). Five stories were identified in Channel belt 1. Individual stories are characterized by multiple sandstone bed sets, which may be capped by mudrock, and are bounded above and below by an erosional surface. The two bottom stories are made up to 3-m thick lenticular sandstone separated by a sharp near-horizontal erosional surface. The lenticular sandstones comprise of convex upward, fine to medium-grained, massive (Sm), and trough cross-stratified (St) sandstone with sharp, curved bedding plane at the base. The massive and trough crossstratified sandstone beds are overlain by horizontal laminated sandstone (Sh) beds with a sharp horizontal bedding plane. The Sh is either onlaped by Sm or St beds with sharp, horizontal, or curved bedding planes. The Sm, St, and onlaped Sh beds are overlain by St and Sp beds with sharp horizontal bedding plane. The Sh may be overlain by interlaminated claystone, siltstone, and poorly sorted, ripple laminated sandstone in some places. The Sh may be overlain by interlaminated claystone, siltstone, and poorly sorted, ripple laminated sandstone in some places.

The three upper stories are made up of up to 2-m thick tabular sandstone bounded below by near-horizontal erosional surfaces. The tabular sandstones are made of tabular, fine to mediumgrained trough cross and planar cross-stratified (Sp) sandstones with sharp, horizontal bedding plane. The Sp overlies the St in the tabular sandstone. The Sp may be overlain by horizontal laminated sandstone beds in some places. The Sp beds are up to 1-m thick in some places. The uppermost tabular sandstone is erosionally truncated and overlain by sandstones from Channel belt 2. The sandstone of Channel belt 1 contains abundant coal and siderite intraclast, and fossilized plant fragments. Channel

belt 1 overlies the No. 5 Block coal bed. Paleocurrent data from trough cross-bedded sandstone indicate northeast to southwest direction of paleoflow. The lenticular sand bodies, which are overlain by Sh and onlaped by Sm and or St beds, are interpreted as mid-channel bar deposits while the tabular sand bodies are downstream accretion compound strata (Miall, 1996). The lack of a lateral accretion bar suggests low translation by channel. The abundance of coal and siderite intraclast suggest abundant vegetation and wet environment common in distal coastal plain depositional environments (Miall, 1996; Allen et al., 2014). Combined these features suggest Channel belt 1 are deposits of a distal, low sinuosity fluvial system.

#### Channel Belt 2: Braided Fluvial System

Channel belt 2 is made up of up to 5-m thick, multistory, amalgamated, medium-grained sandstone bounded above and below by sharp, undulating, horizontal, and curved erosional surfaces. Three stories were identified based on discontinuous, sub-horizontal, basal erosional surface, and channel lag deposits. Channel lag deposits, which comprise pebble size coal clast and iron-rich claystone clast and veins, were used to infer the base of the story where the basal erosional surfaces were not apparent. Individual stories are characterized by up to 0.3 m thick, amalgamated, medium-grained, compound through cross-stratified (St) sandstone beds with sharp or gradational, horizontal, or trough-shaped bedding plane. The St are rarely overlain by horizontal, fine-grained, ripple laminated sandstone (Sr) beds. Where the Sr is absent St may be overlain by up to 0.3-m thick, medium-grained, planar cross-stratified sandstone beds (Sp), or St. Channel belt 2 sandstones contain coal intraclast and petrified plant stems in places and is overlain by deltaic, lake, and well-drained floodplain deposits. The deltaic

deposits are characterized by coarsening upward, interlaminated shale, and very fine-grained sandstone, the lake deposits are characterized by laterally continuous, tabular, massive sandstone beds, while the well-drained floodplain deposits characterized by discontinuous, lens-shaped, coarsening upward ripple laminated sandstone beds and laterally continuous interlaminated siltstone, mudstone, and shale, which are overlain by the Upper No. 5 Block coal bed. Paleocurrent data from the trough cross-stratified beds indicate both northwest and southwest paleoflow direction. However, the dominant paleoflow is toward the northwest. Neither lateral nor downstream accretion macroforms were observed in Channel belt 2. The abundance of compound trough cross-bedded facies suggests the system was dominated by 3D dunes. The presence of bed sets bounded by curved and/or horizontal bedding planes and the absence of a clear macroform such as lateral or downstream accreting deposits suggest that Channel belt 2 is dominated by compound bars common in braided channel fills (Miall, 1996; Bridge, 2009; Allen et al., 2014). Combined all these features lead to the interpretation of Channel belt 2 as deposits of a braided fluvial channel.

#### Channel Belt 3: High-Velocity Channel

Channel belt 3 is characterized by a fine to medium-grained, single-story tabular sandstone body bounded above and below by erosional surfaces. The facies association of Channel belt 3 is made up of poorly sorted, planar cross-stratified, trough cross stratified, massive, horizontally stratified, low angle cross stratified, and convex upward sandstone strata (**Figure 5**). This channel belt is composed of two distinct sandstone units: A lower unit dominated by interbedded poorly sorted and ripple bedded sandstone that is 0.5–1-m thick, and an upper unit dominated by sandstones with upper flow regime structures (Allen, 1982) such as horizontal and low angle cross stratified and convex upward sandstone (Miall, 1996; Fielding et al., 2009). In some areas, the lower sandstone unit has some trough cross-bedded facies at the base, which transitions abruptly into the poorly sorted ripple laminated facies locally. The low angle cross-stratified and convex upward facies are the most dominant bedform in this channel belt. Channel belt 3 overlies the Upper No. 5 Block coal bed. The abundance deposits with upper flow regime structures and the abrupt transition in facies succession suggest deposits by a high velocity flooding event therefore Channel belt 3 deposits may be deposits of a high velocity channel. The presence of low angle cross-beds and convex upward strata suggest supercritical flow event (Bridge, 2009; Allen et al., 2014; Miall, 2014). The presence of very coarse horizontal, ripple, and poorly sorted bedded sandstone is indicative of a fluvial system with substantial erosive power.

#### Channel Belt 4: Sinuous Meandering Channel

Channel belt 4 is characterized by fine to coarse-grained, multistory inclined tabular and lenticular sandstone bodies bounded by sharp and erosional surfaces. A typical Channel belt 4 story is composed of poorly sorted sandstone, ripple stratified sandstone, planar cross-stratified sandstone, trough crossstratified sandstone, and massive sandstone facies (**Figure 5**). Each story has poorly sorted and massive sandstone beds overlying an erosional base. The poorly sorted and massive sandstone beds are overlain by inclined, trough cross, and planar cross strata, which may be draped by ripple stratified sandstone locally. The sandstone bodies of Channel belt 4 stack vertically and extend laterally to form lateral accretion bars. The lower bounding surface for this channel belt is undulating erosional. The upper bounding surface is covered by soil and vegetation in the study area. Three stories were identified in the study location. The deposits of the first story of Channel belt 4 are overlain by floodplain deposits characterized by laterally continuous carbonaceous shale and claystone. The floodplain shale and claystone are overlain by the LKC bed. The Channel belt 4 sandstones below the LKC are deformed and have abundant root traces. The second and third stories do not have coal beds and are not deformed. The inclined geometry, vertical and lateral succession of massive, trough cross-stratified, planar cross-stratified, and ripple laminated facies are lateral accretion (point bar) deposits, which are common in the sinuous meandering channel. This led to the interpretation of Channel belt 4 sandstone bodies as deposits of a sinuous meandering fluvial channel system. The deformed sinuous channels were interpreted as water escape features caused by the oversaturation of the sinuous channel deposits (Plink-Björklund, 2015).

#### Machine-Assisted Approach to Dune Height Estimation From Cross-Bed Height

#### Multiple Regression

The machine-assisted model was developed using multiple regression analysis and the SVR algorithm on the flume experiment data used by Leclair and Bridge (2001) to derive the empirical equations. A p-value of 0.005 was selected to test the significance of the statistical relationship (Davis and Sampson, 1986). Backward elimination showed that dune height (hm), as well as the dune length (lm) (i.e., the dune wavelength), had a high level of significance relationship (i.e., p-value < 0.005) with the cross-set thickness (Sm). It is difficult to measure the length of dunes in ancient deposits; hence, the relationship between cross-set thickness and dune height was further analyzed using SVR with the goal of developing a more efficient model for predicting dune height from the cross-set thickness in ancient channel deposits. These results from the multiple regression analysis highlighted the statistical relationship between cross-set thickness and dune height.

#### Support Vector Regression vs. Polynomial Regression

Dune height was set as the variable to be predicted from cross-set thickness (**Figure 6**) in order to compare SVR dune height predictions with the predictions from the empirical equation derived from polynomial regression (Eq. 4; Leclair and Bridge, 2001). The SVR model was used to predict dune heights from cross-set thicknesses derived from the flume experiment that Leclair and Bridge (2001) used to develop a polynomial regression model for dune height prediction from cross-set thickness. The Gaussian kernel, which used normal

distribution analysis on data points in order to highlight the bestfitted hyperplane in the dune heights and cross set thicknesses regression plot, produced the best predictions in the SVR model. The MSE and root MSE (RMSE) were used to compare the accuracy of predictions from both methods. The RMSE of predicted mean dune height from using the Leclair and Bridge (2001) model (Eq. 5) was 16.8 mm, whereas the RMSE of predicted mean dune height from the SVR model was 9.3 mm, indicating the SVR model estimates were closer to the actual dune height produced in the flume experiment.

#### Paleochannel Depth Estimates

Paleochannel depth was estimated from measured bar thickness corrected for compaction in Channel belts 1 and 4 using Eq. 1; however, only incomplete channel bars were observed in Channel belts 2 and 3. Measurements of these incomplete bars were recorded to constrain minimum paleochannel depth (**Table 1**). Paleochannel depths estimates were also determined from dune heights predicted from cross-set thicknesses for all of the MPAF channel belts using Eq. 5 (Bridge and Tye, 2000; **Tables 1, 2**). Paleochannel depth was estimated from the dune height (hm) based on the observation that the ratio of bankfull depth to dune height is commonly between 6 and 10 (Bridge and Mackey, 1993; Bridge and Tye, 2000). Paleochannel depth estimated from cross-set thickness was recorded for individual channel stories in Channel belts 1 and 4 (**Table 2**). Channel stories were not identified in Channel belts 2 and 3 because of the amalgamated nature of deposit from Channel belt 2 and the absence of channel stories in Channel belt 3. Overall, depth estimated from cross-set thickness measurements is greater than those estimated from bar thickness, suggesting that the bars of the MPAF have been largely subjected to substantial erosional truncation (**Figure 7**).

Bar thicknesses for Channel belt 1 were acquired from two channel bar deposits with a roll-over top, which is indicative of non-eroded channel bar deposits (Chamberlin and Hajek, 2015). The bar thickness was acquired from story 4 of Channel belt 1 using the same method as shown in **Figure 8**. The uncompacted thicknesses of these bars are 2.2 and 4.4 m. The bankfull paleochannel depths estimated from cross-set thicknesses using the SVR model were determined from measurements of 32 cross-set thicknesses. Depth estimates range from 6.1 to 13.9 m (**Table 2**).

Depths could not be estimated from measured bar thicknesses for Channel belt 2 because of the compound nature of bars, erosional truncation, and partial preservation of the bar deposits. However, these partial bars can constrain minimum paleochannel to >0.83–2.2 m. The bankfull paleochannel depth estimated from the mean of 45 cross-set thicknesses ranges from 3.4 to 5.6 m.

Depths could not be determined from bar thickness for Channel belt 3 because of the lack of macroform scale sand bodies, which are used to interpret fluvial channel bars. Twentyeight cross-set thicknesses were used to calculate the mean cross-set thicknesses used for SVR dune height prediction and paleochannel depth estimation. This yielded a depth estimate of 2.8–4.7 m.

Depth was estimated from bars with roll-over tops, which indicate they are fully preserved fluvial channel bar deposits (**Figure 8**; Chamberlin and Hajek, 2015). Bar thicknesses, which were measured from two bars in story 1 are 4.6 and 5.5 m, and one bar in story 2 of Channel belt 4 is 3.3 m. The bankfull paleochannel depth estimated from 37 cross-set thickness measurements from Channel belt 4 ranges from 3.6 to 9.6 m. The similarity in SVR predicted and bar thickness predicted bankfull paleochannel depths from stories 1 and 2 of Channel belt 4 indicates that Paleochannel depths predicted using SVR are accurate (**Tables 1**, **2**).

#### Paleochannel Width Estimates

The apparent width of the channel boundaries were measured from the boundaries of incised channel deposits. The incised paleochannel widths represent minimum channel widths as the

#### TABLE 1 | Results of estimated paleochannel geometry.

feart-07-00361 January 14, 2020 Time: 15:32 # 11

TABLE 2 | Results of estimated paleogeometry (channel depth and width), paleoslope, and paleodischarge.

NA, not applicable due to insufficient data from outcrop; dm, mean bankfull depth. \*\*\*Estimated using bed-form diagram. \*\*Calculated using data from thin section and grainsize card. \*Data from grain size card.

upper channel margins were commonly truncated. Channel widths were determined from nine sand bodies in the channel belts (**Table 1**). Apparent channel width was corrected for true width using paleocurrent direction. Paleocurrent direction was determined by 116 measurements of the orientation and dip of trough cross-beds. Measured channel width estimates range

FIGURE 8 | Example of how direct measurement of preserved bar is measured from outcrop data in Channel belt 4.

from values greater than the 13.4 to 40.3 m measured widths from the outcrop.

Channel widths were also estimated from the scaling relationships of channel width (w) to depth (dm) ratio defined using modern systems by Gibling (2006). Estimated SVR channel depths were used to determine channel width. Channel width was estimated for individual stories in Channel belts 1 and 4 based on their interpreted channel styles. The paleochannel widths acquired from Channel belt 1 outcrop were from preserved channel boundaries, so the width was recorded as actual width. The paleochannel widths, which were determined from the measurement of five channel boundaries from Channel belt 1 deposits, are up to 40.3 m (**Table 1**). The measured paleochannel widths determined from the deposits of Channel belt 1 fall within the channel width ranges determined from the scaling

relationship. The width range determined from the scaling relationship is 31–695 m (**Table 2**). The width was determined by using the minimum and maximum w/d of 5 and 50, and SVR bankfull flow depth estimates of 6.1–13.9 m (**Tables 1**, **2**). The common w/d ratios for fixed channels in distal humid environments (Gibling, 2006) were used for widths analysis in Channel belt 1.

The measured paleochannel widths of the fluvial channel incision from Channel belt 2 are 15.2 and 28.28 m (**Table 1**). Two instances of channel incision were observed in Channel belt 2. The deposits of the incised channel were eroded therefore the upper boundaries of the paleochannels could not be determined. The w/d scaling relationship of 50–1000 for braided fluvial channels (Gibling, 2006) and mean flow depths estimated from the SVR were used for estimating widths for Channel belt 2. The paleochannel width range derived from w/d scaling relationship is 168–5600 m (**Table 2**).

Channel belt 3 does not have enough sedimentary features to determine channel width. It was difficult to measure any form of channel boundaries because deposits of Channel belt 3 did not have macroforms such as channel bars that could be used in identifying channel boundaries. The w/d scaling relationship of 50–1000, for braided and low sinuosity fluvial channels (Gibling, 2006) and mean flow depths estimated from the SVR were used for width estimate in channel belt 3. The paleochannel width range derived from w/d scaling relationship is 141–4700 m (**Table 2**).

The paleochannel incision widths measured from channel belt 4 are up to 27.6 m (**Table 1**). There were three instances of channel incision in Channel belt 4. The deposits in the channel incision were also eroded at the upper section, which made it impossible to estimate the true channel widths. The w/d scaling relationship of 30–250, for meandering fluvial channels (Gibling, 2006) were used for estimating widths in Channel belt 4. The paleochannel width range derived from w/d scaling relationship is 108–2400 m (**Table 2**).

#### Paleoslope Estimation

The paleoslope was estimated using grain size determined from thin-section petrography and a grain size card with a graphical representation of Wentworth grain size classes (grain size table in Channel Paleohydrology **Supplementary Material**). The overall grain sizes of the MPAF channels deposits range from pebble to clay sizes, which is common in fluvial channel deposits (**Table 2**). The clastic sediment of the MPAF fluvial channels are moderately sorted. The sediment grain sizes of Channel belt 1 ranges from coarse to fine-grained sand (0.5–0.17 mm) with the fine to medium-grained sand (∼0.25 mm) being the most dominant mode. The grain sizes of Channel belt 1 from thin section analysis yielded a D50 grain size value of 0.23 mm, categorized as finegrained sandstone. The sediment grain sizes of Channel belt 2 vary from coarse to fine-grained sand (0.1–1.05 mm). Thin section analysis of Channel belt 2 resulted in a D50 grain size value of 0.33 mm, categorized as medium-grained sandstone. Channel belt 3 has sediment that varies from pebble to mediumgrained sand (>2–0.25 mm). The D50 value of Channel belt 3 is 0.3 mm, categorized as medium-grained sandstone. Channel belt 4 sediment clast size varies from pebble to mud (>2 to <0.088 mm). Channel belt 4 grains have a D50 value of 0.4 mm, categorized as medium-grained sandstone.

Paleoslope of the MPAF channels estimated using Eq. 8 ranges from 0.00007 to 0.0004 (**Table 2**), which suggests a low paleoslope comparable to slope ranges for the Amazon, Mississippi and Niger Rivers (slope range ∼0.00002–0.0005; Blum et al., 2013). The estimated paleoslope for Channel belt 1 0.00007–0.0001. The estimated paleoslope for the other channels are: Channel belt 2 are 0.0002–0.0003, Channel belt 3 are 0.0002–0.0003, and Channel belt 4 are 0.0001–0.0004, which are an order of magnitude steeper than the estimated slope of the low sinuosity channel. The estimated lowest slope for Channel belt 1 agrees with the dominant fine grain size observed in the sand body.

#### Paleohydrology

The MPAF channel belt flow velocities, which were estimated using the bedform phase diagram, ranges from 0.6 to 2 m/s (**Figure 9** and **Table 2**). The flow velocity of Channel belt 1 is in the range of 0.85–1.5 m/s. The estimated flow velocity of Channel belt 1 was determined by using the bankfull depth range of 6.1–13.9 m, the dominant sedimentary structure, which is trough cross-stratification produced by dune bedforms, and the fine-grained sand bed form diagram. The estimated flow velocity of Channel belt 2 ranges from 0.6 to 1.6 m/s. The flow velocity of Channel belt 2 was determined using the bankfull depth range of 3.4–5.6 m, and plotting the chart area for the dominant sedimentary structure, which is dune scale cross-sets, on the medium-grained sand bedform phase diagram. Channel belt 3 deposits are dominated by lamination produced from lowamplitude, upper flow regime bedform; therefore, the bankfull depth range of 2.8–4.7 m was used to estimate a velocity range of 1.5–2 m/s. The velocity of Channel belt 4 was determined using the dominant sedimentary structure, which is the dunescale cross-sets and the estimated channel depth range of 3.6– 9.6 m to estimate the velocity, which ranges from 0.6000 to 1.7 m/s. Paleodischarge for the channel belts range from 206 to 38658 m3 /s. Channel belts 2 and 3 with paleodischarge values of 621–34,496 and 695–38,658 m3 /s are the highest paleodischarge. The other paleodischarge ranges are 206–11351 m3 /s for Channel belt 1 and 408 to 27648 m3 /s for Channel belt 4 (**Table 2**).

# Errors and Uncertainties Associated With Numerical Analysis

Measurement of channel fill structures in the outcrop is subject to bias in interpretation sedimentary features from outcrop data, which represents an initial source of error. Detailed architectural analysis of outcrop aided the identification of channel bar (Bridge, 2009; Holbrook and Wanas, 2014). Other problems associated with measuring thickness data from channel bar include compaction and erosion. The degree of compaction is dependent on several factors such as original packing, original void ratio, shape of grains, degree of roundness of grains sand composition, and size grading (Ethridge and Schumm, 1977). Equation 1, which takes into account factors that affect the degree of compaction, was used to compensate for the 10% error due

in measured bar thickness to compaction. It is impossible to estimate channel depth from an eroded bar; hence, data from eroded bar were recorded to give a minimum estimate of channel bar thickness. Errors associated with identifying and measuring channel bars can be up to 60% (Holbrook and Wanas, 2014). Errors identifying channel bar thickness using story thicknesses can be up to 25% (Holbrook and Wanas, 2014). Mean channel depth estimated from cross-set thickness is subject to the bias of preferential preservation of dunes during waning flow, which may not be representative of actual bankfull flow. This results in the estimation of mean dune heights that are not representative of actual flow conditions, which cause errors in the depth estimates. Mean channel depth estimated from dune height and cross-set thickness by using Eq. 5 have errors up to 25% if the depth range is averaged (Holbrook and Wanas, 2014). The accuracy of dune heights estimated using SVR is 84%. The score was calculated using the coefficient of R 2 . Errors associated with estimating channel depth with Eq. 5 will also be encountered (Pedregosa et al., 2011). Measurement of channel width is also subject to the bias of interpretation, which may lead to errors in data acquired. Accuracy of channel width measurement is dependent on the identification of channel banks, which proved to be extremely difficult in the study outcrop. Channel width estimated form scaling factors have error ranges by a factor of ±4 (Blum et al., 2013). Channel dept and width estimated from the rock record are representative of extreme events, which may have resulted in extreme geomorphology and discharge than is normal to the depositional system (Gibling, 2006). Error associated with slope estimates is the assumption that bed shear stress required to move the sediment load is constant across the channel. The bankfull shield number used as a constant is for slope estimation varies by ±2 (Holbrook and Wanas, 2014). The grain size used for estimating slope is an average value and was measured with a Wentworth grain size calibrated card, which has errors of ∼1/2 phi (Lin and Bhattacharya, 2017). The errors associated with paleodischarge estimates are up to an order of magnitude (Holbrook and Wanas, 2014; Lin and Bhattacharya, 2017). The instantaneous paleodischarge estimates

are not representative of the annual or seasonal discharge in the paleochannel. The instantaneous discharge equation is a function of velocity of bankfull floodwaters, cross-sectional area of a channel, assuming the sediment supply is constant. The assumption of constant sediment supply adds more errors to the paleodischarge estimates as sediment type and sediment load varies. Another error associated with the paleodischarge estimate is based on the assumption that the cross-sectional area of a channel is the same, which is not so (Holbrook and Wanas, 2014). The cross-sectional area was also calculated using depth and width estimates, which means the errors from those numerical estimates are reflected in paleodischarge estimates. Additionally, the estimates from the instantaneous discharge have errors because the cross-sectional area for the paleochannels is estimated from mostly eroded channel deposits, which may not be representative of the actual sediment load.

## DISCUSSION

#### Machine-Assisted Paleohydrological Analysis

The machine-assisted SVR analysis increased the accuracy of dune height prediction from cross-set thickness resulting in a higher accuracy of estimated channel depth. The results of the SVR analysis were compared to the widely used polynomial regression model developed by Leclair and Bridge (2001) using MSE and overall predictions showed comparable and in most cases better performance. The result from MSE analysis showed that the SVR model had a RMSE of 9 mm while the polynomial model had an RMSE of 16.8 mm.

#### Advantages of SVR Over Polynomial Regression

The SVR model considers outlier data when estimating dune cross-strata thickness. Furthermore, SVR analysis is done in a higher dimension so it allows for a comparison of additional variables that may improve prediction. For example, cross-set length may be added to the SVR model to improve the accuracy of predicted dune height. In geology, all data are important because they reflect the variability of the conditions during bed formation. The model of Leclair and Bridge (2001) was based on the data along the best line of fit. The predictions of the SVR model are non-linear due to the use of the Gaussian kernel, i.e., the model uses normal distribution to predict the zone of best fit. The result of utilizing outlier data and the Gaussian kernel is a more accurate prediction as evident by the RMSE of SVR compared to the previous model.

#### Channel Belt Evolution

Facies architecture, channel dimensions, paleoslope, and paleodischarge of the MPAF channel belts reveal an evolution of fluvial channel form in response to changing paleohydrological conditions. Of note, the channel depth shows variability among MPAF channels: 6.1–13.9 m for Channel belt 1, 3.4–5.6 m for Channel belt 2, 2.8–4.7 m for Channel belt 3, and 3.6–9.6 m for Channel belt 4. Channel width estimated from scaling relationships showed that the low sinuosity channels had a lower range of width, while the braided channels had the largest width range. The independently estimated paleoslope and paleodischarge, which were compared among the MPAF channel belts (**Table 2**), showed that all the channel belts at the study location had a low slope (0.00007–0.0004) and variable paleodischarge (**Figure 10**). The estimated paleoslope of the lower MPAF channel belts is similar to slope ranges for the Amazon, Mississippi, and Niger Rivers (slope range ∼0.00002– 0.0005; Blum et al., 2013). This indicates that the paleoslope estimates obtained for the MPAF fluvial systems are consistent with physiographic models that suggest MPAF were deposited in low-gradient environments caused by unloading type relaxation of the foreland basin (Cecil, 1990; Cecil et al., 2003; Greb et al., 2008). Given the relatively low paleoslope estimated for all MPAF channel belt, the formation of upper plane beds and low amplitude bedforms from high-velocity flow is considered to reflect flooding events caused by intense precipitation (Kosanke and Cecil, 1996; Miall, 1996, 2014; Cecil and Dulong, 2003; Cecil et al., 2003).

Results indicate Channel belt 1 deposits formed from a low gradient, fine-grained, low-sinuosity channel form in an anastomosing fluvial system. Channel belt 1 has the thickest channel depth (6.1–13.9 m), which results in a relatively low w/d. This low w/d combined with the grain size, slope, and an abundance of coal intraclast and large plant fragments may reflect channel confinement due to bank stabilization by vegetation. The abundance of coal intraclasts and large plant fragments suggests abundant vegetation, which flourished in the humid climate, surrounded areas in the low sinuosity fluvial system (Cecil, 1990; Cecil et al., 2003; Allen et al., 2014). The estimated range for the paleoslope of the low-sinuosity channel of 0.00007–0.0001 is an order of magnitude lower than other MPAF channels in the study area. The low slope and high amount of fine-grained deposits may have contributed to Chanel belt 1 having the relatively highest channel depth and a low w/d (Gibling, 2006). The estimated paleodischarge for Channel belt 1 is 206–11351 m3 /s. This range of paleodischarge rate is the lowest of all the MPAF channel belts in the study area and is consistent with the abundance of fine-grained sediments deposited (**Figure 10**; Miall, 1996, 2014; Catuneanu, 2006).

Channel belt 2 directly overlies Channel belt 1 (**Figures 2**, **3**) and contains features that indicate deposition by a braided fluvial system in a humid climate, with low channel confinement, low paleoslope, and high paleodischarge. Channel belt 2 deposits are characterized by an abundance of trough cross-stratified, medium-grained sand, with a paleochannel depth range of 3.4– 5.6 m. The abundance of trough cross-stratification suggests that the channel system was dominated by 3D dunes. The low paleochannel depth, higher slope values (0.0002–0.0003) and higher paleodischarge (621–34496 m3 /s) compared to other MPAF channel belts (**Figure 10**), suggests that the braided channel style of Channel belt 2 is due to mainly to an increased slope gradient. Channel belt 2 is interbedded with some plant trunk and coal clast, which suggests an abundance of vegetative material in the fluvial system; this suggests that the area was vegetated, but this vegetation wasn't sufficient for bank stabilization.

Overlying the No. 5 block coal bed (**Figures 2**, **3**), Channel belt 3 sandstones were deposited by a low gradient, high-velocity fluvial system in a seasonal semi-arid climate. Channel belt 3 is characterized by an abundance of medium-grained, low amplitude bedforms and a relatively high width to depth ratio. The abundance of low amplitude bedforms in Channel belt 3 was attributed to high-velocity flooding event(s), which is supported by high paleodischarge values (695–38,658 m3 /s) compared to other MPAF channel belts. Channel belt 3 exhibits the shallowest paleochannel depths estimates compared to the other MPAF channels. This combined with the low paleoslope for Channel belt 3 (0.0002–0.0003) suggests that high-velocity flow resulted from a control other than paleoslope (Cecil et al., 2003; Cecil and Dulong, 2003; Plink-Björklund, 2015). The Burdekin River is an analog system developed in a seasonally wet-dry climate that experiences high velocity flows due to monsoonal precipitation events during wet seasons (Fielding and Alexander, 1996). The Burdekin River, just like CB 3, is dominated by upper plane beds and dunes, which represent fluctuating periods of extreme and moderate flow events.

Channel belt 4 directly overlies Channel belt 3 and contains features that indicate deposition by a low-gradient sinuous fluvial system. The high-velocity channel (Channel belt 3) is overlain by Story 1 of the sinuous channel (Channel belt 4). Story 1 is characterized by coarse to medium-grained, inclined sandstone beds with a channel depth range of 4.1–6.8 m (**Table 2**). The occurrence of convoluted beds and root structures in story 1 of the sinuous fluvial system may have been due to changes in water levels of the fluvial system brought about by seasonality in rainfall due to increasing aridity (Cecil et al., 2003; Cecil and Dulong, 2003; Fielding et al., 2009; Allen et al., 2014; Plink-Björklund, 2015). Overlying Story 1 is the LKC bed, which is overlain by Stories 2 and 3. Stories 2 and 3 are characterized by coarse to medium-grained, inclined sandstone beds with a channel depth range of 3.6–9.6 m (**Table 2**). The paleoslope (0.0001–0.0004) and paleodischarge (408–27,648 m3 /s) values of Channel belt 4 indicate an order of magnitude increase in the maximum paleodischarge rate despite the low channel gradient.

The effect of paleochannel geometry, paleoslope, and paleohydrology on the fluvial channel architecture and depositional style of the MPAF varies. The geometry of the MPAF channel varies in the study area. The channel width and depth, which has been used to compare fluvial channels of arid to humid differing climatic regimes (Fielding et al., 2009; Allen et al., 2014), was used to compare the MPAF channels. The depth ranges for the MPAF channels showed variability among MPAF channel belts: 6.1–13.9 m for Channel belt 1, 3.4–5.6 m for Channel belt 2, 2.8–4.7 m for Channel belt 3, and 3.6–9.6 m for Channel belt 4. The independently estimated paleoslope and paleohydrology, which were compared among the MPAF channel belts (**Table 2** and **Figure 10**), showed that the variation in channel depth of all the channel belts at the study location had developed on a low slope and with variable paleodischarge (**Table 2** and **Figure 10**). The low paleoslope observed in all MPAF channel belt also indicates that fast-flowing events that resulted in the formation of upper plane beds and low amplitude bedforms may have been due to flooding events caused by fluctuation in precipitation intensity

(Kosanke and Cecil, 1996; Cecil et al., 2003; Cecil and Dulong, 2003; Miall, 2014).

# Possible Controls on Fluvial Channel Geometry and Paleohydrology

Eustatic rise and fall of sea level may have controlled paleoslope by changing fluvial base level (Blum and Törnqvist, 2000) across the basin in which a reduction in fluvial gradient due to sealevel rise should result in a reduction of channel flow velocity, while an increase in fluvial gradient should increase channel flow velocity. Glacio-eustatic models suggest fluvial channel sandstones are mainly deposited during glaciation when baselevel is low in the basin and floodplain mudrocks are deposited during interglacial periods, when the base-level is high (Falcon-Lang, 2004; Greb et al., 2008; Haq and Schutter, 2008; Falcon-Lang and Dimichele, 2010). The effects of base-level rise and fall have a direct influence on the water levels and hence the flow depth of the fluvial channel. Sediment accommodation in a fluvial channel is limited by the water level and this is reflected in the thickness of preserved channel sand bodies of the MPAF (Shanley and McCabe, 1994; Currie, 1997; Blum and Törnqvist, 2000; Bhattacharya et al., 2016). The thicker channel depth observed in CB 1 and CB 4 may be due to eustatic base-level rise, while the lower thickness values of CB 2 and CB 4 may be due to eustatic base-level fall. However, the occurrence of water escape structures and rooting features in CB 4 sandstones suggests a variable water level in the fluvial system, which is not consistent with the eustatic base-level fluctuations. This suggests that glacioeustasy played an important role in long-term accommodation succession of the fluvial system, but other factors overprinted the glacio-eustatic control.

Tectonic controls on accommodation and physiography of the Alleghenian foreland may have affected the evolution of MPAF channel belts. Tectonic subsidence and uplift may lead to the increase or decrease of slope and hence fluvial gradient. Relaxation and uplift of a subsided Alleghenian foreland basin would have resulted in an increase of slope (Holbrook and Schumm, 1999; Holbrook et al., 2006), which we interpret to have caused an increase of the fluvial gradient and change in fluvial styles reflected in the CBs 1–3 (**Table 2**). The effect of tectonic subsidence and uplift on the base level is similar to the eustatic effect on base level in a sedimentary basin. However, tectonic processes have third-order cycles (i.e., >1 my) which do not fit the higher frequency, fourth-order cycles (i.e., 0.1–1 my) of the MPAF channel belt (see MPAF age estimates, **Figure 2**). This implies that tectonic influence on MPAF geometry and paleohydrology has been masked by other controls. Additionally, tectonic uplifts which resulted in the formation of the Pangean Mountains may have resulted in the formation of a rainshadow zone, which led to extended periods without precipitation (Greb et al., 2008). However, the duration of tectonic processes does not match the higher frequency processes of the MPAF (DiMichele et al., 2010; Gibling et al., 2014).