Evaluating alluvial stratigraphic response to cyclic and non‐cyclic upstream forcing through process‐based alluvial architecture modelling

Formation of alluvial stratigraphy is controlled by autogenic processes that mix their imprints with allogenic forcing. In some alluvial successions, sedimentary cycles have been linked to astronomically‐driven, cyclic climate changes. However, it remains challenging to define how such cyclic allogenic forcing leads to sedimentary cycles when it continuously occurs in concert with autogenic forcing. Accordingly, we evaluate the impact of cyclic and non‐cyclic upstream forcing on alluvial stratigraphy through a process‐based alluvial architecture model, the Karssenberg and Bridge (2008) model (KB08). The KB08 model depicts diffusion‐based sediment transport, erosion and deposition within a network of channel belts and associated floodplains, with river avulsion dependent on lateral floodplain gradient, flood magnitude and frequency, and stochastic components. We find cyclic alluvial stratigraphic patterns to occur when there is cyclicity in the ratio of sediment supply over water discharge (Qs/Qw ratio), in the precondition that the allogenic forcing has sufficiently large amplitudes and long, but not very long, wavelengths, depending on inherent properties of the modelled basin (e.g. basin subsidence, size, and slope). Each alluvial stratigraphic cycle consists of two phases: an aggradation phase characterized by rapid sedimentation due to frequent channel shifting and a non‐deposition phase characterized by channel belt stability and, depending on Qs/Qw amplitudes, incision. Larger Qs/Qw ratio amplitudes contribute to weaker downstream signal shredding by stochastic components in the model. Floodplain topographic differences are found to be compensated by autogenic dynamics at certain compensational timescales in fully autogenic runs, while the presence of allogenic forcing clearly impacts the compensational stacking patterns.


| INTRODUCTION
Alluvial deposition is controlled by both autogenic and allogenic controls (Abels, Kraus, & Gingerich, 2013;Hajek, Heller, & Schur, 2012), which are difficult to be disentangled since they act at overlapping spatial and temporal scales (Bridge, 1993;Stouthamer & Berendsen, 2007). Moreover, it remains to be answered whether autogenic variability can produce cyclicity that resembles allogenic stratigraphic products (Feng et al., 2019;Hajek & Straub, 2017). Allogenic controls on a fluvial depositional environment refer to changes in upstream and downstream conditions, such as climatically driven water discharge and sediment supply variation, tectonically driven basin subsidence and exhumation, and climatically-and/or tectonically driven base-level fluctuation. Some of these controls may have cyclic origins when they are related to Earth's orbital cycles that play important roles in past climate changes. Astronomically driven climate cycles may therefore cause cyclicity in alluvial depositional records (Abels et al., 2013). Autogenic controls refer to internally generated dynamics that can produce complex, relatively large-scale stratigraphic patterns even while allogenic forcing remains constant Hajek & Straub, 2017). Practical examples are channel migration, channel bifurcation and avulsion, crevasse splay formation, levee breaching, consolidation of previously deposited sediments, and interaction with local flora and fauna.
Regularly-distributed paleosols alternating with avulsion belt deposits in the Bighorn Basin, Wyoming, USA, have been ascribed to astronomical climate forcing (Abdul Aziz et al., 2008;Abels et al., 2013Abels et al., , 2016Van der Meulen et al., 2020). Abels et al. (2013) produce a qualitative depositional model for the observed cyclic floodplain successions within the frame of astronomical climate forcing, which is in line with previous sedimentary models (Kraus & Aslan, 1993). In this model, each cycle consists of an avulsion phase featured by thick, laterally-extensive heterolithic intervals and an overbank phase characterized by a few stable channels and fine sediments on which strong palaeosols have developed. These intervals of floodplain stability during the overbank phase alternating with episodes of large-scale fluvial system reorganization during the avulsion phase could be entirely autogenic; however, the remarkable regularity and the close match in timescales with precession forcing indicate these sedimentary cycles to be likely paced by allogenic, astronomically-forced climate changes (Abels et al., 2013(Abels et al., , 2016Van der Meulen et al., 2020). In this case and other similar cases, uncertainties remain, such as how an alluvial system reacts to cyclic upstream climate changes, and in which conditions these climate changes have higher chances of being preserved as sedimentary cycles (Abels et al., 2013). In this context, a process-based numerical alluvial architecture model may help to elucidate the response of an alluvial system to cyclic upstream climate changes through vividly illustrating the interaction between allogenic and autogenic forcing.
Process-based numerical models are suitable for exploring, explaining and predicting river behaviour in response to changing environmental parameters. They are employed to compare the arguably distant relationship between numerical modelling and field observation, namely sedimentary processes and sedimentary products (Hajek & Wolinsky, 2012). Existing basin-scale models for fluvial stratigraphy simulation include Mackey and Bridge (1995), its follow-up Karssenberg and Bridge (2008), Jerolmack and Paola (2007) and Dalman and Weltje (2008), amongst others. Some models that highlight short-scale fluvial dynamics incorporate complex physical processes, such as Delft3D (Lesser, Roelvink, van Kester, & Stelling, 2004). However, these are not yet able to simulate a basin-scale alluvial system at the long timescale we wish to analyse. In this context, we choose the Karssenberg and Bridge (2008; hereafter referred to as KB08) model as an appropriate tool to solve our research questions. The KB08 model is primarily diffusion-based, in which the rate of diffusion is a function of discharge of individual channels and thus varies in space and time. This, together with stochastic components relating to channel bifurcation and avulsion, produces dynamic and unpredictable model behaviour that could be to a large extent comparable to autogenic processes in reality (Karssenberg & Bridge, 2008).
The objectives of this study are four-fold: (a) to elucidate impacts of cyclic and non-cyclic upstream climate forcing, including their amplitude and wavelength, on the modelled alluvial deposition, (b) to investigate conditions that can lead to cyclic stratigraphic records, (c) to gain insights into compensational timescales and mechanisms related to interaction

Highlights
• Cyclicity in Q s /Q w ratio can produce floodplain cyclicity that consists of an aggradation phase and a non-deposition phase. • The non-deposition phase occurs after the Q s / Q w ratio starts to decrease, while the aggradation phase is initiated after the Q s /Q w ratio begins to increase. • Compensational stacking at specific time scales occurs in the fully autogenic scenario, whereas such compensational stacking is strongly interfered when there is allogenic forcing present. • Upstream Q s /Q w signals get destructed and even totally shredded in the downstream propagation process, and only those with large magnitudes and long, but not very long, wavelengths have chances of preservation in the stratigraphic records.

|
between allogenic and autogenic controls, and (d) to illustrate allogenic signal propagation downstream in the alluvial basin.

Karssenberg and Bridge (2008) model
The process-based numerical model presented by Karssenberg and Bridge (2008) is applicable to target basinscale stratigraphy. It considers sediment transport, erosion and deposition within a network of channel belts, floodplain and topographic slope based on the earlier work by Mackey and Bridge (1995). Boundary conditions at the upstream inlet are mean annual water discharge (Q w ) and sediment load (Q s ) that includes both bedload and suspended load. Boundary condition in the downstream is the base level change over the width of the modelled basin. Basin subsidence is absent in the model, but is mimicked by base-level rise.
In the model, sediment and water are transported through a network of channel belts according to the diffusion equation, the rate of which is a function of discharge of individual channels and thus varies in space and time. The width of a channel belt increases over time, while the depth of it depends on water discharge. Channel-belt aggradation and degradation are calculated using the sediment continuity equation. Floodplain aggradation is a function of channel-belt aggradation rate and distance from channel belts, and it occurs only when the topographical elevation of a floodplain cell is less than the that of the nearest channel-belt cell. Floodplain degradation is treated using a sediment diffusion-advection approach. Sediment conservation occurs in the channel belt system, but not in the floodplain, where sufficient sediment is assumed to be available.
There are two ways for channel belt bifurcation to occur in the model. The first is stochastic bifurcation, both location and frequency of which are determined randomly. The second is dependent bifurcation, which is a function of cross-valley slope (S cv ) advantage over the down-valley slope (S dv ) as well as return period of flood discharge exceeding a certain threshold value (Q a ). The probability of a dependent bifurcation, P(B), is calculated at each grid (Equation 1) at each time step: where k s is the slope proportionality constant, e Q is threshold discharge component, and e s is slope exponent (Karssenberg & Bridge, 2008;Mackey & Bridge, 1995).
Whether a bifurcation is initiated or not depends on Equation 2: where P(B) x,y,t is the bifurcation probability at a location with the coordinate of (x, y) at time t, and U (0,1) is a random number chosen from a uniform distribution of 0 to 1.
At bifurcations, water is distributed over two channels according to Equation 3: where q 1 and s 1 are, respectively, the water discharge (m 3 /year) and the gradient of the first bifurcating channel, q 2 and s 2 are the same variables for the second bifurcating channel, and q 0 is the water discharge (m 3 /year) of the channel upstream of the bifurcation point. n 1 1, 2 and n 2 1, 2 generate random variables that have normal distributions with mean zero and variance 2 . These noise terms are included to represent local random effects on channel gradients at the bifurcation point. If q 1 ∕q 0 . or q 2 ∕q 0 falls below a critical value u crit (herein 0.4, Table S1), a bifurcation turns into an avulsion.
Presence of the random number in Equation 2 and noise terms in Equation 3 results in the fact that outputs of KB08 model are not deterministic. Given this, we run 10 realizations for each scenario to statistically confirm our findings, accompanied by full post-processing for all of them.

| Model set-up and scenarios
The modelled rectangular basin has a down-valley length of 60 km and a cross-valley width of 40 km, and it is discretized by square raster cells with a constant cell length of 200 m ( Figure 1). The initial basin slope is set to be 0.00011 (0.11 m/km), and it matches a constant water discharge (Q w0 ) of 7.9 × 10 10 m 3 /year and a constant sediment supply (Q s0 ) of 1.0 × 10 6 m 3 /year to produce graded fluvial profiles (Karssenberg & Bridge, 2008). The total simulation is 40 kyr with a time step of 1 year. Other parameters not specified here are listed in Table S1.

| The equilibrium scenario
The equilibrium scenario is to generate graded fluvial profiles that result from matching water discharge, sediment supply, and basin slope. The base level is set at 0 m at the downstream boundary of the basin throughout the 40-kyr simulation. (1) if P(B) x,y,t > U (0,1) : A bifurcation is initiated; else: No bifurcation is initiated. (3)

| The autogenic scenario
Based on the equilibrium scenario, we create long-term accommodation through steadily rising the base level at a rate of 0.4 m/kyr, which matches the long-term depositional rate that can be approximated by dividing Q s0 (1.0 × 10 6 m 3 /year) over the basin size (60 km by 40 km). In this scenario, stratigraphy is built up under full autogenic controls, since there are no variations in the water discharge and sediment supply at the basin inlet. This scenario will be referred to as Scenario NC (abbreviated for non-cyclic).

| The allogenic scenarios
The relation between water discharge and sediment flux has been widely studied and discussed. In general, water discharge presents a quasi-instantaneous response to precipitation and evaporation changes, despite that multiple catchment parameters have their impacts, such as ground water levels and vegetation cover. The response of sediment flux to water discharge is more complex, and that could be highly non-linear and different at various timescales (Allen, 2008;Coulthard & Van de Wiel, 2013;Forzoni, Storms, Whittaker, & de Jager, 2014;Metivier & Gaudemer, 1999;Tucker & Slingerland, 1997;Whipple & Tucker, 2002). Larger catchments tend to strengthen the inertia and non-linear response of sediment transport systems to external forcing (Allen, 2008;Castelltort & van den Driessche, 2003;Metivier & Gaudemer, 1999). Our objective is to evaluate the effect of cyclic and non-cyclic upstream forcing on alluvial stratigraphy, and we follow a simple relation between sediment load and water discharge based on Syvitski, Morehead, Bahr, and Mulder (2000): where Q st (m 3 /s) is the sediment load at time step t, a is a constant related to long-term water discharge, basin relief and temperature, Q wt (m 3 /s) is the water discharge at time step t, and b is a constant related to long term-suspended load, basin relief and temperature (Syvitski et al., 2000). Constant a and b are adjusted to 9 × 10 −9 and 0.925, respectively, based on sensitivity tests in our model to have sediment and water volumes in balance with basin size, basin slope and base level rise rate. We first design two scenarios, where Q w is in the simple sinusoidal form and varies in amplitude (A = 0.1 or 0.2, Equation 5), and Q s is dependent on Q w following Equation 4.
where A is 0.1 or 0.2, T is 10 kyr, and Q w0 (m 3 /s) is the water discharge at t = 0, which is 7.9 × 10 10 m 3 /year, as has been mentioned above. These two scenarios with Q s /Q w ratio amplitudes of 0.1 and 0.2 will be hereinafter referred to as Scenarios C10 and C20, respectively.
The Q s /Q w ratio is cyclic when Q w is cyclic, as can be clearly seen if we divide Q st on both sides of Equation 4, yielding: As Q s responds exponentially to Q w changes in the relation proposed by Syvitski et al. (2000), the Q s variability will be extremely large when the amplitude of Q w is large. However, the KB08 model is not able to deal with such extremely large Q s variability. Since one of our goals is to test the impact of cyclic Q s /Q w forcing amplitude on the alluvial stratigraphy, we choose to keep Q s stable while varying Q w to result in cyclic Q s /Q w ratios, which shares high consistency with the model choices by Zhang, Burgess, Granjeon, and Steel (2019). Specifically, we design a set of amplitude-themed scenarios, with constant Q w wavelengths of 10 kyr and variable Q w amplitudes from 0.05 to 0.4 (referred to as Scenarios A5, A10, A20 and A40, respectively, Table S2). Similarly, we design wavelength-themed scenarios with constant Q w amplitudes of 20% and varied Q w wavelengths from 1.25 to 20 kyr (referred to as Scenarios W1.25, W2.5, W5, W10, and W20, respectively, Table S2). Meanwhile, we also design a scenario with constant Q w and variable Q s (referred to as Scenario CW, Table S2), which is to validate the role of Q s /Q w ratio.

| Post-processing
All numerical simulations are run on a Linux computer with an Intel i7 eight-core processor at 2.93 GHz using a RAM of 16 GB, and the operation system is Ubuntu 18.04.2 LTS. For each scenario, it takes 10-15 hr to run one 40 kyr-long realization with a basin size of 40 km by 60 km. A series of topographic maps at different time steps are generated, and they are post-processed using Anaconda Spyder (Python 3.6), including calculation of regionally averaged depositional rates, time series analysis of regionally averaged depositional rates, and identification of compensational timescales. These calculations are all implemented in the upstream realm at the crosssection topography along line BB' with its location shown in Figure 1, if not otherwise specified. This decision is based on consideration of the fact that line BB' has a certain distance from the basin inlet, which prevents from overwhelming upstream influence, while it is also relatively far from the basin outlet, minimizing the effect of the base level variation. For all realizations, base-level locations are monitored to ensure that they remain in the downstream portions of the modelled basin.

| Calculation of regionally averaged depositional rates
We calculate the 1D deposition rate (D x,y,t ) at a certain location using Equation 7: where h x,y,t is surface elevation at the location with coordinate (x, y) at time t, Δt. is the sampling interval and D x,y,t is the depositional rate at the location with coordinate (x, y) during the time interval from t to t + Δt.
The average depositional rate over the floodplain width (D x,y,t ) is calculated in order to analyse regional-scale depositional patterns. This is done on cross-section topography along the line BB' (location shown in Figure 1), unless otherwise specified.

| Time series analysis
The REDFIT protocol in the software PAST 3 version 3.14 (Hammer et al., 2001) for time series analysis is deployed on the D x,y,t calculated in Section 2.3.1. The method is run with n = 1,000 Monte Carlo simulations of the autoregressive (AR1) process, oversampling and segmentation intervals of 1, and a Blackman-Harris window. All 90%, 95%, and 99% level χ 2 and the theoretical AR1 are displayed. To focus on long-wavelength signals (>1,000 years), the upper cut-off frequency is set to be 1 kyr -1 in all analyses. Given the good time control in these numerical exercises, we only regard peaks with their confidence levels higher than 99% as reliable signals present in the stratigraphic record (Kemp, 2016;Vaughan, Bailey, & Smith, 2011).

| Compensational timescale identification
We use the standard deviation of sedimentation/subsidence (σ ss ; Wang, Straub, & Hajek, 2011) to characterize the compensational timescale (Equation 8), which is a dimension-reduced term of that by Straub, Paola, Mohrig, Wolinsky, and George (2009): where r (T;x) is the average deposition rate at a horizontal coordinate of x along line BB' during a time interval of T, which ranges from 100 to 40,000 years in this study, L is the cross-basin length, which is 20 km in this study and r (x) is the local long-term sedimentation (or subsidence) rate, herein 0.4 m/kyr. ss is expected to decrease as T increases, following a power-law trend (Equation 9; Straub et al., 2009;Wang et al., 2011): where a ′ is a coefficient, and the exponent, κ, is termed the compensation index. At a certain timescale when κ exceeds 0.5, the stratigraphic stacking transits from random/anti-compensational form to compensational form, and when it reaches 1.0, the stacking is purely compensational ). In other words, a transport system is able to visit every spot in a basin repeatedly during an interval longer than this timescale Wang et al., 2011).

| Characteristics of alluvial stratigraphy in the equilibrium scenario
Over the 40-kyr simulations, the surface topography of Scenario Equilibrium reveals very small changes. Specifically, the cross-section topography varies in a very small range (Figure 2a,b), and the longitudinal topography indicates long-term graded fluvial profiles ( Figure 2c). Correspondingly, slow deposition sporadically occurs at certain locations and certain time moments, accompanied by very small-scale erosion at certain locations and certain time moments (Figure 2d,e). The REDFIT power spectrum shows no peaks with their confidence levels higher than 99% (Figure 2f).

| Characteristics of alluvial stratigraphy in the autogenic scenario
Scenario NC presents long-term topographic build-up (panels a1 and b1, Figure 3), because of long-term accommodation created by base-level rise. Depositional and erosional patterns occur randomly in time and space, with deposition taking the dominance both temporally and spatially (panels c1 and d1, Figure 3). There are no peaks with their confidence levels higher than 99% (panel e1, Figure 3). We apply the method of Straub et al. (2009) in analyzing the cross-section topography (panel b1, Figure 3), and identify a compensational timescale of 2.0 kyr (Figure 4a). Above this compensational timescale, ss decays in the temporal domain following a power law with a nearly identical value of κ.

| Impact of cyclic Q s /Q w ratio
Scenarios C10 and C20 present regularly-alternating episodes of depositional and non-depositional patterns (e.g.  panels a2-d2 and a3-d3, Figure 3), which are not shown in the Scenario NC. Specifically, we refer to the episode with low/negative average depositional rates as the non-deposition phase, as manifested by the near-horizontal average elevation profile (e.g. panel a3, Figure 3), clustered timelines in the cross-section topography (e.g. panel b3, Figure 3), and zero/ negative average depositional rates (e.g. panels c3 and d3, Figure 3). In comparison, we refer to the other episode with positive depositional rates as the aggradation phase, which is embodied by the rapidly increasing average elevation profile (e.g. panel a3, Figure 3), separated timelines in the crosssection topography (e.g. panel b3, Figure 3), and positive and sometimes very high depositional rates (e.g. panels c3 and d3, Figure 3). We observe peaks at periodicities of exact 10 kyr in the power spectra (panels e2 and e3, Figure 3), which coincide with the allogenic forcing wavelength. Using the method of Straub et al. (2009), we found that ss in Scenarios C10 and C20 behaves very differently when compared to Scenario NC. The compensational timescales in Scenarios C10 and C20 are similarly prolonged to ca. 5.1 kyr (Figure 4b,c), if we only look at the σ ss in the < 10 kyr time window. Fluctuation of ss in the > 10 kyr time window presents slightly cyclic features in Scenario C10 ( Figure 4b) and strongly cyclic features in Scenario C20 (Figure 4c), which could be ascribed to the fact that we are estimating compensational timescales across the aggradation/non-deposition phases that comprise allogenic cycles (S. Trampush, pers. comm., June 14, 2019). For Scenario C20, the alluvial system is in a non-deposition phase (Points a to c in Figure 5A) when the Q s /Q w ratio decreases or is low. A few active channel belts are present in the model during this phase, and they constantly widen and gradually build up the topography ( Figure 5B,C). The aggradation phase occurs when the Q s /Q w ratio increases or is high (Points d to f in Figure 5A). This phase is initiated with a first avulsion that results from super-elevation of the channel belt over the floodplain (panel d in Figure 5C). After that, frequent channel avulsion causes rapid aggradation over the entire width of the modelled basin ( Figure 5C). In the aggradation phases, channel belts are narrow due to their short life times, and water discharge in individual channel belts is also relatively low ( Figure 5B).
The KB08 model is stochastic, meaning different outcomes for different realizations of the same scenario. However, based on 10 realizations of Scenario C20, consistency is seen when it comes to regular alternation of aggradation and non-deposition phases ( Figure 6). Their initiation timing lags the (Q s /Q w ) min by 0.6-1.2 kyr in a 10-kyr cycle, while their termination timing lags the (Q s /Q w ) max by 0.3-1.6 kyr in a   10-kyr cycle (Figure 6a). Correspondingly, the relative duration of an aggradation phase is 40%-61% of an allogenic cycle wavelength (Figure 6b). The downstream propagation of the cyclic upstream forcing is analysed via the average depositional rate across the modelled basin in Scenario C20. As the distance from the basin inlet increases, alternation of aggradation and non-deposition phases become less distinct, as reflected by the shortened duration of non-deposition phases (red coloured part in panels a1-e1, Figure 7) and lowered spectra power and confidence levels of the dominant peaks (panels a2-e2, Figure 7). No >99% confidence level peaks are observed at cross-section along line DD' (panel d2, Figure 7), though the most dominant peak is still at a periodicity of 10 kyr. In comparison, no peaks are present at the allogenic wavelength at cross-section along line EE' (panel e2, Figure 7), which indicates full shredding of the upstream Q s /Q w signal. The location where the base level reaches the alluvial plain is shown in the Figure S1.

| Impact of Q s /Q w ratio amplitude
Scenario A5 presents similar depositional patterns to Scenario NC that has a zero Q s /Q w ratio amplitude (panels a1-d1, Figure 8), and there are no peaks with >95% confidence levels at the forcing wavelength in the power spectrum (panel e1, Figure 8). In comparison, cyclic patterns start to emerge in Scenario A10, reflected in the elevation profile, cross-section topography, chronostratigraphic diagram, and average depositional/erosional rate (panels a2-d2, Figure 8), accompanied by a >99% confidence level peak at 10 kyr that coincides with the allogenic forcing wavelength (panel e2, Figure 8). Scenarios A20 and A40 present clearer alternations of aggradation and non-deposition phases (panels a3-d3 and a4-d4, Figure 8), accompanied by peaks at 10 kyr with higher power and confidence levels (panels e3 and e4, Figure 8). Time lines for the non-deposition phase in cross-section topography of Scenario A40 are extremely clustered, indicating very low sedimentation over a long time period, associated with greatly confined channel belts (e.g. red eclipse marked in panel b4, Figure 8) that are later filled (e.g. red eclipse marked in panel c4, Figure 8). Note that the aggradation phase is defined as the moment in time when channel belt avulsion and rapid aggradation occurs on the alluvial plain. Aggradation within the confined channel belts, until the moment when super-elevation is reached and avulsion is triggered, is also defined as part of the non-deposition phase. Scenario CW presents highly similar features to Scenarios C20 and A20 ( Figure S2), and they share very similar Q s /Q w ratio amplitudes and exactly the same Q s /Q w ratio wavelengths.

| Impact of Q s /Q w ratio wavelength
With a Q s /Q w ratio amplitude of 20%, even a very small Q s / Q w ratio wavelength can lead to clear cyclicity in the temporal domain, as reflected in the chronostratigraphic diagram (e.g. panels c1 and 2, Figure 9), the deposition rate time series (e.g. panels d1 and d2, Figure 9), and the >99% confidence level peaks at the forcing wavelength in the power spectra (e.g. panels e1 and e2, Figure 9). However, many cycles with their wavelengths shorter than 5 kyr are overprinted by subsequent cycles (e.g. panels b1 and b2, Figure 9). Interestingly, some

F I G U R E 4 Decay of σ ss with increasing time window for
Scenarios NC (a), C10 (b), and C20 (c). Error bars represent geometric standard deviation of σ ss , while red dots indicate the average σ ss at corresponding time window, blue dashed trend lines represent best-fit linear regression to log-log data, and vertical blue dashed lines, passing through the intersection points of two blue trend lines, indicate the predicted compensational time scales, over which the the stratigraphic stacking transits from random/anti-compensational to compensational form seemingly regular, large-scale aggradation occurs at a longer wavelength than the forcing wavelength in Scenario W1.25 (panels a1 and d1, Figure 9). In other words, there seem to be 4-5 larger-scale cycles in Scenario W1.25, besides 16 smaller-scale cycles with wavelengths of 1.25 kyr.
As the allogenic forcing wavelength increases from 1.25 to 10 kyr, the aggradation and non-deposition phases become more clearly separated topographically, reflected in the average elevation profile, cross-section topography and chronostratigraphic diagram (Figure 9). However, when the allogenic forcing wavelength is as long as 20 kyr, some regional aggradation occurs during its non-deposition phase (e.g. the peak at around 19 kyr in panel d5, Figure 9), and such aggradation is smaller magnitudes than its counterpart in the aggradation phase. This relatively smaller-scale aggradation occurring in the non-deposition phase makes the resultant depositional patterns less obviously cyclic, despite the fact that the spectral power peak at the forcing wavelength crosses the 99% confidence level (panel e5, Figure 9).

| KB08 model as a tool for alluvial stratigraphy simulation
Numerical modelling of alluvial stratigraphy has its clear limits when it comes to representation of true stratigraphic records (Hajek & Wolinsky, 2012). However, it is very useful to test a system's responses to changing boundary conditions, where outcrop analogue studies often lack data and resolution. Specifically, from one field case to another, multiple boundary conditions change, and it is often difficult for outcrop analogue studies to identify the dominant processes that have controlled the observed changes in the alluvial stratigraphy. In contrast, we can observe responses of the system to changing boundary conditions through numerical modelling. Nevertheless, it is crucial to evaluate to what extent the numerical modelling results represent true alluvial processes and stratigraphic responses.
It is challenging for a single modelling tool to incorporate full temporal and spatial scales that autogenic and allogenic forcings operate over. The KB08 model employs a diffusion scheme to illustrate the stratigraphic patterns at a basin-filling timescale, which is much longer than a Blue and red dashed lines indicate (Q s /Q w ) min and (Q s /Q w ) max , respectively. Box and whiskers plots show the range of data (whiskers), the lower and upper quartiles (box boundaries), and median (orange line within the box). Isolated dots represent extreme outliers, defined as greater than two times the interquartile range below or above the first or third quartiles, respectively. Numbers in formats of a ± b represent "mean ± standard deviation" of the time lag of an aggradation phase initiation or termination with respect to (Q s /Q w ) min or (Q s /Q w ) max . (b) Relative duration of aggradation phases. Numbers in formats of a ± b represent "mean ± standard deviation" of the relative aggradation phase duration in one cycle  morphodynamic timescale (Paola, Straub, Mohrig, & Reinhardt, 2009). Therefore, the KB08 model likely underrepresents the stochasticity of the transport system, and the model tends to smoothen topography with the aim of algorithmic simplicity and computational economy. Yet, most morphodynamics occur at a channel scale, where the channels define the channel belt. In this study, we only focus on the channel-belt scale, assuming that intra-channel-belt morphodynamics will not significantly change the channel-belt morphodynamics. Modelling tools with layer-averaged Naiver-Strokes equations coupled to non-linear sediment transport algorithms, such as Delft3D (Lesser et al., 2004), are superior in simulating shorter-time-scale morphodynamics, but they are not yet computationally feasible for simulation over basin-filling timescales that we aim to study. In terms of our scenario set-up about Q w and Q s forms, we always adopt a sinusoidal form of Q w , based on the previous research that has shown that cyclic orbital forcing can result in cyclic precipitation variability, which further induces a similar time evolution of run-off that transmits the climate signal to the geomorphic system (Godard, Tucker, Fisher, Burbank, & Bookhagen, 2013). For Q s forms, Q s signal is dependent on Q w in Scenarios C10 and C20, following a non-linear rule, and this is consistent with the diffusive transport model (e.g. Whipple & Tucker, 2002), which assumes sediment flux to be a function of water discharge and river slope (Armitage, Whittaker, Zakari, & Campforts, 2018; F I G U R E 8 Alluvial stratigraphic patterns of scenarios with different Q s /Q w amplitudes (Scenarios A5, A10, A20, and A40). Curves of Q s /Q w ratios are displayed above panel a of each scenario, and a in the y axis equals 0.8 × 10 −5 and b equals 1.8 × 10 −5 . Red ellipse in panel b4 shows an example of largely confined channel belts, while red ellipse in panel c4 shows an example of rapid aggradation in the channel belt before superelevation is reached and avulsion occurs. For other legends, see details in Figure 2   Forzoni et al., 2014;Nittrouer, Shaw, Lamb, & Mohrig, 2012;Syvitski et al., 2000;Whipple, 2001). In amplitudeand wavelength-themed scenarios, Q s signal is independent on Q w , and this is partially in line with the advective stream power model (e.g. Willett, McCoy, Perron, Goren, & Chen, 2014, among many others), in which sediment is assumed to be supply-or detachment-limited. Specifically, the relationship between available sediment supply and the river transport capacity greatly influences the Q s signal.
When the available sediment amount is always smaller than the lowest river transport capacity, Q s should match the erosion rate estimates that balance estimates of the long-term rock uplift rate (Armitage et al., 2018;Godard et al., 2013;Whipple, 2001;Zhang et al., 2019). However, no essential difference is observed among these scenarios with and without Q s dependence on Q w in our modelling results, in terms of regular alternation of aggradation and non-deposition phases, when they have similar Q s /Q w amplitudes and wavelengths F I G U R E 9 Alluvial stratigraphic patterns of scenarios with different Q s /Q w wavelengths (Scenarios W1.25, W2.5, W5, W10, and W20).
Curves of Q s /Q w ratios are displayed above panel a of each scenario, and a in the y axis equals 1.0 × 10 −5 and b equals 1.5 × 10 −5 . For other legends, see details in Figure 2   (Figures 3 and 8; Figure S2). Therefore, all our interpretations and following discussions are targeting the role of the Q s /Q w ratio. In the KB08 model, different realizations of one scenario with identical boundary conditions show different outcomes. These variabilities are ascribed to the presence of stochastic components embedded in the model, which produce unpredictable model behaviours mimicking effects of autogenic forcing. Similarity between allogenic realizations with the same inputs is, however, robust at the allogenic forcing wavelength scale, in terms of regular alternation of aggradation and non-deposition phases. Therefore, we think the KB08 model helps elucidate the interaction between allogenic and autogenic forcing and thus serves well for the scope of our study.

Q s /Q w signal on modelled alluvial stratigraphy
Cyclic Q s /Q w forcing can drive the alluvial deposition through interaction with autogenic processes. An increase in Q s /Q w ratio results in immediate aggradation in the channel belt, and, depending on preceding floodplain topography, channel-belt avulsions and aggradation over the width of the modelling basin will occur when the superelevation is reached. These aggradation and non-deposition phases are in pace with the imposed Q s /Q w ratio changes with sufficiently large amplitudes and long wavelengths. In scenarios with short Q s /Q w wavelengths, we observe longer-time-scale basin-wide aggradation in addition to the shorter-time-scale aggradation as responses to Q s /Q w changes (panels a1, b1 and d1, Figure 9). In other words, there seem to be four longerwavelength cycles in addition to the 32 short-wavelength (1.25 kyr) cycles that are expected based on the imposed Q s /Q w forcing. However, this longer-time-scale basin-wide aggradation does not occur in the fully autogenic scenario, which means that this behaviour is not intrinsic to the model creating superelevation of the channel belts over the alluvial plain at quasi-fixed timescales. It instead means that these long-term phases are a response to the imposed allogenic forcing. Our model runs are not long enough to statistically determine any stable frequency in these long-term basin-wide aggradation phases. It is beyond the scope of the current work to analyse this model character any further, which would need very long model runs with a range of high frequency forcing. For scenarios with long wavelengths, a similar model behaviour is observed. Basin-wide avulsion and aggradation occurs not only in pace with the forcing, but also intermittently. Specifically, superelevation can be reached over a sufficiently long time period in the non-deposition phase of a long allogenic cycle, given the fact that the depositional rate in this phase is very low. Similarly, it takes shorter time periods for autogenically-controlled aggradation in the nondeposition phases to occur in small-amplitude scenarios (e.g. panels c1 and c2, Figure 8), due to the fact that channel belts are easier to be filled and superelevated by within-channel aggradation in these scenarios. Therefore, we can conclude that, the model prefers basin-wide aggradation to occur under autogenic controls at a certain timescale that is long enough for channel-belt aggradation to reach the superelevation. The exact timescale at which this autogenic behaviour happens depends obviously on the allogenic forcing but also inherent model properties (e.g. basin subsidence rate, basin size, and basin slope).

| Preservation and shredding of allogenic Q s /Q w signals
Preservation of allogenic signals in fluvial records requires either the allogenic signal to be relatively strong or the autogenic forcing to be relatively weak. The autogenic forcing presents a compensational timescale of ~2 kyr (Figure 4a), which act as a low bandpass filter to prevent the allogenic signal from being preserved, known as the signal shredding effect (Hajek & Straub, 2017;Jerolmack & Paola, 2010;Straub & Foreman, 2018;Toby, Duller, De Angelis, & Straub, 2019). As demonstrated in our KB08 model results, Q s /Q w ratios with larger amplitudes can better withstand signal shredding from autogenic dynamics. Based on Scenarios W1.25 and W2.5 (Figure 9), Q s /Q w ratios need to have sufficiently long wavelengths to prevent a cycle from being overprinted by the subsequent cycle. This is consistent with the argument about the maximum surface roughness proposed by Wang et al. (2011), which highlights that the stratigraphic package should be thicker than the maximum channel depth in order to be preserved. Meanwhile, the Q s /Q w wavelength should not be too long, which will render the non-deposition phase subjected to autogenically controlled aggradation (panels a5 and c5, Figure 9). In other words, although the depositional rate in the channel belt during a non-deposition phase is very low, super-elevation can be reached if the time is sufficiently long. The exact amplitude and wavelength thresholds between signal preservation and shredding again depends on the specific model set-up, including basin size (Powell, Kim, & Muto, 2012), basin subsidence rate, bank cohesiveness and also measurement locality. As determined by the long-term accommodation, the long-term aggradation rates to a large extent also determine whether signals will be preserved or shredded (Foreman & Straub, 2017 depict amplitude and wavelength thresholds between signal preservation and shredding for our model and scenario setup ( Figure 10). In general, allogenic signals with large amplitudes and long, but not very long, wavelengths have large chances to be preserved in the alluvial stratigraphy.
However, it is worthwhile mentioning that these threshold values in Figure 10 are relative and they are not universally applicable to another modelled/real basin that has different properties. Moreover, there are many ideal situations in this study that potentially contribute to higher chances of signal preservation in this numerical exercise than those in reality. For instance, there is a perfect balance between average Q s /Q w ratio and basin slope and that between sediment supply rate (1.0 × 10 6 m 3 /year), basin size (60 km by 40 km) and base-level rise rate (0.4 m/kyr). Moreover, the constantly rising base level provides continuously growing accommodation that would favour longterm aggradation, which also enhance the chances of signal preservation. Furthermore, we have to mention that it is relatively more feasible to detect signals out of the modelled stratigraphy than the case in reality, because we have great control of time in the modelled basin and thus spectral analysis has an extremely high accuracy. The absent mass balance for floodplain aggradation plays a certain role in helping signal preservation, but its effect should be minor. Overall, it should be harder for field geologists to detect allogenic signals in the real basins due to complex natural conditions that might not synergy with each other in preserving allogenic signals.

| Compensational timescales and autogenic forcing
The compensational timescale represents the timescale at which the floodplain roughness is levelled out over a certain segment of a sedimentary basin ). There are many factors that can influence the compensational timescale, such as basin size, basin slope, magnitudes of sediment load (Q w0 in our study) and water discharge (Q w0 in our study ;Jerolmack & Paola, 2010;Powell et al., 2012;Straub & Esposito, 2013). As stated above, the compensational timescale in a fully autogenic run can be very different from that in an allogenic run, due to interference of allogenic forcing on the autogenic behaviour of river networks. For instance, during a non-deposition phase in an allogenic scenario that lasts longer than the compensational timescale in a fully autogenic scenario, channels are strictly confined within channel belts, which suppresses the inherent tendency of a transport system to level up topographic differences throughout the basin. This is consistent with the argument by Abels et al. (2013) that super-elevation during the overbank phase with low depositional rates will only be reached at sufficiently long timescales such that cyclic climate forcing may pace the changes between overbank phases and avulsion phases.
In general, a large forcing amplitude will suppress the tendency of a transport system to move laterally during the non-deposition phase, as the higher stream power with respect to sediment load causes confined, incisive channel belts (e.g. panel b4, Figure 8). How exactly the allogenic forcing with varied amplitudes and wavelengths interferes with autogenic forcing in producing different compensational timescales is beyond the scope of the current work, and it will be carefully evaluated in the future research.

| Model-field comparison
In the KB08 model, Q w and Q s are set to represent the annual mean discharge and annual sediment load, respectively, and their ratio controls the fluvial graded profile, which is consistent with previous research (Montgomery, 2001;Van den Berg van Saparoea & Postma, 2008;Willett & Brandon, 2002). Increasing Q s /Q w ratio, either resulting from relatively decreasing water discharge to sediment load or relatively increasing sediment load to water discharge, contributes to overall regional-scale sedimentation in the model. In contrast, decreasing Q s /Q w ratio results in highly organized channel belts eroding into the alluvial plain and bypassing sediments, as can be clearly seen in Scenario A40 (panel b4, Figure 8). This is in line with phenomena as reported by previous studies (Blum & Aslan, 2006;Demko, Currie, & Nicoll, 2004). Deposition occurs in the model after the Q s /Q w ratio begins to increase, first filling the topographically low areas and then creating super-elevation (e.g. Figure 5C). With a time lag with respect to (Q s /Q w ) min , the aggradation phase takes F I G U R E 1 0 Amplitude and wavelength thresholds for Q s /Q w signal preservation. Hollow stars indicate shredded signals, while solid ones suggest preserved signals. Note that these threshold values are dependent on measurement localities and inherent properties of the modelled basin, and they can not be universally applied to all localities in all basins place when lateral topographic gradients increase and superelevation is sufficiently high. This is also in line with previous field studies, such as Mohrig, Heller, Paola, and Lyons (2000) and Slingerland and Smith (2004). In the KB08 model, during the time interval with an decreasing Q s /Q w ratio, non-deposition phases occur on the floodplains with stable channel belts that only locally avulse or splay. Moreover, a longer non-deposition phase, in the field context, contributes to generally more developed soil profiles. The modelled non-deposition phase seems to share similarities with the so-called overbank phase in the qualitative precession-forcing model of Abels et al. (2013) for the lower Eocene series of the Bighorn Basin, Wyoming, USA. The modelled non-deposition phase is characterized by relatively stable channel-belt locations and insignificant aggradation on the floodplain, and local avulsion may still occur. In the precession-forcing model, the overbank phase shows strongly developed paleosols on floodplain fines that indicate long periods of non-deposition or low deposition (Abels et al., 2013). Field data concerning these floodplain cycles have so far only been gathered in floodplain strata. Therefore, it is not known yet whether these overbank phases coincide with the large-scale incision of the channel belt. Most overbank phases on the floodplains of the Bighorn Basin are however characterized by multiple soils that may be cumulative in nature (Kraus, 2002;Kraus & Gwinn, 1997). This means that small amounts of sediment arrive at these distal floodplain positions, also when channel belts have stable positions and possibly only splay and avulse locally.
In the KB08 model, during the time interval with increasing Q s /Q w ratio, sediment filling first takes place in the modelled channel belts, gradually leading to super-elevation and consequent frequent avulsion due to high depositional rates. The modelled aggradation phase resembles the so-called avulsion phase in the precession-forcing model of Abels et al. (2013). In the precession-forcing model, quickly avulsing channel belts at regional scale produce extensive sheet-like heterolithic avulsion-belt deposits (Abels et al., 2013), filling topographic differences between super-elevated channel belts and distal floodplains. In the KB08 model, the floodplain build-up continues as long as the Q s /Q w ratio keeps increasing; stratigraphic build-up with avulsion-belt deposition continues also when super-elevation is levelled out already.
The duration of avulsion-belt formation has been a topic of investigation, as palaeo-climatic data are needed to be set in sub-precession age models (e.g. Bowen et al., 2015;Wang et al., 2017). However, there are no ways to produce age models at a precision within a precession cycle for the Eocene or anywhere in the pre-Quaternary. Constraints on the relative duration of overbank and avulsion phases within the 21-kyr precession cycle remain dependent on present-day estimates of avulsion-belt deposition rates (Bowen et al., 2015). Wang et al. (2017) argued that the avulsion-belt deposition at one single location represents only 5%-10% of the precession cycle. According to Bowen et al. (2015), to produce the approximate 7-m-thick avulsion belt deposits, the duration of avulsion-belt formation should last between 3.5 and 14 kyr, that means, the avulsion duration accounts for 18%-70% of a precession cycle. In our model, the duration of an aggradation phase is negatively correlated with Q s /Q w forcing amplitude, which varies between 40% and 61% of a cycle wavelength for Scenario C20 (Figure 6b). Whether our results provide useful ways to constrain age models for ancient floodplain stratigraphy that are also used as paleoclimatic records (Abels, 2012;Bowen et al., 2015) remains to be resolved by more detailed 3D characterization of field data.

| CONCLUSIONS
In this study, we evaluate the effect of cyclic and noncyclic upstream forcing on alluvial stratigraphy using a process-based alluvial architecture model (Karssenberg & Bridge, 2008). We find that cyclicity is preserved in the sedimentary records when the sediment load over water discharge (Q s /Q w ) ratio has sufficiently large amplitudes and long, but not very long, wavelengths, the absolute values of which depend on inherent properties of the modelled basin. Within one allogenic cycle, there are a non-deposition phase and an aggradation phase. After the Q s /Q w ratio starts to decrease, the non-deposition phase occurs, during which channel belts are confined and stable at certain locations, and the vast floodplain undergoes very low sedimentation. The aggradation phase occurs after the Q s /Q w ratio starts to increase, accompanied by sedimentary filling in the channel belt and creation of channel belt superelevation over the adjacent floodplains. Frequently shifting channel belts cause rapid sedimentation over the entire basin. The non-deposition and aggradation phases are comparable to the overbank-and avulsion-phases as identified in field successions in the Bighorn Basin, Wyoming, USA, which are attributed to precession-paced climate changes. In the KB08 model, larger forcing amplitudes result in longer non-deposition phases due to the fact that longer time is needed to fill the deeper-incised channel belts before superelevation can be reached. The upstream Q s /Q w signal is shredded in the downstream transmission process, and large amplitudes are more favourable for Q s /Q w signal preservation. We identify compensational timescales in the fully autogenic model runs by applying the method of Straub et al. (2009), and we find that the presence of an allogenic forcing will strongly influence the compensational timescale of a transport system. Findings of this study provide insights into the transmission and preservation processes of upstream cyclic Q s /Q w signals, during which allogenic and autogenic forcings interact with each other to produce alluvial stratigraphy.