The impact of clastic syn‐sedimentary compaction on fluvial‐dominated delta morphodynamics

In natural deltaic settings, mixed hydrodynamic forcings and sediment properties are known to influence the preserved delta deposits. One process that has not received much attention yet is syn‐sedimentary compaction of clastic sediment on millennial‐scale delta evolution. To study how compaction interacts with delta morphodynamics and preserved sediment, a modelling approach is proposed. A 1D grain‐size dependent compaction model was implemented into Delft3D‐FLOW, which provides an opportunity to understand the underexplored connection between grain sizes supplied to the deltas and sediment compaction. The compaction model allows deposited sediment to decrease in volume due to the accumulation of newly deposited sediments above or the elapsed time. Differences in morphological trends are presented for scenarios defined by the composition of sediment supply (mud rich and sand rich) and the maximum allowed compaction rate in the model (0–10 mm year−1). The resultant deposits are classified into sub‐environments: delta top, delta front and pro delta. The delta top geometry (e.g. area increase, rugosity and aspect ratio), sediment distribution alongshore and across sub‐environments, and delta top accommodation (e.g. volume reduction and average water depth) are compared. The modelling results show that compaction of the underlying delta front and pro delta deposits increases the average water depth at the delta top, driving morphological variability observed in the mud‐rich and sand‐rich deltas. The morphological changes are more prominent in the mud‐rich deltas, which experience larger compaction‐induced volume reduction for the same scenario. Moreover, higher compaction rates further increase the delta top accommodation, resulting in more deposition and evenly distributed sediment at the delta top. This leads to a less significant area increase and a wider delta top with a smoother coastline. The presented modelling results bridge the knowledge gap on the influence of syn‐sedimentary compaction on long‐term delta morphodynamics and preserved sediment. These findings can be applied to unravel the controlling processes in ancient delta deposits and predict the evolution of modern systems under changing climates.


Abstract
In natural deltaic settings, mixed hydrodynamic forcings and sediment properties are known to influence the preserved delta deposits. One process that has not received much attention yet is syn-sedimentary compaction of clastic sediment on millennial-scale delta evolution. To study how compaction interacts with delta morphodynamics and preserved sediment, a modelling approach is proposed.
A 1D grain-size dependent compaction model was implemented into Delft3D-FLOW, which provides an opportunity to understand the underexplored connection between grain sizes supplied to the deltas and sediment compaction. The compaction model allows deposited sediment to decrease in volume due to the accumulation of newly deposited sediments above or the elapsed time. Differences in morphological trends are presented for scenarios defined by the composition of sediment supply (mud rich and sand rich) and the maximum allowed compaction rate in the model (0-10 mm year −1 ). The resultant deposits are classified into sub-environments: delta top, delta front and pro delta. The delta top geometry (e.g. area increase, rugosity and aspect ratio), sediment distribution alongshore and across sub-environments, and delta top accommodation (e.g. volume reduction and average water depth) are compared. The modelling results show that compaction of the underlying delta front and pro delta deposits increases the average water depth at the delta top, driving morphological variability observed in the mud-rich and sand-rich deltas. The morphological changes are more prominent in the mud-rich deltas, which experience larger compaction-induced volume reduction for the same scenario. Moreover, higher compaction rates further increase the delta top accommodation, resulting in more deposition and evenly distributed sediment at the delta top. This leads to a less significant area increase and a wider delta top with a smoother coastline. The presented modelling results bridge the knowledge gap on the influence of syn-sedimentary compaction on long-term delta morphodynamics and preserved sediment. These findings can be applied to unravel the controlling processes in ancient delta deposits and predict the evolution of modern systems under changing climates.

| INTRODUCTION
Compaction occurs in all deltas during and after their evolution, causing volume reduction of the deposited sediments. The leading cause of compaction is the expulsion of over-pressured pore fluid induced by increasing overburden weight due to ongoing deposition (von Terzaghi, 1943). This dewatering process is syn-sedimentary, occurring during active deposition. Compaction also occurs postdepositionally under constant overburden stress, caused by time-dependent pore fluid expulsion in fine-grained sediments (van Asselen et al., 2009;Zoccarato et al., 2018).
Delta evolution is thought to be controlled by hydrodynamic forcings and sediment properties (Galloway, 1975;Orton & Reading, 1993). As these controls differ between deltas, many studies have been conducted to decipher their process interactions using the preserved delta deposits (Bhattacharya & Willis, 2001;Goodbred et al., 2003;Rice et al., 2020). However, delta deposits are not straightforward to interpret, which is further complicated by autogenic processes in the delta, such as syn-sedimentary compaction. Very little is known about the impact of synsedimentary compaction on millennial-scale delta evolution, mainly due to a lack of field measurements and limited possibilities to acquire such measurements considering the spatial (> 10 3 km 2 ) and temporal (> 10 3 years) scale of delta evolution. Compaction datasets from contemporary deltas are available in the literature (Aly et al., 2012;Becker & Sultan, 2009;Gebremichael et al., 2018;Higgins et al., 2014;Minderhoud et al., 2018;Saleh & Becker, 2018;Stanley, 1990;Teatini et al., 2011). However, sedimentation is low in the measurement sites due to reduced sediment supply imposed by dams and artificial levee construction. Therefore, these compaction measurements are most likely post-depositional rather than syn-sedimentary. Compaction datasets are also available from laboratory experiments, which provide a first-order approximation of compaction behaviour (Bjerrum, 1967;Merckelbach & Kranenburg, 2004;Mesri, 2003;Mesri et al., 1997;Taylor & Merchant, 1940;Toorman, 1996Toorman, , 1999Townsend & McVay, 1990;von Terzaghi, 1943;Winterwerp & van Kesteren, 2004). Despite some shortcomings, these are the most reliable published compaction datasets.
Syn-sedimentary compaction occurs in active depositional areas where the freshly deposited sediments dewater and compact by subsequent deposition. This creates additional accommodation to store more sediments (Colombera & Mountney, 2020;Paumard et al., 2020;van Asselen et al., 2011), leading to a feedback loop between sedimentation and accommodation. This shows the potential impact of compaction on delta morphodynamics and the sediment archive, which is still underexplored. Therefore, It is important to include compaction in the simulation of deltas to better understand its influence on natural delta evolution.
To this end, this paper has two objectives: (1) formulate, implement and evaluate a new compaction algorithm into the existing open-source coupled hydrodynamicmorphodynamic model (Delft3D) (Lesser et al., 2004), and use the updated code to (2) demonstrate how compaction interacts with syn-sedimentary delta evolution by analysing delta top morphology (i.e. area increase, rugosity and aspect ratio), sediment distribution alongshore and across sub-environments (i.e. delta top, delta front and pro delta), and delta top accommodation.
This work builds on previous research using Delft3D that simulates delta development and characterises the impact of grain size on delta evolution (Burpee et al., 2015;Caldwell & Edmonds, 2014;Edmonds & Slingerland, 2007, 2008, 2010Geleynse et al., 2010Geleynse et al., , 2011Hillen et al., 2014;van der Vegt et al., 2016van der Vegt et al., , 2020. Previous compaction studies using Delft3D, which focus on a shorter timescale model evolution, were used to guide the implementation of the compaction algorithm into Delft3D (Winterwerp et al., 2018;Zhou et al., 2016). The modified Delft3D was used to generate progradational and aggradational fluvial-dominated deltas. These datasets were then analysed to gain insight into the role of compaction on delta morphodynamics.

| METHODOLOGY
Laboratory experiments show that compaction occurs in two phases: primary and secondary compaction (Bjerrum, 1967;Mesri, 2003;Mesri et al., 1997;Taylor & Merchant, 1940;von Terzaghi, 1943). The primary compaction is a syn-sedimentary process, imposing volume reduction due to the weight of depositing sediment. This phase is characterised by a large magnitude of volume reduction that occurs over a short timescale. Subsequently, the secondary compaction continues post-depositionally with a lower magnitude over a longer timescale.

K E Y W O R D S
accommodation, delta morphology, preserved sediment, process-based forward models, synsedimentary compaction Once overburden pressure exceeds pore pressure, the pore fluid is expelled out of the pores, ultimately flowing to the surface as a groundwater flow. This means that local porosity and permeability will interfere with the actual compaction rates of the sediments. If the fluid cannot escape the pores, an overpressure condition will occur. This study does not incorporate complex 3D groundwater flow and assumes that the expelled pore water can always escape vertically to the surface (Meckel et al., 2006(Meckel et al., , 2007. Furthermore, it was assumed that different pore fluid types (i.e. air and water) have a negligible impact on compaction over geological timescales (>10 3 years). Therefore, there is no differentiation between compaction and consolidation, which is important in geotechnical applications.
For this purpose, new 1D primary and secondary compaction formulas were developed for clastic sediments (e.g. sand and mud) based on methods described by Terzaghi and Mesri (Mesri, 2003;Verruijt, 2010;von Terzaghi, 1943), which exclude the permeability variation of sediments. Subsequently, the compaction formulas were implemented into Delft3D to simulate compaction's impact on seaward building deltas that actively aggrade. This is not the first time compaction has been incorporated into a process-based numerical simulation. Previous studies have used a simpler model with uniform compaction that occurs regionally to study delta evolution (Liang et al., 2016a(Liang et al., , 2016b. It is important to note that the secondary compaction is included in the modelling because this compaction type operates over a long timescale (Bjerrum, 1967), comparable to the timescale over which the simulated deltas evolve (>10 3 years). However, the primary compaction is still the dominant mechanism when the simulated deltas are actively depositing, triggering compaction due to the increased weight of the deposited sediment. Below, the compaction formulas are described together with their implementation into Deft3D FLOW2D3D version 6.02.08.62644 (Deltares, 2021a).

| Primary compaction
von Terzaghi (1943) shows that the subsidence of a sediment bed (ΔH p in m) due to primary compaction depends on the inverse of Young's modulus for the sediment ( 1 E in kg −1 m s −2 ) and the effective stress imposed on the bed, which is the difference between overburden stress ( in kg m −1 s −2 ) and the opposing pore fluid pressure (p in kg m −1 s −2 ), as shown in Equation 1.
It is assumed that the sediment grains and pore fluid are incompressible and insoluble. Therefore, the subsidence only occurs by porosity reduction due to pore fluid expulsion. Given a sediment bed with thickness H (m), the subsidence can be estimated using Equations 2 and 3: For geological timescales simulation, detailed modelling of pore fluid expulsion rate, indicated by the second term on the right-hand side of Equation 3, is simplified by using a constant rate C rp (m year −1 ), excluding the effects of the sediment's permeability variation by assuming the porosity is uniformly connected. Consequently, pore fluid can always escape during compaction, which allows compaction to be modelled with less complexity than in previous studies (Winterwerp et al., 2018;Zhou et al., 2016). The simplified Equation 3 is shown below, which was implemented into Delft3D.
where Δ H p,t is the subsidence due to primary compaction at each simulation time interval, C rp the primary compaction rate, H t−1 represents the bed thickness at the previous simulation time interval, and t − t−1 indicates the net increase in overburden stress between two successive simulation time intervals, which is positive during deposition. However, if the net increase in overburden stress is zero or negative due to hiatus or erosion, respectively, the primary compaction that occurs previously is deactivated, whereas the secondary compaction is activated because only one compaction type is active at each simulation time interval. The primary compaction is applied to coarse-grained and fine-grained sediments.

| Secondary compaction
After the pore pressure set up by the overburden stress has dissipated due to pore fluid expulsion, the secondary compaction continues at a slower rate (Mesri et al., 1997). This compaction type occurs for an extended period (on the order of thousands of years) (Bjerrum, 1967). The secondary compaction was formulated as a function of time, assuming an independent relationship with overburden stress, as shown by Equation 5. (1) Here, the subsidence imposed by secondary compaction (ΔH s,t in m) depends on the bed thickness at the previous simulation time interval (H t−1 in m), the bed porosity at the previous simulation time interval (e t−1 ), secondary compaction rate (C rs in m year −1 ) and elapsed time relative to a reference time ( t t o ). The secondary compaction occurs in finegrained sediments due to their low hydraulic permeability, resulting in delayed pore fluid dissipation Zoccarato et al., 2018). The secondary compaction only occurs after the onset of the primary compaction. Therefore, if the deposited sediment never experiences any overburden to turn the primary compaction on, the sediment will never experience the time-dependent dewatering (secondary compaction). The compaction formulas (Equations 4 and 5) are based on laboratory studies measured over a short timescale (on the order of days). To model delta evolution, which occurs over geological timescales (>10 3 years), the primary and secondary compaction rate (C rp and C rs ) was adjusted to the compaction rate of natural deltas, which is explained in more detail in Section 3.2. In addition, the morphological changes are accelerated by upscaling the modelling time (t and t 0 ) using an acceleration factor called the morphological scaling factor (MORFAC) (see also Section 2.1). Therefore, the resultant subsidence due to primary and secondary compaction (ΔH p and Δ H s ) is over geological time.
It is important to note that the actual compaction rate varies locally in the model due to spatial and temporal variations of sedimentation and erosion. This leads to complex feedback between compaction and morphodynamics and preserved sediment. The compaction formulas (Equations 4 and 5) are validated with measurement data in Text S1.

| Implementation of compaction formula into Delft3D
In Delft3D, sediments smaller than or equal to 64 μm are assumed cohesive and transported as suspended load. Sediments larger than 64 μm are partially transported in suspension, adding to the suspended load, and partially through saltation and creep, constituting bedload. In this study, the cohesive suspended load is transported using a cohesive transport formulation, calculated by solving a depth-averaged (2DH) advection-diffusion (mass balance) equation for the suspended sediment (Galappatti, 1983). The non-cohesive suspended and bedload are transported using the Engelund-Hansen formulation (Engelund & Hansen, 1967), allowing the partitioning of sand into suspended and bedload, for which the transport is calculated separately. For information about the governing equations behind processes in Delft3D and their implementations, refer to Delft3D documentation, which is available online (Deltares, 2021b).
The deposited sediment is stored in a stratigraphic column, schematised according to the multi-layer concept (Ribberink, 1987), consisting of a transport layer, underlayers and a base layer ( Figure 1 and Text S2). The compaction formulas (Equations 4 and 5) were implemented in the underlayers of Delft3D. During compaction, the bed porosity in each underlayer decreases over time from its initial to critical values (Equation 6). The upper and lower limit of bed porosity is determined from the depositional and compacted porosity of mud and sand based on their proportion in the supply (Text S2).
Here, the bed porosity (e t ) is updated for each time interval by multiplying the porosity at the previous time interval (e t−1 ) with the ratio of bed thickness between two successive time intervals (H t−1 and H t ). In case of no compaction (H t equals H t−1 ), no change in porosity will be recorded (e t equals e t−1 ). Parameter H t is derived by subtracting bed thickness at the previous time interval (H t−1 ) with subsidence due to primary or secondary compaction (ΔH p∕s ), as shown by Equation 7. It is important to note Multi-layer administration in Delft3D, consisting of a transport layer, n numbers of underlayers, and a base layer. The transport layer has a fixed thickness of 0.2 m (th TL ). Seventy five underlayers were used in the simulation (th UL ), each having a maximum thickness of 0.3 m. The base layer has a flexible thickness (th BL ). the second-order impact of compaction, which influences the erodibility of the deposited sediment, was not included in the modelling. Therefore, the erodibility did not change regardless of how much the deposited sediment experiences compaction-induced subsidence.

DELTAS
This section discusses the model geometry and parameters used to create 4D synthetic deltas using Deft3D. Subsequently, the modelling scenarios are explained.

| Model setup
Following a similar approach as previous studies (Geleynse et al., 2010(Geleynse et al., , 2011van der Vegt et al., 2016van der Vegt et al., , 2020, the modelling setup consists of a channel debouching into a sloped basin with a 0.1° gradient ( Figure 2). The parameters used in the simulation are shown in Table 1. The stratigraphic column consists of a transport layer, 75 underlayers and a base layer containing sediment whose grain size consists of the initial bed size (sand = 100 μm) and the grain size of the supplied sediment (mud = 20 and 50 μm and sand = 100 μm). Grain sizes are homogeneously mixed in the stratigraphic column. The bed porosity of deposited sediment has initial and critical values, which are set to 0.75 and 0.08 and 0.52 and 0.25 for mud-rich and sandrich deltas, respectively (see also Text S2).
The channel conveys a constant discharge of 1,600 m 3 s −1 , which should be considered a continuous bankfull discharge. The sediment concentration in the discharge is 0.15 kg s −1 , which agrees well with modern deltas' suspended load measurements (Milliman & Farnsworth, 2011). Although the simulated deltas are fluvially dominated, it is necessary to stir up smaller-grained sediment deposited in the basin to simulate the sediment distribution similar to active coastal systems (van der Vegt et al., 2020). This was performed by including a semi-diurnal tide with a 1 m amplitude arriving perpendicular to the initial coastline ( Figure 2). No waves were implemented in the model to limit the offshore reworking force and reduce the computation time.
The total simulation time interval represents 57 days of hydrodynamic time. The simulated deltas development is extended to the geological timescale by accelerating their morphology development using the morphological scaling factor (MORFAC) of 120 and bankfull discharge (Li et al., 2018). The geological time can be computed by multiplying the hydrodynamic time with the acceleration factors, which yields 3,363 years, assuming the bankfull discharge occurs for 2 days a year (Li et al., 2018). Therefore, each time interval indicates approximately 60 years of change. Other modelling parameters are shown in Table 1.

| Scenario setup
Twenty two modelling scenarios were defined by varying the sediment supply compositions and primary and secondary compaction rates (parameters C rp and C rs in Equations 4 and 5) ( Table 2). The supplied compositions consist of 15% sand and 85% mud for mud-rich deltas, contrasting with 70% sand and 30% mud for sand-rich deltas.
The primary compaction rate (C rp ) is based on the postdepositional compaction rate of recent Holocene deposits in muddy and sandy deltas, such as the Nile and Ganges-Brahmaputra, respectively (Aly et al., 2012;Becker & Sultan, 2009;Gebremichael et al., 2018;Higgins et al., 2014;Saleh & Becker, 2018;Stanley, 1990;Steckler et al., 2022). The adjustment of the primary compaction rate using the post-depositional compaction rate is based on the following reasons: (1) No syn-depositional compaction rates measured over the delta scale (>10 3 km 2 and > 10 3 years) are available in the literature. (2) The post-depositional rates measured in the Holocene deltas are still relatively high, reaching 18 mm year −1 , indicating that compaction might still be transitioning from a syn-depositional to a F I G U R E 2 The model setup consists of a channel debouching into a sloped basin with a 0.1° gradient. The boundary condition at the channel inflow (horizontal thin black line) is water and sediment discharge, whereas the open boundaries in the basin (horizontal and vertical thick black lines) are the Neumann at the west and east combined with water level along the north boundary. The initial water depth ranges from 4 m at the delta apex (horizontal thin red line) to 21.5 m at the basin edge.
post-depositional rate. The measured rates show that sandrich deltas compact slower than mud-rich deltas. This is supported by the fact that sand only experiences primary compaction with higher compacted porosity than mud (Text S2), thereby potentially reaching compaction equilibrium before compaction measurements. Therefore, although the reported post-depositional rates of sand-rich deltas are lower than mud-rich deltas, sand compacts faster than mud. As a result, C rp for mud was set from 0 to 10 mm year −1 , while sand was set from 0 to 100 mm year −1 ( Table 2). In contrast, mud is more compacted than sand over a longer timescale, which undergoes both primary and secondary compaction. In addition, different grain geometry between mud (sheet-like) and sand (rounded) contributes to more efficient mud compaction (Revil et al., 2002).
The secondary compaction rate (C rs ) was determined by a sensitivity analysis using two criteria: (1) C rs is lower than C rp , and (2) the model is numerically stable. The C rs 3 × 10 −4 times were set lower than the C rp of mud (Table 2). It is important to note that, as deposition and erosion vary across the model domain, the actual compaction rates may also vary locally. However, these local compaction rates never exceed C rp , which is considered the maximum allowed compaction rate in the model. Figure 3 shows an example of a modelling scenario for a non-compacted mud-rich delta. The bathymetry development indicates that the delta experiences different stages of evolution, starting with a rapid basin infilling that occurs during the early simulation period (simulation time intervals 15-25), which causes seaward delta progradation. As the delta front progrades and the shoreline of active deposition reaches deeper water depth, the progradation slows down. Delta top dynamics influence the morphology evolution of the delta, shown by avulsion that occurs at simulation time intervals 30 and 50, promoting land building towards the avulsion directions. The influence of delta top dynamics decreases due to compaction, as shown by the bathymetry evolution of compacted mud-rich deltas (Text S3). In contrast, compaction has less impact on the morphology evolution of sand-rich deltas (Text S3).

OUTPUT DATA
The impact of syn-depositional compaction on the simulated deltas was analysed using their spatial and temporal characteristics. For this purpose, the Delft3D outputs containing sediment mass, elevation, water depth and porosity were post-processed for each time interval, resulting in metrics that show delta top geometry (area increase, rugosity and aspect ratio), delta top accommodation and delta-wide sediment distribution. An explanation of how to compute these metrics for mud-rich and sand-rich deltas is provided below.

| Delta top geometry
The deposited sediments were characterised into sub-environments: delta top, delta front and pro delta ( Figure 4A,B), based on their elevation. All grid cells with an elevation higher than the brink point, which indicates the transition between upstream and downstream of the delta, were classified as the delta top (van der Vegt et al., 2016(van der Vegt et al., , 2020. Distributary channel and floodplain are included in the delta top. The delta front is below the brink point and contains at least 1% sand, while the pro delta is located below the wave base. A more detailed algorithm for classifying sub-environments is explained in Text S4. The delta The aspect ratio (S dt ) was quantified by dividing the delta width (W dt in m) measured parallel to the initial shoreline with two times delta top length (L dt in m) measured perpendicular to the initial shoreline (Caldwell & Edmonds, 2014) ( Figure 4C), as shown in Equation 10. The aspect ratio of 1 indicates an ideal semi-circular or triangular delta top. In contrast, the aspect ratio less or higher than 1 (S dt <1 and S dt >1) shows an elongated or wide delta top, respectively.
The results calculated by morphology metrics were plotted for all compaction scenarios. As a gradual process, the influence of compaction on delta morphology enhances over time. To account for this effect, a trendline was calculated for each compaction scenario accounting The table shows the modelling scenarios for mud-rich (MS01-MS11) and sand-rich deltas (CS01-CS11), developed under different compaction scenarios (C rp and C rs ).

Supplied sediment concentration (kg m −3 )
C rp (mm year −1 ) for temporal variation in the morphology trend (Text S5). It was used to estimate a linear equation. Then, the area increase, rugosity and aspect ratio were determined using this equation at the end of the simulation. The calculation was repeated for all compaction scenarios. Finally, the results were plotted over compaction scenarios.

| Distribution of deposited mass
Two metrics were developed to analyse the deposition trends in simulated deltas. The first metric indicates sediment distribution alongshore, calculated by mapping the simulated deltas to a radial coordinate system consisting of three segments: western (−90° to −30°), central (−30° to 30°) and eastern (30° to 90°) ( Figure 5). Next, the cumulative deposited mass was calculated in each segment, normalised by the cumulative deposited mass in the delta at the end of the simulation. However, it is still unclear how evenly the deposition occurs alongshore. Therefore, an additional calculation was performed by taking the ratio of cumulative deposited mass in the central segment (−30° to 30°) to an average cumulative deposited mass in western and eastern segments (−90° to −30° and 30° to 90°) at the end of the simulation. The second metric indicates the deposition trends across delta sub-environments. This metric was measured by calculating the ratio between the cumulative deposited mass at each sub-environment (i.e. delta top, delta front and pro delta) and the delta's cumulative deposited mass.

| Accommodation
Compaction-induced subsidence increases relative water depth, influencing the accommodation available for sediment deposition. Therefore, two metrics were developed to estimate the accommodation at the delta top: volume reduction and water depth.
To calculate volume reduction, the subsidence due to compaction in each grid cell was first estimated by subtracting the compacted from the uncompacted bed thickness. The uncompacted thickness was measured by multiplying the porosity ratio between two successive simulation time intervals with the compacted thickness, assuming the subsidence only occurs due to porosity F I G U R E 4 (A) Plan-view bathymetry of the non-compacted mud-rich delta at the last simulation time interval. (B) The delta regions were classified into delta top, delta front and pro delta. (C) The delta top characteristics, which are shoreline length (L shore ), width (W), length (L) and total grid cells (ngrid), were used to calculate the delta top geometry metrics (area increase, rugosity and aspect ratio).

F I G U R E 5
The non-compacted mud-rich delta at the last time interval was mapped into radial coordinates centring at the delta apex, which consists of three segments relative to the delta apex from −90° to −30°, −30° to 30° and 30° to 90°. reduction (Section 2.1). This subsidence calculation concomitantly neglects the erosion effect, as porosity does not change during erosion. Next, the volume reduction was calculated by multiplying the total subsidence with the grid cell dimension, normalised by the uncompacted volume for the area within the delta top.
The second metric is water depth, which indicates the average water depth of the delta top for each simulation time interval. The calculation only includes all wet grid cells with water depths larger than 0.1 m that experience deposition of at least 15 mm. The deeper water depth reflects a more significant accommodation (Muto & Steel, 2000).

DELTA EVOLUTION
The temporal characteristics of the simulated deltas were analysed by first showing the morphological metrics (area increase, rugosity and aspect ratio). Subsequently, the sediment distribution and accommodation metrics are presented. These metrics were analysed after simulation time interval 20 when the simulated deltas reached dynamic equilibrium, indicated by relatively stable development.
To simplify the description of the modelling results, the compaction scenarios refer to the values of the maximum allowed primary compaction rate of mud (C rp ), which varies from 0 to 10 mm year −1 .

| Mud-rich deltas
The largest area increase occurs in the non-compacted scenario (0 mm year −1 ), reaching 5.6 × 10 5 m 2 /dt ( Figure 6A,B). When the compaction of 0.01 mm year −1 is added, the delta top area expands with a significantly slower rate of up to 3.4 × 10 5 m 2 /dt. The increasing compaction scenario to 10 mm year −1 leads to a gradual decrease in the area expansion rate (R 2 = 0.72). As a result, the compaction scenarios 5-10 mm year −1 have the lowest rate of area increase at the end of the simulation.
The rugosity changes between compaction scenarios, with the non-compacted scenario having the most rugose shoreline ( Figure 6C,D). The rugosity decreases linearly with increasing compaction scenarios (R 2 = 0.71). Consequently, the compaction scenarios 5-10 mm year −1 have the smoothest coastline at the end of the simulation.
The aspect ratio decreases from compaction scenario 0 to 0.075 mm year −1 , reaching the most elongated delta shape for compaction scenario 0.075 mm year −1 . The trend reverses for compaction scenarios higher than 0.075 mm year −1 , leading to a more semi-circular delta top ( Figure 6E,F). Two distinct aspect ratio trends lead to a low coefficient correlation (R 2 = 0.01). Compaction scenarios between 0.075 and 0.1 mm year −1 show the largest spread in aspect ratio values.

| Sand-rich deltas
Although less apparent than in mud-rich deltas, a slower area increase over compaction scenarios was also observed in sand-rich deltas ( Figure 7A,B). The non-compacted scenario (0 mm year −1 ) has the highest rate of delta top expansion, reaching 5.2 × 10 5 m 2 /dt, which is smaller than its mud-rich delta counterpart ( Figure 6A,B). The delta top expands at a slightly slower rate as the compaction scenario increases to 10 mm year −1 (R 2 = 0.14).
The rugosity shows a gently decreasing trend as the compaction scenario increases to 10 mm year −1 (R 2 = 0.28), indicating a slightly smoother coastline ( Figure 7C,D). No clear aspect ratio trend can be observed as the compaction scenario increases to 10 mm year −1 (R 2 = 0.02) ( Figure 7E,F). This suggests that compaction does not influence the ultimate shape of sand-rich deltas.
It is important to note that the aspect ratio is more sensitive to the delta top dynamics than other morphological metrics. This metric assumes that net deposition occurs parallel and/or perpendicular to the initial shoreline (xaxis and y-axis) (Section 4.1). Therefore, it is less responsive to a net land building that occurs diagonally to the axis, which is the case for sand-rich deltas with compaction scenarios 0.01 and 0.1 mm year −1 (Text S3). This leads to a less obvious aspect ratio trend than other morphological metrics in sand-rich deltas ( Figure 7E,F).

| Distribution of deposited mass
The observed morphological characteristics are influenced by the sedimentation patterns alongshore and downdip. First, the deposition patterns alongshore were analysed using the normalised cumulative deposited mass in three segments and the ratio of cumulative deposited mass in the central segment to an average of two side segments (Figure 8). The results show that the deposition mainly occurs in the delta centreline area for mud-rich and sand-rich deltas, which increases for higher compaction scenarios. The increasing compaction also leads to more evenly distributed sediment alongshore (R 2 = 0.18), which is particularly evident for mud-rich deltas with compacted scenarios higher than 0.1 mm year −1 ( Figure 8B). However, delta top dynamics (e.g. avulsion), which tends to redirect sediment delivery to one side of the delta, can dampen the compaction impact, leading to asymmetric deposition alongshore, such as in non-compacted to low-compacted mud-rich deltas (0-0.1 mm year −1 ) and all sand-rich deltas (Figure 8).
The distribution of sediment downdip is shown by the cumulative deposited mass at each sub-environment relative to the delta (Figure 9). The non-compacted mud-rich and sand-rich deltas show a steady increase in sediment deposition being concentrated in the delta top through time, whereas deposition in the delta front decreases through delta growth. Deposition in the pro delta remains relatively constant. This shows that the deltas gradually transition to being more delta top-dominated throughout their growth. Compaction influences this trend by causing further sediment retention in the delta top, limiting sediment delivery to the delta front and pro delta. This promotes faster transitioning to the delta top-dominated.

| Accommodation generated by compaction
The degree of morphological alteration of the simulated deltas depends on how much compaction occurs in the F I G U R E 6 The first column shows box plots, which indicate a collapse time series of the area increase, rugosity and aspect ratio for mud-rich deltas with compaction scenarios ranging from 0 to 10 mm year −1 (A, C and E). Each box plot shows a median value (second quartile, horizontal green lines) between the first and third quartiles (interquartile range, grey boxes). The vertical black lines indicate the maximum and minimum data range. The second column indicates the area increase, rugosity and aspect ratio at the end of the simulation, accounting for temporal variation of morphology trends (B, D and F). Trendlines with linear equations and coefficient correlations are also included for each morphology trend. Note that the x-axis is in the log scale. model, which differs for mud-rich and sand-rich deltas. This relationship was analysed using plots showing the volume reduction and average water depth at the delta top (Figures 10 and 11). The result shows that both mud-rich and sand-rich deltas experience increased volume reduction due to compaction through time. The trend changes over compaction scenarios, showing up to three times more volume reduction experienced by mud-rich deltas than sand-rich deltas for the same compaction scenario (Figure 10).
The delta top in natural deltas is mainly dry (except for the channels), as deposition occasionally occurs due to flooding. In contrast, the delta top in the simulated deltas receives deposition continuously due to flooding through the entire simulation time (bankfull discharge parameter in Section 3.1), which only occurs in wet grid cells activated once the threshold water depth (0.1 m) is reached.
The non-compacted mud-rich and sand-rich delta show decreasing water depth at the delta top through time. The increasing compaction scenario leads to more severe flooding at the delta top, resulting in increasing water depth through time, reaching 0.14 and 0.07 m in mud-rich and sand-rich deltas, respectively ( Figure 11). Mud-rich deltas have a greater water depth influenced by more significant volume reduction than sand-rich deltas.

F I G U R E 7
The first column shows box plots, which indicate a collapse time series of the area increase, rugosity and aspect ratio for sand-rich deltas with compaction scenarios ranging from 0 to 10 mm year −1 (A, C and E). Each box plot shows a median value (second quartile, horizontal green lines) between the first and third quartiles (interquartile range, grey boxes). The vertical black lines indicate the maximum and minimum data range. The second column indicates the area increase, rugosity and aspect ratio at the end of the simulation, accounting for temporal variation of morphology trends (B, D and F). Trendlines with linear equations and coefficient correlations are also included for each morphology trend. Note that the x-axis is in the log scale.

IMPACT OF COMPACTION ON DELTA MORPHODYNAMICS
Morphological trends have been identified for the simulated mud-rich and sand-rich deltas in response to synsedimentary compaction rates. This section combines the understanding gained from analysing the simulation results to explain the evolution of the simulated mud-rich and sand-rich deltas under compaction influence.
Compaction occurs more efficiently in the simulated mud-rich deltas as the porosity can be reduced by up to 90% compared to only 65% in the simulated sand-rich deltas for the same compaction scenario (Section 3.1 and Text S2). This is also the case in natural systems due to the plate-like shapes of the mud grains' geometry, allowing a smaller intergranular space when compacted (Revil et al., 2002). As a result, there is considerably more volume reduction in the simulated mud-rich deltas (Figure 10), leading to greater average water depth than in the simulated sand-rich deltas, a proxy for accommodation ( Figure 11). The additional accommodation in the delta top is created by compaction of the underlying delta front and pro delta deposits. More compaction results in larger accommodation available for deposition, driving the observed morphological variability (Figures 6 and 7).
Limited compaction leads to delta top-accommodationlimited deltas (e.g. the non-compacted mud-rich and sand-rich delta) (Figure 12). The supplied sediment fills the limited accommodation in the delta top. This leads F I G U R E 8 The left hand figures indicate the distribution of deposited sediment alongshore, shown by the ratio of cumulative deposited mass in each delta segment (−90° to −30°, −30° to 30° and 30° to 90°) and the delta's cumulative deposited mass at the end of the simulation for mud-rich and sand-rich deltas (A and C). To show how evenly the deposition occurs alongshore, the ratio of cumulative deposited mass in the delta centreline (−30° to 30°) is calculated to an average of two side segments (−90° to -30° and 30° to 90°) at the end of the simulation, shown in the right hand figures for mud-rich and sand-rich deltas (B and D). Trendlines with linear equations and coefficient correlations are also included. The calculation of delta segments is explained in Section 4.2. Different coloured lines and dots indicate compaction scenarios, from 0 to 10 mm year −1 .
to a rapidly aggrading delta top, indicated by decreasing water depth over time ( Figure 11). For such conditions, a relatively small portion of sediment will be stored at the delta top, whereas most of the supplied sediment will be transported downdip (Figure 9). The sediment bypassing the delta top builds the delta front and pro delta. As the down-dip deposition continues, the delta front and pro delta elevation increases, which will be classified as the delta top once their elevation reaches the brink point elevation (Text S4). This promotes rapid delta top area expansion ( Figures 6A, 7A and 12). The deposition along-shore is mostly influenced by the delta top dynamics (e.g. avulsion), resulting in unevenly distributed sediment on the western and eastern segments of the delta (Figure 8). This leads to an elongated delta top ( Figures 6F and 7F) with a rugose coastline (Figures 6D and 7D).
More compaction results in delta top-accommodationsupply balance deltas (e.g. mud-rich deltas with compaction scenarios 0.01-1 mm year −1 and sand-rich deltas with compaction scenarios 0.1-10 mm year −1 ) ( Figure 12). The deltas respond to the increasing accommodation by promoting more deposition than the delta top-accommodation-limited deltas (Figures 9 and 11) (Colombera & Mountney, 2020). This reduces sediment bypassing to the delta front and pro delta, leading to a slower delta top area expansion ( Figures 6A, 7A and 12). The influence of delta top dynamics on the deposition alongshore is dampened as the accommodation increases.

F I G U R E 9
The deposited mass at each sub-environment (i.e. delta top, delta front and pro delta) relative to the delta at each simulated time interval for mud-rich (A, C and E) and sand-rich deltas (B, D and F). Different coloured lines and dots indicate compaction scenarios, varied from 0 to 10 mm year −1 . Each simulation time interval equals approximately 60 years of delta development.
This concentrates the deposition in the delta top centreline, while the western and eastern segments of the delta receive comparable sediment, leading to a more elongated delta top ( Figures 6F and 7F) with a smoother coastline ( Figures 6D and 7D) than the delta top-accommodationlimited deltas.
Large compaction leads to delta top-accommodationdominated deltas (e.g. the compaction scenarios 5-10 mm year −1 of mud-rich deltas) (Figure 12), characterised by insufficient deposition to fill the additional accommodation, shown by the ever-increasing average water depth over time ( Figure 11A). This is remarkably similar to previous compaction studies using numerical simulation, which shows an up-dip increase in sediment trapping efficiency in response to the increasing compaction scenario (Liang et al., 2016a(Liang et al., , 2016b. As a result, limited sediment bypassing the delta top to delta front and pro delta (Figure 9), leading to the slowest rate of the delta top area increase at the end of the simulation (Figures 6A,  7A and 12). Furthermore, as the accommodation is created everywhere in the delta top, the deltas distribute the sediments more evenly across the delta top, resulting in the most semi-circular delta top ( Figures 6F and 7F) with the smoothest coastline at the end of the simulation ( Figures 6D and 7D).

| DISCUSSION
This study, for the first time, has been able to show the impact of compaction on a geological timescale evolution of delta morphodynamics. This is because the interaction between sedimentation and compaction cannot be understood by only studying the end product: the preserved sediments. In addition, experimental studies in a flume lack the sediment volume required for compaction. As compaction evolves slowly over geological timescales, its impacts on delta morphodynamics are often overlooked in observations that typically cover much shorter timescales. Therefore, it is essential to assess how compaction and delta morphodynamics interact over a longer timescale (>10 3 years), which is F I G U R E 1 0 The cumulative volume reduction in the delta top at each simulation time interval for mud-rich and sand-rich deltas (A and B). Each simulation time interval equals approximately 60 years of delta development.

F I G U R E 1 1
The average water depth in the delta top at each simulation time interval for mud-rich and sand-rich deltas (A and B). Each simulation time interval equals approximately 60 years of delta development.
F I G U R E 1 2 Down-dip topography of the delta top-limited-accommodation, delta top-accommodation-supply balance and delta topaccommodation-dominated deltas (A, B and C). The deltas respond to the increasing accommodation by promoting deposition at the delta top, limiting sediment delivery to the delta front and pro delta. This slows down the delta top expansion rate, which is mainly the case in the accommodation-dominated delta. The size of additional accommodation depends on the compaction rate of the deposited sediment, indicated by porosity values. Larger compaction resulted in thinner, denser and less porous deposited sediment. The method for classifying delta sub-environments can be seen in Text S4. an advantage in studying compaction using simulation models.
The simulations were conducted using simplified boundary conditions, such as time-invariable fluvial discharge and tides, while waves are absent in these simulations. Waves will remove the sediment from the delta top edge. The sediment is then redeposited and recompacted along the direction of wave propagation or carried offshore. This mechanism is particularly important in sediment distribution in sand-rich deltas (van der Vegt et al., 2020). Moreover, this model does not consider the second-order impact of compaction, which decreases sediment erodibility, allowing the sediment to withstand erosion (Grabowski et al., 2011;Winterwerp et al., 2012Winterwerp et al., , 2018Wu et al., 2018;Zhou et al., 2016). Nevertheless, the simulated deltas teach us that syn-sedimentary compaction is an important process that influences the morphology and preserved sediment trends. Therefore, we should be aware of the fundamental impact of this process on delta morphodynamics and preserved sediment as a whole.
The delta simulations presented in this paper are simplified versions of reality because many aspects of natural deltas are not represented, such as vegetation. Vegetated delta top tends to slow down flow velocity in the water column, allowing the sediment to settle out of suspension (Fagherazzi et al., 2012;Nardin & Edmonds, 2014). This mechanism is particularly important in mud-rich deltas exposed to strong waves, which contributes to retaining fine-grained sediment at the delta top rather than being transported offshore (e.g. the Mahakam delta) (Storms et al., 2005). Incorporating the vegetation algorithm in the simulation will enhance compaction impact by enhancing sediment retention in the delta top, which is required for land building, allowing the delta top to develop. This is an exciting topic for future studies. The simplification of the modelled deltas is also shown by using two grain-size compositions (mud rich and sand rich) supplied to the basin, which are considered very fine-grained compared to natural and previously modelled deltas. Therefore, this study should be used to help interpret the deltas with similar supplied compositions, and additional work would need to be done to understand any influence of compaction on coarse-grained deltas.
The strength of this study is to show that synsedimentary compaction that occurs at different rates does affect delta evolution, which can be placed in the context of previous grain-size studies (Caldwell & Edmonds, 2014;van der Vegt et al., 2020). The previous studies show an elongated delta top shape with a rugose coastline developed in very-fine-grained deltas. Had the compaction been considered, the delta top would have been more semicircular with a smoother coastline as the deltas respond to compaction by distributing sediment more evenly across the delta top. In the context of previous compaction modelling studies (Liang et al., 2016a(Liang et al., , 2016b, this study corroborates earlier findings that increasing compaction rate enhances sediment trapping efficiency in the delta top, leading to a smaller delta top area. Therefore, this study can help unravel the controlling processes in ancient delta deposits and predict the evolution of modern systems. The compaction algorithm resulting from this study can also be applied to predict the morphological evolution of non-deltaic coastal environments (e.g. lagoon and estuary), where long-term empirical data are unavailable. As these areas are exposed to a rising sea level, amplified by compaction (Nicholls & Cazenave, 2010), an accurate relative sea-level rise prediction is essential. This prediction relies on compaction-induced subsidence estimation of the subsurface coastal deposits, mainly consisting of mixed clastic sediments (sand and mud). This compaction algorithm allows for such prediction, combined with process-based forward modelling, which can reconstruct the coastal morphodynamic changes in response to sealevel rise. The ability to estimate the long-term evolution of the coastal regions, where millions of people reside, is a powerful tool to help regulators make better policy and management decisions.

| CONCLUSION
New compaction formulas were derived for primary and secondary compaction, validated using laboratory data, and implemented into Delft3D. The updated Delft3D code was used to generate prograding deltas with varying sediment supply compositions (mud rich and sand rich) and maximum allowed compaction rates (0-10 mm year −1 ). The simulated deltas were analysed by post-processing Delft3D outputs to derive the delta top geometry, deltawide sediment distribution and delta top accommodation metrics. It is concluded based on the analysis that: • The interaction between compaction-influenced accommodation and delta top dynamics drives the observed morphology trends. • The increase in the maximum allowed compaction rate to 10 mm year −1 imposed to simulated deltas leads to more significant additional accommodation creation, resulting in more sediment deposition in the delta top, which occurs more evenly alongshore, leading to a more semi-circular delta top with a smoother coastline. The increasing deposition in the delta top also leads to lower sediment delivery downdip, slowing the delta top expansion, resulting in a smaller area. • The influence of compaction on morphological trends (e.g. area increase, aspect ratio and rugosity) in mud-rich deltas is more significant than in sand-rich deltas. This occurs due to more considerable volume reduction, leading to greater accommodation than sand-rich deltas.

ACKNOWLEDGEMENT
We acknowledge the full-time PhD scholarship from the Indonesia Endowment Fund for Education (LPDP), which the first author receives. The authors would also like to thank Associate Editor Miquel Poyatos-Moré, and two reviewers, Elisabeth Steel & Rebecca Caldwell, for their constructive comments and advice, which helped improve the manuscript.