Controlling factors on source rock development: implications from 3D stratigraphic modeling of Triassic deposits in the Western Canada Sedimentary Basin Open Access

In the past decade, with the development of shale plays, there was a renewed interest in describing and characterizing in greater details fine-grained and organic-rich deposits (Bruns et al., 2015; Chalmers and Bustin, 2012; Jarvie 2012a, b; Løseth et al., 2011; Passey et al., 2010). Recent studies show that the location of organic accumulations results from the interplay of three controlling factors: primary productivity, preservation and dilution (Bohacs et al., 2005; Slatt and Rodriguez, 2012) and that these factors are strongly impacted by the regional geodynamic settings (Trabucho-Alexandre et al., 2012). Understanding the complex and dynamic interactions between these controlling factors requires a full 3D quantification of the primary productivity, the dilution and the preservation of organic matter at basin scale. Stratigraphic forward modeling provides a powerful tool to achieve this goal.

The aim of the present work is to use the recent development in stratigraphic forward modeling of organic-rich rocks (DionisosFlow and DORS) (Chauveau et al., 2013; Granjeon, 2014) to better understand the impact on organic accumulation of the spatial and temporal evolution of the primary productivity, sedimentation rates and oxygen levels at basin-scale in the Montney and Doig Formations.

Our work focuses on the Lower and Middle Triassic Montney and Doig Formations from the Western Canada Sedimentary Basin. Recent studies (Crombez et al., 2016a) proposed a 3D stratigraphic architecture and conceptual models for the formation of organic rich heterogeneities in the Montney and Doig Formations. Based on these previous works, the present study aims at testing different scenarios of primary productivity, basin restrictions and organic matter transportation through stratigraphic modeling, in order to predict the location of organic-rich heterogeneities and discuss allocyclic controls on primary productivity, dilution and preservation.

The Montney and Doig Formations are present in the foreland basin of the Canadian Cordillera in the subsurface of Alberta and British Columbia (Fig. 1a). These Formations were deposited during the Early and Middle Triassic (Golding et al., 2015b; Orchard and Zonneveld, 2009) on the western margin of Pangaea and recent studies suggest a fore-arc basin geodynamic setting (Beranek and Mortensen, 2011; Rohais et al., 2016). At the time of the deposition, the basin was located in latitudes inferior to 30°N (Fig. 1b). These low latitudes imply a dry and arid climate over the basin during the deposition of the studied interval (Davies et al., 1997). The Lower and Middle Triassic strata mainly consist of siltstone and fine-grained sandstone deposited in wave dominated environments (Armitage, 1962; Davies et al., 1997; Zonneveld et al., 2010). Detailed analysis of the sedimentology highlighted an interval of turbiditic deposits (Moslow and Davies, 1997) and an interval of phosphatic shale deposits (Golding et al., 2014). Recent work showed that the studied interval was deposited in a marine basin with a maximum estimated bathymetry of 300 m (Crombez, 2016; Crombez et al., 2016a). Study on the organic content shows that these formations can be considered as a type II/III marine source rock with actual hydrogen index up to 600 mgHC/gTOC and TOC up to 12 wt% (Riediger, 1997; Riediger et al., 1990). In the Montney and Doig Formations, two intervals present significant organic matter accumulation, the upper Montney (initial TOC of non-thermicaly mature rocks (TOC_ini) up to 8 wt% and average TOC_ini of 2.2 wt%) and the basal Doig phosphates zone (TOC_ini up to 14 wt% and average TOC_ini of 1.1 wt%) (Crombez et al., 2016b).

Recent studies on the stratigraphic architecture of these formations (Crombez, 2016; Crombez et al., 2016a) show that the Montney and Doig/Halfway Formations respectively represent the expression of two second order stratigraphic cycles (sequence A and B). Moreover, the Montney Formation (sequence A) is subdivided in three third order sequences (sequence 1, 2 and 3) that respectively represent the second order TST (transgressive systems tract), HST (highstand systems tract) and FSST (falling stage systems tract); whereas the Doig and Halfway Formations (sequence B) were only subdivided into 2 system tracts: the transgressive system tract 4 (TST4) and the highstand system tract 4 (HST4, Fig. 2) that are interpreted as second order TST and the HST-FSST. The sequence boundary between the two second order cycles is marked by a major erosion that is the result of a second order lowstand sea-level (Fig. 2). In the Montney Formation, preserved strata only present sediment inputs from the East and the South-East whereas in the Doig and Halfway Formations sediment inputs are also coming from the West (Fig. 2).

Nutrient supply in this interval is interpreted to be linked to both continental (sequence A) and marine sources (sequence B). Connection to upwelling cells, and marine nutrient supply in sequence B are interpreted by Crombez et al. (2016b) as resulting from the regional geodynamic evolution. Over this period, anoxic layer development is interpreted to be linked to basin restriction (sequence A3) or to high primary productivities (TST4).

In the Lower and Middle Triassic strata, two significant organic accumulations are present respectively in the sequence 3 and in the TST4. A recent study (Crombez et al., 2016b) proposed three conceptual models (Fig. 2) of basin types to explain the organic matter accumulations in the Montney and Doig Formations:

• type 1: an oxic basin (sequences 1, 2 and HST-FSST 4), with no restriction and a moderate primary productivity that lead to limited accumulation of organic matter (TOC_ini < 3 wt%);

• type 2: a restricted basin (sequence 3), with high basin restriction and moderate productivity that lead to moderate concentration of organic matter in the central part of the basin (1 wt% < TOC_ini < 5 wt%);

• type 3: a high primary productivity basin (TST4), with no restriction and an increased primary productivity that lead to accumulation of very rich interval on the margins of the basin (TOC_ini up to 14 wt%).

The sequence 3 and the TST4 contain organic-rich deposits and they represent a perfect case study to investigate the regional controls on primary productivity, dilution and anoxia.

The isopach, bathymetric and paleogeographic maps used as input data for the model were derived from Crombez (2016). The depositional architecture developed in this study resulted from the correlation of over 400 wells, using the sequence stratigraphic method and concepts defined in Hunt and Tucker (1992) and Catunaenu et al. (2009). It provided maps for the area where the Lower and Middle Triassic strata are still present (Fig. 1). A reconstruction of the missing parts of these formations was carried out, based on our understanding of the regional stratigraphic architecture and of the paleogeographic and geodynamic (Beranek and Mortensen, 2011; Crombez, 2016; Rohais et al., 2016) evolution of the basin.

The aim of this study is (1) to simulate the stratigraphic architecture of the Montney and Doig Formation and quantify the controlling factors (accommodation, sediment supply and transport parameters), (2) to test the effects of different primary productivity values and basin restriction conditions on organic matter distribution and (3) to test different scenarios of the organic matter transport.

DionisosFlow is a process-based multi-lithology diffusion driven forward modeling tool developed by IFPEN (Granjeon, 2014). Based on a diffusion algorithm that simulates the transport of different class of sediments at basin scale, the model allows for testing the relative impact of numerous parameters on the basin fill architecture (e.g. sediment flux, location of the sources, subsidence history, eustatic variations) (Csato et al., 2014; Gvirtzman et al., 2014; Hawie et al., 2015; Leroux et al., 2015). Diffusion of the sediments in DionisosFlow mainly depends on two factors: (1) the gravity: the conversion of potential energy to kinetic energy, (2) the water: the transport capacity of the water, that are calibrated using the diffusion coefficients $Kigravi$⁠, $Kiwater$ respectively linked to the slow creeping transport and water driven displacement. As this diffusion only account for transportation induced by gravity or water fluxes, it can be completed by a transport equation that incorporates wave induced diffusion (associated with $Kiwave$). It is important to notice that theses equations are only relevant over long time period and at basin scale, DionisosFlow cannot simulate single events (e.g. floods, crevasse splay, storm), sedimentary architectural elements (e.g. levees, dunes…) or hydrodynamic processes.

Recently, a new model, DORS, has been plugged into DionisosFlow. The DORS model aim at simulating the production, exportation, transportation and burying of organic matter (Chauveau et al., 2013, 2016). As a matter of simplification, the numerical modeling does not aim at representing the complex nutrient transportation and delivery to the basin, but only the primary productivity variations resulting from these processes. In modern oceans, at the scale of the present study chlorophyll concentration (Sigman and Hain, 2012) tends to show that regardless of the nutrient sources (continental runoff or upwelling cells), the primary productivity is correlated to the distance to shore. In DORS, the primary productivity is therefore defined as a function of the distance to the shore, decreasing toward the basin center. Once produced, the organic matter is degraded while settling through the water column, following the law proposed by Martin et al. (1987). Then the organic matter is transported along the sea floor as silt particles. Lastly, the organic matter is buried and degraded using burial efficiency laws (Burdige, 2007). Here the preservation of the OM is a function of the sea-bottom oxygen concentration, and sedimentation rates. This equation globally represents the different processes (chemical and biological) that occurs on the sea bed (Chauveau et al., 2013, 2016).

The simulated area is 1000 km by 500 km with a 10 × 10 km² grid cells. It is centered on the present-day Montney and Doig subcrop area (Fig. 1). The simulation starts at the base of the Lower Triassic (252.2 Ma) and ends at the top the Middle Triassic [237 Ma, in the present work chronostratigraphic boundaries are taken from Cohen et al. (2013)]. It covers a time interval of 15.2 My with 0.2 My time step (76 time steps).

The model takes into account five sediment classes: sand, silt, shale, bioclastic carbonate and organic matter. In order to constrain the simulation, six paleo water depth/topography maps and five thickness maps were used. The first five paleo water depth maps and four thickness maps were used to define subsidence. The last elevation and thickness maps were used to compute the compaction at present day using compaction laws as a function of burial depth for each sediment class.

The eustatic variations used in the model is presented on the Figure 3. In order to fit the stratigraphic evolution, the eustatic curve from Hardenbol et al. (1998) was modified (Crombez, 2016, long term and medium term on Fig. 3) in order to simulate the regional subsidence/uplift due to the tectonic activity associated with the geodynamic evolution of the basin.

Lastly, the location and rates of the sediment inputs were defined based on paleogeographic reconstructions from Crombez (2016) and Rohais et al. (2016). During the deposition of Montney Formation (sequence A), the sediments mainly come from the eastern, southern and northern part of the model, whereas during the deposition of the Doig-Halfway Formations (sequence B), the sediments mainly come from the East, South and West (Fig. 4). The diffusion coefficients (K) used in the present study are reported in Table 1. These coefficients are calibrated in order to reproduce realistic depositional slopes and sediment partitioning at basin scale (Alzaga-Ruiz et al., 2009; Burgess et al., 2006; Granjeon, 2014).

Three scenarios corresponding to the three aforementioned basin types were developed to investigate the impact of diffusion coefficient variation on organic matter distribution (Tab. 2). In these scenarios, the diffusion coefficient of the organic matter (K_OM) was successively set equal to diffusion coefficients of sand (K_sand), silt (K_silt) and shale (K_shale). This screening test aimed at selecting the optimal diffusion coefficient to be used in the rest of the study.

In a second step, scenarios were designed to test the impact of primary productivity and basin restriction on organic matter accumulation (Tab. 3). In these scenarios, the primary productivity is a function of the distance to the shore: end-values are defined at the shoreline and at 150 km from the shore. The primary productivity is interpolated in between these two extreme values and is not affected by the bathymetry or the emerged area topography. This gradient aims at representing the trend of the primary productivity along continental margin. Primary productivity values defined in the models range from 225–112.5 gC/m²/a to 112.5–50 gC/m²/a. In DORS, the water stratification due to the basin restriction is simulated through a basin restriction coefficients (R_water). This coefficient affects in association with the bathymetry the oxygen renewal in marine environments. It ranges from R_water = 1 (open marine) to R_water = 0.05 (highly restricted) (Chauveau et al., 2013, 2016). The objective of this work is to estimate and discuss the controls on primary productivity and basin restriction conditions required to reproduce the organic accumulation and oxygen conditions described in the three basin types.

The best-fit scenario is presented on Figure 5. On this figure, the color scale represents the variation of a synthetic gamma ray, which is computed based on the relative sediment proportion GR = 4 × %shale + 2 × %silt + 0.1 × %sand). As this synthetic gamma ray is based on sediment proportion, it will not show high gamma ray values that are induced by authigenic uranium enrichment. It will therefore not represent the high gamma ray values linked to high uranium concentration in sequence 3 and TST4 (Crombez et al., 2016b).

The model of Figure 5 reproduces the trends of the stratigraphic architecture of the Montney and Doig Formations highlighted in previous studies: (1) the maximum backstep of the shoreline during the early sequence 2; (2) a long regression from the highstand 2 up to the end of sequence 3; (3) the unidirectional filling of the sedimentary basin during sequence 1, 2 and 3 and the bi-directional filling during sequence 4. Moreover, the Figure 5 highlights typical sedimentary geometries present in the Montney and Doig Formations: (1) prograding clinoforms; (2) sandy turbiditic interval (Davies et al., 1997; Moslow and Davies, 1997).

Errors between the best-fit scenarios and the case study are presented on Figures 6 and 7. The Figure 6 shows a relative error that is generally inferior the 30% in the Montney and inferior to 50% in the Doig and Halfway. The increased error in this last interval is linked to the non-calibration of the Halfway thickness that is not considered as a potential organic rich rock. However, thickness in the Doig Formation was better calibrated and therefore resulted in a lower relative error (lower part of sequence B, Fig. 7). In addition to the maps, Figure 7 illustrates the comparison between real data and pseudo-wells extracted from the simulation at specific well locations. On the wells, the depths of the stratigraphic surfaces are from Crombez (2016), on the synthetic wells the surface are placed over the time step that corresponds to the age of the sequence boundary.

Sedimentation rates have a direct impact on the preservation and the dilution of the organic content (Tyson, 2001), it is therefore important that for each sequence, at basin scale, thickness errors are inferior to 25 m. On this figure, it is apparent that some of the sedimentological details present on the well logs are not represented in the model. Despite these differences, the best fit model successfully reproduces overall stratigraphic trends and sediment repartition; with reasonable thickness errors. It can therefore be used to further investigate the controls on organic matter accumulations.

The results of the experiments summarized in Table 2 are presented on Figure 8. This figure shows similar impact of the diffusion coefficient changes on the three basin types: the organic matter tends to concentrate in the center of the basin. On Figure 8, the concentration of the organic matter in the central part of the basin is more pronounced in the type 3 configuration and weaker in type 1 configuration. The Figure 9 presents a dip-oriented profile of organic matter concentration (location on Fig. 8) together with associated oxygenation and bathymetry for each configuration.

The scenarios with high diffusion show a relative depletion (compared to low diffusion scenarios) of organic matter in areas located at less than 140 km of the shoreline and a concentration of the organic deposits beyond 140 km. The amplitude of organic matter variations is higher in type 3 (TOC_ini up to 4.5 wt%) than in type 2 (TOC_ini up to 2.5 wt%) and type 1 (TOC_ini up to 0.5 wt%) basins. However, in the deepest part of the basin, the relative increase of organic matter concentration is close to 70% in both type 2 and type 3 basins, which shows that the relative impact of an increase of organic matter diffusion is the same regardless of the primary productivity. Here, the increase of the transport coefficient induced a displacement of the maximum accumulation area (Fig. 9): between K_sand and K_shale, the maximum accumulation area, shifted 50 km seaward. Lastly, on the three configurations, the Figure 9 shows no change in the bathymetry profile and only minor change in the oxygen concentration.

Various studies focused on the sedimentary processes that drive the formation of mudstone and others fine-grained sediments (Ghadeer and Macquaker, 2011; Macquaker and Bohacs, 2007; Macquaker et al., 2010; Schieber and Southard, 2009; Schieber et al., 2007). All these studies highlight that mudstone and others fine-grained sedimentary rocks get mainly deposited in environments that still present unidirectional or oscillatory currents. In marine source rocks, the organic fraction that reaches the water sediment interface is mainly composed of faecal pellets, flocculates or other aggregates (Ghadeer and Macquaker, 2011) called the “marine snow” (Alldredge and Silver, 1988). Therefore, the organic fraction presents on the sea bed is composed of various particles with various sizes and densities (Alldredge and Silver, 1988; Wakeham and Lee, 1993). Recent study (Schieber et al., 2007) suggested that the organic sediments may have similar transport characteristics than fine-grained sandstone, but it does not take into account the wide variety of particles that compose the sedimentary organic matter.

In DORS, the organic matter is considered as a sediment class with its own transport coefficient. In this study, we tested various diffusivity coefficients in order to properly reproduce the organic matter accumulations (from K_sand to K_shale, Figs. 8 and 9). In the center of the basin, the organic matter concentration was approximately 70% higher when using K_sand than when using K_shale, regardless of the primary productivity or oxygenation of the water (Fig. 9). Although DORS does not reproduce explicitly the flocculation or aggregation processes, we found that using a transport coefficient of silt to fine-grained sand is a good approximation for modeling basin-scale distribution of organic matter.

The values of the sediment and water fluxes used in the best-fit model are presented in the Figure 9A. Average sediment inputs are higher in the Montney Formation (sequence A: 28.10³ km³/Ma) than in the Doig-Halfway Formations (sequence B: 5.10³ km³/Ma). Within the Montney Formation, the total sediment input increases from sequence 1 up to sequence 3 (from 16.10³ km³/Ma to 43.10³ km³/Ma) and the limit between the Montney and the Doig Formation is marked by a sharp drop of total sediment inputs. In sequence A, the relative proportion of sediment delivered from the three sources (N, E, S) is not constant (Fig. 9A). During sequence 1, the maximum sediment input is coming from the North, whereas in sequence 2 the maximum sediment input is coming from the East. In the Montney Formation, the sediments coming from the South only contribute 15 to 20% of the total sediment supply, whereas the source from North or from East provides up to 60% of the sediments inputs. During the deposition of the Doig Formation, sources from West represent more than 40% of the sediment input whereas the sources from East and South respectively carry 35% and 25% of the sediments. Similar trends are also visible in water discharge (Fig. 9B), the Montney Formation (sequence A) shows water discharge from 1.10³ m³/s to 1.10⁵ m³/s. The Montney/Doig boundary presents a major turnover with a sharp drop in the total water discharge (6.10² m³/s).

In order to investigate the impact of the primary productivity and basin restriction on anoxia and organic matter distribution, nine scenarios were tested on the best fit model. Bottom water oxygen concentration and organic matter maps were extracted from the model at intervals representative of type 1, type 2 and type 3 basins. The bathymetries corresponding to those intervals are presented on Figure 11. It shows that the bathymetry in type 1 and 2 are deeper than in type 3. It also shows a major change in the physiography of the basin: type 1 and 2 intervals only present sediment supply from the East whereas the interval of type 3 presents sediment supply from the East and the West. The Figures 12–14 present the benthic oxygen and organic matter concentration map resulting from the sensibility analysis.

On Figure 12, it is apparent that all scenarios with high restriction and the model with medium restriction and medium to high primary productivity present anoxic conditions in the central part of the basin. On Figure 12, only LP-0.15 (low organic productivity – basin restriction coefficient = 0.15 ∼ moderately restricted) and the three low restriction scenarios present oxic conditions. Here, only scenarios with high or medium primary productivity, associated with medium to high restriction present organic accumulation higher than 1.5 wt% (with a maximum of 4 wt% in HP-0.05 scenario). The Figure 12 shows that the model with MP-0.15 and HP-1 present the same distribution of organic matter than the simulation with MP-1 and LP-0.05. Among the nine scenarios tested, only the three models with low restriction resulted in a fully oxic environment with low organic accumulations (Fig. 12).

The Figure 13 shows similar results than Figure 12 due to the similar shape of the basin (see Fig. 10): the three models with high restriction present anoxia as well as the scenarios with medium restriction associated with medium to high primary productivity. Oxic conditions are prevalent in models with low restriction and in the model with moderate restriction and low productivity (LP-0.15 scenario). The organic content of the nine scenarios showed in Figure 13 are significantly lower than the one displayed in Figure 12 because of the higher sedimentation rates and the resulting dilution of organic matter (Fig. 9). On Figure 13 only two scenarios with high restriction and high to moderate productivity (HP-0.05 and MP-0.05) result in anoxic deposits and TOC_ini up to 1 wt% in scenario HP-0.05.

On Figure 14 only two scenarios show evidences of anoxia (HP and MP-0.05) and two others show local anoxia and widespread dysoxia (HP-015 and LP-0.05). In type 3 basin (Fig. 14), the range of organic richness is larger than in type 1 and 2 basins (Figs. 12 and 13). Here, the maximum TOC_ini reaches 10 wt% in the HP-0.05 scenario. Four scenarios (MP-0.05, HP-0.15, MP-0.15 and HP-1) result in significant organic deposits (> 2 wt% up to 10 wt%) and three only present limited organic accumulations (MP and LP-1 and LP-0.15). In all the scenarios, organic rich areas are located along the Western margin of the basin. Among the 9 scenarios only two present the organic distribution and oxygen condition consistent with our definition of a type 3 basin: HP-0.05 and MP-0.05.

Dilution of organic matter by non-organic particles is considered to be one of the main controlling factors of source rock accumulation (Bohacs et al., 2005; Katz, 2005; Tyson, 2001). In the present study, the use of DionisosFlow provides a way to investigate the sediment supplies variations through time at basin scale (Hawie et al., 2015; Leroux et al., 2015). The relations between water discharge, sedimentation load, drainage area and topography of modern rivers has been extensively studied (Dai and Trenberth, 2002; Milliman and Meade, 1983; Milliman and Syvitski, 1992; Syvitski et al., 2003; Syvitski et al., 2005).

The Figure 15 illustrates the relation between the drainage area and the water discharge from more than 160 rivers (Dai and Trenberth, 2002). On this figure is plotted the average water discharge of the sediment sources used in the best-fit model. As no major climatic changes from dry to humid is known to occur during the Lower and Middle Triassic on Western Canada (Davies, 1997a; Davies, 1997b; Golonka et al., 1994; Hallam, 1985; Sellwood and Valdes, 2006), the variations of water discharge are most likely linked to change of the drainage area. During sequence A, the increase of water discharge and sediment supply from the North, East and South (Fig. 10) suggests an expansion of the drainage areas from the continent (Fig. 15) (Golding et al., 2015a). At the transition between sequence A and B, a sharp drop of water discharge and sediment supply suggest a significant reduction of the continental drainage area to the east (Milliman and Syvitski, 1992). However, the drop of water discharge (divided by 370) is of larger magnitude than the drop of sediment supply (divided by 15), resulting in a higher sediment concentration in rivers during sequence B. This may be explained by a topographic rise of this area, resulting in a steeper topography and higher erosion rates of the drainage area during sequence B. Concomitantly, the exhumation of sediment sources from the proto Canadian Cordillera to the west, further indicates significant paleogeographic and geodynamic transformations at that time.

Edwards et al. (1994) show the occurrence of Triassic strata in two basins separated by the Sweetgrass arch (Kent and Christopher, 1994): (1) the Alberta basin and (2) the Williston basin. In the Alberta basin, recent publications suggest that the Triassic deposition took place in a fore-arc basin (Beranek and Mortensen, 2011; Rohais et al., 2016). The DionisosFlow model presented in this study suggests that during the sequence A (Lower Triassic) a significant part of the sediment supply derived from the continent, whereas during the sequence B (Middle Triassic) sediments came from the proto Canadian Cordillera and from the Canadian continent, but with a narrower and steeper drainage area. In the Williston basin, detailed stratigraphy (Cumming, 1956) shows that the continental Lower Watrous Formation was deposited during the Middle and Upper Triassic and rest upon a major unconformity (Carlson, 1968; Fowler and Nisbet, 1985; Kent, 1994). Here, the deposition of continental strata on top of an unconformity highlights a major change of the regional subsidence rates.

In our case study, the geodynamic evolution seems to have a significant impact on the sources and sinks of Lower and Middle Triassic strata. During sequence A, the Triassic strata are sourced by a drainage area that may extend eastward the Canadian Shield (Carlson, 1968; Kent, 1994). During sequence B due to the evolution of the proto Canadian Cordillera (Rohais et al., 2016), the basin is most-likely sourced by the cordillera and to a lesser extent by the continent. Here, the reduction of the drainage is interpreted to be linked to the rise of a structural high that will split a widespread drainage area into two catchment areas: one flowing to the Alberta basin and a second flowing to the Williston basin. The rise of this topographic high will increase the erosion and induce a drop in the sediment supply and water discharge due to the reduced drainage area. In the same time, the geodynamic evolution of the proto Canadian Cordillera affected the subsidence rates in Western Canada by: (1) reducing the subsidence rate in the Alberta Basin, (2) provoking the rise of a relative structural high between the Alberta basin and the Williston basin, that reduced the sediment and water delivery to the Alberta basin and (3) increasing the subsidence in the Williston basin, allowing the deposition of continental red beds.

Huc et al. (2005) highlighted the first order control of the geodynamic settings on source rock accumulation because of their link to volcanic activities increases that affect primary productivity. The present study shows another important control of the geodynamic and physiographic evolution on the sedimentary delivery in the basin at the scale of second order stratigraphic cycles:

• in sequence A, the extensive drainage area resulted in high sedimentation rates (> 150 m/Ma) that caused a high dilution of the organic matter (Fig. 13). In anoxic environment, this dilution results in a decrease of the organic concentration whereas in oxic environment it significantly increases the preservation (Tyson, 2001);

• in sequence B, the reduced drainage area resulted in low sedimentation rates (< 50 m/Ma). In this scenario, anoxia is required to generate significant organic matter accumulation (Fig. 14).

In sedimentary basins, anoxia can be due to high productivity (e.g. the Peruvian margin or the Namibian margin), to a restriction of the basin (e.g. the Black Sea) or to a combination of both (Arthur and Sageman, 1994; Brumsack, 2006; Sageman and Lyons, 2003). In the Montney and Doig Formations, three types of basin are present (Crombez et al., 2016b):

• an oxic basin;

• a restricted basin;

• a high productivity basin.

In these three types of basin, the organic matter distribution is not the same:

• in oxic basin, little organic matter is preserved (TOC_ini < 3 wt%);

• in restricted basin, moderate accumulation of organic matter is preserved in the central part of the basin (TOC_ini < 5 wt%);

• in high productivity basin, large accumulation of organic matter occurs on the edge of the basin (TOC_ini < 14 wt%).

The sensitivity tests performed with the new DORS model (Figs. 12–14) allows for a quantification of organic productivity in each basin types:

• type 1: the results of the modeling show that even a moderate restriction induces dysoxia (Fig. 12). The results also show that a productivity of 150–75 gC/m²/y is too low to generate 1.5 wt.% of TOC_ini, and that a productivity of at least 225–125 gC/m²/y associated with basin restriction is needed to decrease the oxygen content in the bottom water (Fig. 12). In order to obtain the correct organic content described by previous studies, the primary productivity must be increased to 250–125 gC/m²/y which is a moderate rate for coastal areas (Saliot, 1994) and is not high enough to induce dysoxia;

• type 2: the modeling results show that for primary productivity higher than 150–75 gC/m²/y associated with strong and moderate restriction, the basin becomes anoxic (Fig. 13). As for type 1, the primary productivity in our simulations was underestimated. An increase to 250–125 gC/m²/y is therefore necessary to fit with the organic content in type 2 basin;

• type 3: the results of the stratigraphic modeling (Fig. 14) show that even the highest primary productivity (225–112.5 gC/m²/y) does not induced dysoxia if not associated with strong restriction. Crombez et al. (2016b) suggested the possibility for a restriction in type 3 basin, which appears to be necessary. The results of the modeling show that if a moderate restriction with a production of 225–112.5gC/m²/y generates anoxia locally, it mostly results in dysoxia.

As a perspective, additional scenarios could be tested: in type 1 and 2 basins primary productivity values were underestimated in our simulations and in type 3 basin higher primary production scenarios could to be used (400–200 and 350–175 gC/m²/y). If these new scenarios do not induce anoxia in the basin, then scenarios with strong restriction will be required for type 3 basin. Crombez et al. (2016b) did not discuss a possible increase of the primary productivity between sequence 1, 2 and 3. The present study emphasizes the increase of sediment supply from sequence 1 to sequence 3. If the sediment sources remain the same throughout sequence A, the increase of the sediment delivery will therefore induce an increase of the total volume of nutrient delivered. This increase of nutrient delivery may induce increase of primary productivity from the sequence 1 to sequence 3. Moreover, Figure 14 shows that with a primary productivity of 225–112.5 gC/m²/y, TOC_ini values are up to 10 wt.% which is close to the value published in previous works. Therefore, the increase of the primary productivity to 400–200 gC/m²/y might be too high inducing an over-accumulation of organic matter.

A previous study on the stratigraphic architecture of the Montney and Doig Formations show that this interval is composed of two second order sequences (Crombez, 2016). In these two sequences, organic rich intervals are not located in the same systems tracts: in sequence A organic rich layers are located in the FSST, whereas in the sequence B they are located in the TST.

In the present work, the stratigraphic modeling emphasized a major drop of the fluvial discharge and sediment inputs at the boundary between sequence A and sequence B. This reduction is interpreted to be linked to a reduction of the drainage area caused by the relative rise of an arch on the continent. In the Lower and Middle Triassic strata of the Alberta basin, it is apparent that the regional geodynamic evolution of Western Canada and particularly the evolution of the Proto Canadian Cordillera (Beranek and Mortensen, 2011; Rohais et al., 2016) controls the dilution rates in sequence A and sequence B. Based on our simulations, we can estimate that during sequence A the average sediment fluxes was 23.10³ km³/Ma whereas during sequence B it was only 5.10³ km³/Ma. This important drop of detrital input resulted in drastically different dilutions between the two second order sequences: organic particles in sequence A were 4.5 times more diluted than in sequence B.

Crombez et al. (2016b) concluded that in sequence A primary productivity is controlled by continental nutrients delivery whereas in sequence B, the high productivity of the TST is linked to marine nutrients brought through upwelling. In this study the primary productivity profile are assumed to have similar trends but different values. The present work highlighted an increase of the sediment fluxes from sequence 1 to sequence 3. In sequence A, as the primary productivity is controlled by continental nutrient inputs (Crombez et al., 2016b) the sediment fluxes increase probably induces a small primary productivity increase. In this interval, the potential moderate increase of the primary productivity is therefore connected to the regional geodynamical evolution of the area that controlled the size of the drainage area. In sequence B, the increase of nutrient inputs is interpreted to be related to upwelling activities (Crombez et al., 2016b). The sequence boundary between sequence A and sequence B is interpreted to be associated with a major rearrangement of the basin paleogeography (Rohais et al., 2016). Consequently, the high primary productivity in TST4 is linked to the regional geodynamic evolution that connected the basin to upwelling cells. In sequence B, the fact that high productivity only occurs in the TST suggests that during the geodynamically induced transgression, adequate water exchange made it possible to deliver the deep marine nutrient in the basin. Wignall and Newton (2001) proposed a model of “transgressive nearshore black shale” where high primary productivity is controlled by transgressive trends. The controls on the organic accumulation in the TST4 of our model are probably more complex as it encompasses longer time scale and thicker deposits.

In the present case study we propose two different models for generating anoxia. In the FSST of sequence A, the anoxia is interpreted to be triggered by a major basin restriction (Crombez et al., 2016b) whereas this study shows that the anoxia in the TST of sequence B is linked to both primary productivity and a moderate basin restriction. In both models, the development of anoxia is associated to the regional geodynamic evolution of western Canada:

• the major regression at the end of sequence A is induced by a regional tectonic uplift;

• the physiography of the basin during sequence A and sequence B is linked to the geodynamical contex;

• the primary productivity that generate anoxia in the TST of sequence B is linked the geodynamic evolution between sequence A and B.

This interpretation contrasts with the model of Algeo et al. (2008) that shows a major control of the climate and water delivery on anoxia, but is more similar to the model proposed by Tyson et al. (1979) and Wignall (1991) where anoxia is controlled by the basin physiography.

Based on previous works on the Montney and Doig Formations and on a 3D stratigraphic forward modeling this study shows that:

• a major reduction of the continental sources drainage area occurs between sequence A and B. It induces a major drop in the sediment inputs during the deposition of the sequence B. This shift is interpreted to result from a major geodynamic change that occurs in western Canada at the Lower/Middle Triassic boundary (Rohais et al., 2016). During sequence A, the increasing sediment supply results in higher nutrient delivery and organic matter dilution;

• an increase of the basin restriction is needed in type 3 basin (with high primary productivity) to trigger anoxia. Here both restrictions are induced by the basin physiography that is linked to the geodynamic evolution of the area;

• the estimated productivity profile is 250–125 gC/m²/y in type 1 (oxic basin) and 2 (restricted basin) and 400–200 gC/m²/y in the type 3 basin (Doig phosphates zone). These values are respectively characteristic of passive margin and upwelling areas. Here additional scenarios are required to test the accuracy of these values;

• in this case study, the distribution of the organic rich layer is controlled by the regional geodynamical evolution: it controls the sediments fluxes variations, the basin restriction and the nutrient delivery to the basin. In the present study, we show that geodynamic have an impact on marine (upwelling) and continental (runoff) nutrient delivery. These conclusions are concordant with Trabucho-Alexandre et al. (2012) but drastically differ from the transgressive model of Slatt and Rodriguez (2012). In the present study, as no climatic changes occur during the deposition of the Lower and Middle Triassic strata, geodynamic framework is the first control on source rock development, the transgressive and regressive trends only play a second role.

Acknowledgments go to IFP Energies nouvelles and the ISTeP for funding this project and to all those who have contributed to it: IFP Technologies (Canada) Inc; and Alberta Energy regulator and the Geological Survey of Canada, BC Oil and Gas Commission for the access to core data and samples. Authors also would like to thank the anonymous reviewer for his helpful comments.