Sedapp: a non-linear diffusion-based forward stratigraphic model for shallow marine environments

The formation of stratigraphy in shallow marine environments has long been an important topic 15 within the geologic community. Although many advances have been made in the field of forward stratigraphic modelling (FSM), there are still some shortcomings to the existing models. In this work, the authors present our recent development and application of Sedapp: a new non-linear open-source R code for FSM. This code uses an integrated depth-distance related function as the expression of the transport coefficient to underpin the FSM with more along-shore details. In addition to conventional 20 parameters, a negative-feedback sediment supply rate and a differentiated deposition-erosion ratio were also introduced. All parameters were implemented in a non-linear manner. Sedapp is a 3D (2DH) tool while also capable of 2D (1DH) scenarios. Two simplified case studies were conducted. The results show that Sedapp can not only assist in geologic interpretation, but is also an efficient tool for internal architecture predictions. 25


Introduction
Shallow marine areas are among the most active environments for sedimentation, where sea level, tectonism, climate all influence the interactions between land and sea. Stratigraphic formation in this 30 environment has long been an important topic within the geoscience community, which has directly resulted in the emergence of sequence stratigraphy (Haq et al., 1987;Li et al., 2015;Catuneanu, 2019).
Traditional qualitative methods on this subject have made great advances in the past half century, but it is difficult to test the validity and internal consistency of a new concept and are less likely to raise counter-intuitive ideas which may actually be true (Burgess, 2012;Burgess and Prince, 2015).
FSM deals mainly with long-term geomorphologic/stratigraphic dynamics (Paola, 2000). It is slightly different from sediment fluid-flow models, which deal more with the fluid-flow dynamics by solving modified Navier-Stokes equations within a full study domain (generally shallow water equations, e.g. HydroSedFoam of Zhu et al., 2019or Delft3D in Ramos et al., 2019. FSM models can 45 be classified into two types, i.e. rule-based models (based on geometric or fuzzy logic) and equation-based models (Paola, 2000;Syvitski and Hutton, 2001;Burgess et al., 2012;Sacchi et al., 2015). The first type easily captures essential features and is less time-intensive, while it does a relatively poor job of demonstrating predictability and revealing the physical processes (Strobel et al., 1989;Kendall et al., 1991;Burgess, 2012). The latter type is also known as deductive models, which 50 are process-based and solve governing equations (Kaufman et al., 1991;Rivenaes, 1992;Granjeon and Joseph, 1999;Griffiiths et al., 2001;Hutton and Syvitski, 2008;Li et al., 2020). For these long-term processes, sediment flux is usually assumed to be proportional to the topographic gradient. Thus, a diffusion equation like Eq. (1) is generally used to derive the governing equations in FSM models (Salles et al., 2018;Ding et al., 2019). 55 https://doi.org/10.5194/gmd-2020-256 Preprint. Discussion started: 17 November 2020 c Author(s) 2020. CC BY 4.0 License.
where h denotes the topography, t denotes the time and Γ denotes the transport coefficient.
Diffusion-based models are good at modeling scaled stratigraphic sequence (relative larger scales, e.g. clinoform formation) processes. Γ in Eq. (1) can be defined using different values for different environments (Zhang et al., 2020). Various Γ types are used based on different needs and environments 60 (Rivenaes, 1997;Zhang et al., 2020). Models with constant Γ values are usually called linear models; otherwise, they are known as non-linear models.
Although the sediment diffusion assumption is considered a practical representation of long-term slope processes, it is still too simplistic when the Γ is used as a constant because natural agents such as air and water, phenomena such as mass wasting, and biological agents actually move sediment at rates 65 that are not determined solely by slope (Salles et al., 2018). This severely limits the application scope of linear diffusion-based models. On the contrary, non-linear models are relatively more flexible. Many non-linear models define the transport coefficient using water depth-related functions (e.g., Clarke et al .1983;Kaufman,1991). These water depth models work well in general coastal zones. However, in the shallow marine environments with river injection, the water depth models are usually not applicable.

70
Depositional processes around the river mouth are more active than those at a distance, even when they are at the same water depth. Hence, it is difficult for water depth models to reveal along-shore variability, especially in 3D scenarios (actually 2D-H, with two horizontal dimensions and the elevation H). Because of this reason, many hybrid hydrodynamic-diffusion models were proposed. For example, river plume was introduced to differentiate the suspended sediment fluxes along the coastline (Syvitski 75 and Hutton,2001;Hutton and Syvitski, 2008;Dalman and Weltje., 2012). This treatment realized the modelling of the convex shapes of delta out from the river mouths. However, the computational load was also significantly increased because of the newly added advection-diffusion process. The introduction of many hydrodynamic parameters also increased the difficulty of its usage.
In addition, many existing models are not free or open-source, making it difficult for people to 80 reproduce and improve them. Some codes, although free, can only be used in Linux systems, which makes them inconvenient to use in most PC terminals.
In this paper, we propose a new non-linear model, which is expected to overcome the shortcomings of the existing models. Along with some other features, this model is integrated into a https://doi.org/10.5194/gmd-2020-256 Preprint.
where Fi is the fraction of the ith class of lithology, h is elevation, t is time, ∇ is the nabla operator, Der is a user-defined parameter denoting the ratio of deposition to erosion (it can be a scalar, vector or tensor value depending on its temporal and spatial variability), Γi is the diffusion coefficient for the ith 95 class of lithology, and q is the source term that is a function of coordinates and time (the source term is used only for endogenetic sedimentation, especially carbonates. If endogenetic sedimentation is ignored, the source term can be left out). Among them, h and Fi are the primary unknowns.
Note that Γi cannot be outside the parentheses, because they are not constants but rather functions of spatial coordinates and time. The expression of a general Γ can be expressed as: where α/αwd are preexponential factors (L 2 /T) , η/ηwd are distance indexes (no dimension), β/βwd are spatial scale factors (  L or wd  L ), and ε is an adjustment factor (L 2 /T) reflecting the environment energy. In particular, distance function D=D(x,t) and water depth function Wd(x,t) change with spatial coordinates and time, and they work for the marine portion only.
where x and y are spatial coordinates. This is especially suitable for cases dealing only with two classes 110 of lithology for simplicity, where Γ1 is the transport coefficient for sand and Γ2 is the transport coefficient for mud.
For 2D (1D-H) scenarios, especially along the section line through the river mouth, the distance related term is generally larger than the water depth related term, so the latter term within the max function in Eq. (4) is usually omitted. For convenience in coding, also ignoring the endogenetic 115 sedimentation, Eq. (5), Eq. (6) and Eq. (4) can be simplified into: The joint effect of c and E in Eq. (9) is equivalent to that of β in Eq. (4). The variable c here, with 120 a dimension of L -1 , is mainly used to facilitate the scale of distance and differentiate the transport characteristics of different sediment types (e.g., sand and mud). E is a dimensionless constant that represents hydraulic characteristic energy.

Code Implementation
Sedapp was written in the R language and its solution procedure was based on the finite volume 125 method (FVM), which has the desired property of local mass conservation and has a clear physical meaning (Versteeg and Malalasekera, 2007;Moukalled et al., 2016;Liu P. et al., 2017). The cell-centered variable arrangement method was used to store the unknowns at the grid element centroids. The non-linearity was implemented through stepwise iteration (Fig.1).
The brief work-flow within a single time step is as below: 130 1) Implement user-defined tectonic subsidence and update the topography; 2) Implement user-defined sea level and identify/update the shoreline location; 3) Solve the differential deposition/erosion function; 4) Implement the compaction and isostatic subsidence. https://doi.org/10.5194/gmd-2020-256 Preprint. Discussion started: 17 November 2020 c Author(s) 2020. CC BY 4.0 License.
Step 3) is a major step. According to the hypothesis of diffusion-based FSM models, the change 135 rate (by either deposition or erosion) is proportional to the gradient of the slope (Fernandes et al., 1997;Pelletier, 2013). If we use the diffusion equation/law directly without any differential treatments between deposition and erosion (in other words, Der is held at 1), it will be contrary to the geological knowledge that deposition and erosion processes are two distinct processes with different rates. Hence, the max() function is used as in Eq. (2). Generally, the erosion process occurs at a different rate than 140 deposition (also called erosion constraints, see Galy and France-Lanord,2001), so Der is usually not equal to 1. For example, if we wanted the erosion rate to be only 1/100 of the deposition rate, the Der can be set to 100. In this case, if it is a deposition process (namely the non-erosion case is desired, Der can be set to a very large value.

145
Generally, sediment supply rate cannot be directly defined through boundary condition settings, since the latter can only determine the boundary slope. Therefore, Sedapp uses a negative-feedback strategy to define the sediment supply rate. At each time step, the total amount of deposition within a step is first calculated using the previously defined test  , and then the adjusted mod α is calculated by Eq. (10): where αmod denotes the modified α of this time step; Vexpected denotes the expected sediment increment, namely the sediment supply rate; and Vtest denotes the computed sediment increment with αtest.

155
The nonlinear transport coefficient is a feature of Sedapp. Sedapp's transport coefficient uses a function of both the distance from the estuary and the water depth. This feature makes it easier to simulate fluvial-deltaic processes in 3D scenarios, which can reflect changes along the shore. Even in 2D cases, this feature also has some advantages (see the discussion section for details).
range when the total amount of sediment is fixed. For example, the c of mud is usually set to 50%-85% of sand, thus reflecting the differential deposition of sand and mud. In addition, the environment energy  can also influence the sediment travel distance that a larger  can make the sediment travel further.
As sedimentation progresses, the position of the estuary may change, so the distance from the estuary is updated at each time step to achieve the nonlinearity of  .

Differential and customizable deposition/erosion rate
During the actual deposition process, the properties of the lower strata (such as compaction degree, lithology, and age, etc.), as well as some external environmental factors (such as temperature, humidity and pH value, etc.), will affect the erosion rate. Therefore, the customized treatment of erosion rate is another Sedapp characteristic.

170
In Sedapp, the deposition rate is a parameter that can be specified directly (for the adjustment process see section 2.2). Furthermore, the Der parameter is a user-defined parameter that controls the ratio of deposition rate to erosion rate. When Der is 1, the deposition rate is equal to the denudation rate ( Fig.2a), and when Der value is 10 or 100, denudation is significantly weakened (Fig.2b). Theoretically, if the value of Der is large enough, it is equivalent to completely eliminating the denudation effect. Der 175 values should be customized according to the actual situation.

Customizable compaction
Compaction is an important geological process after sediment deposition, especially when the sediment thickness is very high. In Sedapp, the compaction process can be easily realized by setting the composition of lithology and porosity curves.

180
In this paper, we designed a pyramid-shaped mountain simulation commonly used by other researchers (as shown in Fig.3, see Rivenaes,1992 andYuan et al.,2019 for reference). The Der value was set to 1. The sediment supply ratio of sand and mud was set to 1:1, and the porosity curve was set as shown in Fig.3d. After simulation, the top of the pyramid was denuded and the foot of the pyramid had deposited sediment of a given thickness.

185
To illustrate the effect of compaction, Sedapp introduces a scale factor that can enlarge the longitudinal scale. Fig.3a shows the original compression scale (that is, the scale factor was equal to 1), and the scale factors in Fig.3b and Fig.3c were 100 and 1000, respectively. It can be seen that sediment https://doi.org/10.5194/gmd-2020-256 Preprint. Discussion started: 17 November 2020 c Author(s) 2020. CC BY 4.0 License.
thickness at the foot of the pyramid in Fig.3c was significantly smaller than that in Fig.3a. The factors that caused these differences were not only depth but also the proportion of sand and mudstone and the 190 shape of depth-porosity curves, which can be easily adapted to different scenarios by modifying the lithologic proportion and porosity-depth functions in Sedapp.

Verification of Sedapp
To identify how well the algorithm works within geological context, some simple benchmark simulations are given below.

Typical stacking patterns
Typical stacking patterns including forced regression, normal regression, and transgression can be formed (Fig.4) by fixing sediment supply while controlling the adjusted sea level rise rate.
During the period of sea-level decline, the shoreline moved seaward, and the onlap points also moved seaward and form the offlap and downlap stratigraphic termination structures (Fig.4a). During 200 slow sea-level rise, the shoreline continued to move seaward, but the onlap points started to move landward, forming an onlap termination structure. At the other end, the downlap structure continued to exist. During rapid sea-level rise, the shoreline started to move landward and the onlap points also moved landward. At this time, downlap structure did not exist above the slope break, but may have existed below the slope break.

Typical two-cycle scenario
To demonstrate the complete base level changing process, this paper designed a simulation with two full sinusoidal cycles as shown in Fig.5. In the first cycle, the shoreline dropped and moved seaward. Then it slowly rose and gradually moved landward until it reached the highest point and tended to stabilize. The water depth of deposition in the strata gradually deepened from left to right on 210 the marine side (Fig.5a), and the sandy content reached a maximum around the shoreline (Fig.5b) near the shoreline. In the strata on the land side, the sand content was stratified. The sand content was relatively large during the early transgression and subsequently relatively small. The second cycle was located above the first cycle and continued the same characteristics as the first cycle, but the deposition range was enlarged and the average single layer thickness was thinner. The projection of the simulation results on the x-y plane clearly shows the variation characteristics of the along the shore. When t = 0, the shoreline was a straight line, and the channel was in the middle of the shoreline. As time went on, the river mouth continued to move forward. From 0 to 2 Ma, the channel first swung to the north, then to the south, and shoreline began to bulge slightly 230 towards the sea side. From 2Ma, the channel continued to swing southward, until the time approached 4Ma and the river mouth began to turn north slowly. From 4Ma to 6Ma, the channel continued to swing northward, and the convex part towards the sea side became more and more obvious. From 6Ma to 8Ma, the channel continued the previous trend, while the convex shoreline became asymmetry (an increasing skewness to the north). From 8 Ma to 10 Ma, the principal line of the channel moved 235 southward, and the convex shoreline gradually returned to symmetry (Fig.7).
The simulation results also show some interesting features on longitudinal sections. Two sections (y = 75m and y = 125m) perpendicular to the shoreline direction are selected (see Fig.7f for the position of the section line). The two sections are located on the north and south sides of the main channel. The distance between the channel and the two sections is varying. In the southern profile (y = 240 75m), from 4 Ma to 10 Ma, the isochronous lines of the formation changes from sparse to dense, and then from dense to sparse (i.e., the thickness of a single clinoform changes from thick to thin first and then from thin to thick) (Fig.8). This is completely contrary to that observed in the northern profile (y = 125m). From 4Ma to 10Ma, the isochronous lines first changes from dense to sparse, and then from https://doi.org/10.5194/gmd-2020-256 Preprint. Discussion started: 17 November 2020 c Author(s) 2020. CC BY 4.0 License.
sparse to dense, reflecting that the deposition rate first increases and then decreases (Fig.9).

245
Under the parameters shown in Tab. 1, due to the existence of estuaries, shoreline will bulge towards the sea side. A closer distance to the river mouth could result in a higher sedimentation rate and a greater shoreline advancing speed. From 2 Ma, the convex shape of the shoreline towards the sea side became more and more apparent, similar to the morphology of some real-world Deltas (Fig.10).

2) Model 2 250
This code can be applied not only to marginal marine environments but also to the continental fault basins. Taking the 3 + 4 sand groups of the third member of Shahejie Formation in the Gaobei slope belt of Nanpu Sag in Bohai Bay Basin as an example, we conducted a simplified 2D real case study. The basic geological background is as follows: During the deposition period of this set of strata, the normal fault tectonic movement in the north of the sag was active, which was the main controlling 255 factor leading to the increase of accommodation space. At the same time, the terrigenous clasts came from the north is sufficient, and the basin was in a balanced state (Li et al., 2018). According to the geological background, a simplified reconstruction model (Model 2) was designed, which assumed that the subsidence rate of the boundary fault and sediment supply rate is constant, neglected the effect of isostasy, and considered the effect of sediment compaction.

260
The simulation results are shown in Fig.11. From the perspective of temporal and spatial stratigraphy, the shoreline mainly moved towards the sag center during the early stage, and then moved back to the land side. The deepest water depth occurs in the middle south part at 2 Ma (Fig.11a). This shoreline phenomenon is usually called autoretreat (Muto and Steel, 2002). The sand fraction section shows that the steep slope belt in the north is richer in sand content than the south (Fig.11b). The 265 porosity section shows that the porosity generally decreases from bottom to top. The porosity also varies horizontally, especially when the depth is deeper than 800 m. The porosity in the north is larger than that in the south.
Due to the over-simplified assumptions, the simulation results are not necessarily be consistent with every practical borehole. However, the general trends are revealed through the simulation, which 270 can strengthen or improve our existing understanding and guide us to seize the main direction. Also, the facies simulation results were in good agreement with the Sedpak results used in Li et al., 2018. https://doi.org/10.5194/gmd-2020-256 Preprint. Discussion started: 17 November 2020 c Author(s) 2020. CC BY 4.0 License.

Discussion
Sedapp is a diffusion-based model, and its transport coefficient is a function of both distance from estuary and water depth. Compared with most existing diffusion models based only on water 275 depth, this modification has great advantages in fluvial-deltaic environments, especially for 3D scenarios. Sedapp not only simulates some surface landscapes, but it also reveals some interesting internal features. In the sections beside the channel in Model 1, the formation rate of the clinoforms has close relationship with the distance between the channel and the section. This may be of great significance to the analysis of ancient strata. Considering the resolution of seismic data, it is easier to 280 observe the changes in the density of the foreset than to directly find a channel. This may provide some important supplementary information in areas with less borehole data.
Sedapp also showed strong simulation ability in 2D scenarios. It is not only competent for the shallow sea environment of continental margin, but also competent for the simulation of continental fault basin (Fig.11). The simulation results have strong comparability with previous studies (Li et al., 285 2018). In addition, Sedapp can avoid some potential problems that the water depth models may meet.
The simulation results of Sedapp and water depth models are not very different where the original slope is gentle (Fig.12a, Fig.12b). However, when the slope is steep, the differences are obvious: due to the steep slope and the sharp increase of water depth, the slope break trajectory simulated by the water depth based model increases significantly, even if the sea level remains unchanged at 6m (Fig.12c).

290
This is seriously contrary to the common sense, especially in estuary or delta front environments. In contrast, Sedapp does not face such a problem. As long as the sea level is constant, the slope break line will remain in a straight line and the clinoforms will also move smoothly to the ocean (Fig.12d).
The transport coefficient is a relatively long-term geomorphologic physical quantity, while wave, tidal, and current energy are relatively short-term hydrodynamic quantities. However, they are closely 295 related. A river entering the sea is a type of jet flow phenomenon. The flow velocity decreases rapidly from the river mouth to the sea, which also has a strong negative correlation with the distance to the mouth of the river. The contour map of water flow velocity is fan-shaped. At the same time, the decrease of velocity is also an important cause of sediment deposition, which also explains the close fan-shaped morphology of a delta front. Correspondingly, an increase in water depth will also decrease 300 the flow velocity. For the open coast without river injection, a model based on water depth seems to be https://doi.org/10.5194/gmd-2020-256 Preprint. Discussion started: 17 November 2020 c Author(s) 2020. CC BY 4.0 License.
reasonable. However, for a coast with river injection, it is difficult to explain the formation of the fan-shaped morphology of a delta. Therefore, it can be concluded that, in more general cases, the transport coefficient should be a function of short-term water energy, which is related to both the estuary distance and the water depth. When there is river injection, the river process is dominant and 305 the estuary distance function is a reasonable proxy for the transport coefficient. When there is no river injection, the water depth plays the main role. In addition, the particle size is also one of the decisive factors (Nash 1980;Andrews and Bucknam 1987). Hence, a choice function (see Eq. (9)) and differentiated α's are used to adapt different environments and lithologies. Although the current results of Sedapp seem plausible, these settings for transport coefficient are still empirical. Due to the complex 310 nature of the tranformation from short-term processes to long-term ones, it is difficult to build an accurate bridge between sediment hydrodynamics and stratigraphic formation, while it may be the focus of the next step.

Code availability
The current version of model is available from the project website: The exact version of the model used to produce the results used in this paper is archived on Zenodo.
Input data and scripts of the case studies are also presented in this site. For more details about Sedapp, please contact Jingzhe Li via email lijingzhe@qust.edu.cn.