Quantifying the spatial distribution of sediment transport in an experimental gully system using the morphological method

Whilst time‐series of sediment transport in gullies in both laboratory experimental and field settings can be determined through instrumentation, quantifying the spatial distribution of transport rates remains challenging. The morphological method, which was proposed for estimating bed‐material transport in both one‐ and two‐dimensions in rivers, provides an alternative. Here, we developed this method for gully systems. A laboratory catchment was used to simulate gully erosion. High‐resolution topographical data were acquired by close‐range digital photogrammetry. Morphological changes were determined using high‐resolution topographic data and an associated level of detection. Based on measured morphological changes, one‐dimensional (1D) and two‐dimensional (2D) sediment transport rates were calculated via cross‐section by cross‐section routing (1D) and cell by cell routing (2D). The 1D application provided a general trend of longitudinal variation of sediment transport for the whole gully system, increased gradually from zones of headward extension to a zone downstream where erosion and deposition were in balance, and sediment transport rates less variable in space. For the 2D application, hydrological and blended hydrological‐hydraulic routing solutions were compared. We found that the level of negative transport was insensitive to whether or not a blended hydrological‐hydraulic routing was used and that results from applying the hydrological routing throughout were not significantly degraded. We also found that consideration should be given to spatial and temporal resolution of the topographic data. The 2D application provided spatial patterns of sediment transport that vary with gully evolution. The main gully remained a high transport corridor but branch transport became more important through time. The framework we report provides an additional tool for both experimental and field quantification of the spatial patterns of sediment transport in gullies; and quantification of how these patterns change under different forcing factors.

provided a general trend of longitudinal variation of sediment transport for the whole gully system, increased gradually from zones of headward extension to a zone downstream where erosion and deposition were in balance, and sediment transport rates less variable in space. For the 2D application, hydrological and blended hydrologicalhydraulic routing solutions were compared. We found that the level of negative transport was insensitive to whether or not a blended hydrological-hydraulic routing was used and that results from applying the hydrological routing throughout were not significantly degraded. We also found that consideration should be given to spatial and temporal resolution of the topographic data. The 2D application provided spatial patterns of sediment transport that vary with gully evolution. The main gully remained a high transport corridor but branch transport became more important through time. The framework we report provides an additional tool for both experimental and field quantification of the spatial patterns of sediment transport in gullies; and quantification of how these patterns change under different forcing factors.

K E Y W O R D S
1D and 2D transport rates, gully erosion and deposition, morphological change detection, sediment routing, spatial pattern of sediment transport, the morphological method

| INTRODUCTION
Gully erosion is one of the main sources of sediment supply to many drainage basins (Valentin et al., 2005). For instance, in the loess plateau of China it is thought to contribute to between 60% and 70% of all eroded sediment (Li et al., 2003) and this figure has been reported to be > 85% in a subtropical catchment in south China (Luk, 1997). Estimating the controls on sediment transport in gully systems is crucial for understanding the mechanisms by which gullies respond to driving factors like climate, vegetation and land-use; but also for developing predictive formulae that may help with land-use management; and for validating numerical models that can predict gully erosion in its own right or as part of wider soil erosion modelling.
A number of soil erosion models now include gully erosion at both the hillslope and catchment scales, including the Annualized Agricultural Non-Point Source (AnnAGNPS) model (Gordon et al., 2007); the Revised Universal Soil Loss Equation (RUSLE) (Fu et al., 2005); the Limburg Soil Erosion Model (LISEM) (De Roo et al., 1996); and the European Soil Erosion Model (EUROSEM) (Morgan et al., 1998). These models allow prediction of sediment budgets, sediment transport capacities and sediment transport rates; although they can also have considerable uncertainty and may be highly demanding in terms of data for input and for calibration (e.g., Hessel et al., 2003). Evidence also suggests that predictions can vary significantly between models (Wang et al., 2013). The development and application of such models is likely to be enhanced via better measurement of sediment transport rates through time and in space in gully systems. Whilst time-series of sediment transport in both laboratory experimental (e.g., Cui, 2002) and field settings (e.g., Bhatti et al., 2011) may be determined through instrumentation, quantifying the spatial distribution of transport rates remains a more challenging goal.
A promising alternative involves inferring sediment transport in space from measured morphological change through time. In studies of rivers, this has been named the 'morphological method' (Ashmore & Church, 1998;Lane et al., 1995). The method quantifies the sediment transport necessary to conserve mass following from measured morphological change (Vericat et al., 2017). If data on morphological change are distributed in space, the sediment transport rate estimates are also spatially-distributed. Initial applications of the morphological method were limited by the spatial and the temporal resolution with which topographic data could be acquired (Ashmore & Church, 1998;Fuller et al., 2003;Lane et al., 1994aLane et al., , 1995. With ongoing development of data acquisition technologies, this limit is now being addressed. Laser scanning has been used for very highresolution gully mapping (Baruch & Filin, 2011;Hegeman et al., 2014;James & Quinton, 2014;Schneider et al., 2012;Shellberg et al., 2013).
Although the morphological method has been applied to quantify sediment transport rates in rivers Ashmore & Church, 1998;Lane et al., 1995;Vericat et al., 2017) and to estimate sediment yield on hillslopes in formerly glaciated terrain (Heckmann & Vericat, 2018), it has yet to be applied to gully erosion and deposition.
Gully systems may be particularly suitable for application of this method because, unlike in rivers, gullies are mostly dry between sediment transport events, overcoming the complications in rivers of obtaining data from inundated areas. It is also possible that applications in gullies avoid a second problem. In rivers, there are normally two boundaries across which sediment can flux, one upstream and one downstream. If one of these is known, we can route it in upstream or downstream direction using the morphological method and then the absolute sediment transport rates are obtained; otherwise, the morphological method only gives the change in sediment transport due to erosion and deposition and not the absolute sediment transport rates . For gullies, their spatial extent is much smaller than rivers so that field surveys can easily cover all gully heads. Provided the spatial extent of the survey is sufficient to include all gully heads, the flux at gully boundaries (upstream boundary condition) is zero. The boundary condition problem is eliminated and absolute transport rates can be determined throughout the measured domain, as well as the total sediment yield estimated.
Given the reduced boundary condition complication and the ease with which morphological data can now be measured in gully systems, the morphological method clearly has potential. There is one challenge compared with rivers. Gullies contain a much greater range of relief, including an important transition from the hillslopes, where sediment transport is likely to be hydrologically-controlled according to a flow routing defined by topographic slope (e.g., hillslopes in Heckmann & Vericat, 2018), to the gully bottoms where sediment transport is likely to be hydraulically-controlled by flow routing defined by interactions between the water surface slope and the gully bed slope during a sediment transporting event (e.g., rivers in Antoniazza et al., 2019). Application of the morphological method to gullies needs to be able to represent these two differing controls on sediment routing.
The aim of this article is to develop the morphological method for determining the spatial patterns of sediment transport in gully systems (required to conserve mass) using approaches now wellestablished for rivers. This is done in both one-and two-dimensions.
To minimize the uncertainties associated with data collection, the focus is upon a laboratory-simulated loess catchment, rather than field data. We include tests on the key parameters that influence results from the morphological method.
The methodology is presented in Figure 1. First, we used a laboratory loess catchment to simulate gully erosion, which allowed ready quantification of morphological change as well as precise determination of hydrological parameters. Second, close-range digital photogrammetry was used to construct digital elevation models (DEMs) and DEMs of difference (DoD), and to determine their uncertainty so as to specify the magnitude of detectable changes as a level of detection. With this level of detection, significant morphological changes were identified.
Third, cross-section by cross-section routing was developed for application of the morphological method in one dimension. Fourth, the same was undertaken in two dimensions. The two-dimensional (2D) method requires a routing treatment and this was done in two different ways: (1) hydrologically, using topographic forcing of hillslope routing throughout the gully system; and (2) blended hydrological-hydraulic, in which routing in the gully bottoms was replaced with a 2D hydraulic treatment. These treatments were compared. Finally, we undertook simulations to identify the key influences upon application of the method.

| Experimental design and data acquisition
The experimental work was undertaken at the rainfall erosion test facility at the State Key Laboratory of Soil Erosion and Dryland district, Xianyang, China. During filling, the loess was compacted layer by layer every 5 cm until the bulk density was within the range 1.36-1.4 g/cm 3 , and the final average bulk density was 1.39 g/cm 3 .
The initial topography (Figure 2a) was designed to represent an initial valley system with no gullies present.
A total of 25 rainfall events were applied to the experiment over a 2.5 month period. Due to equipment limitations, we designed three types of rainfall intensity: 0.5, 1, and 2 mm/min, which represents light, moderate, and heavy rain on the Loess Plateau of China (Changxing et al., 1995), and account for 44%, 36% and 20% of the designed rainfall events, respectively. The details of the simulated rainfall events are shown in Table 1. During the rainfall events, we used a sump at the outlet to collect all sediment (both suspended sediment and bedload) and runoff ( Figure 1). The trapped sediment was sampled, dried, and measured and the average sediment transport (export) rate for each rainfall event was determined (Table 1). We also measured the volume of runoff for each rainfall event and calculated the average discharge (Table 1). Topographic surveys were undertaken with a lower frequency (Table 1). In order to link change between topographic surveys to rainfall and runoff data, we calculated the average rainfall intensity and accumulated rainfall duration between surveys (actually only 23 rainfall events affected the topographic surveys) ( Table 1).
The DEMs were produced using close-range digital photogrammetry. A JPL Carl Zeiss SMK 120 Stereo Photogrammetry Camera (Carl Zeiss AG) was used to acquire imagery. The JX-4 Digital Photogrammetry Workstation (Beijing Geo-Vision Tech. Co., Ltd) was used for the photogrammetry. The latter is a conventional photogrammetric platform with a built-in camera calibration system. Unlike recent photogrammetry software (such as Pix4D and AgiSoft Metashape), the JX-4 Digital Photogrammetry Workstation produced a gridded stereo model and then generated DEMs. Eighteen control points and 20 check points were evenly distributed in the study area, and a local coordinate system was employed to survey the control and check points using an electronic theodolite. The horizontal and vertical root mean square positioning errors (RMSEs) of these points were better than ±0.3 mm. The control points were used in the bundle adjustment.
DEMs were generated to a resolution of 10 mm. Check points were used to assess DEM quality in terms of mean check point error and standard deviation of error (Table 1). The vertical mean errors were between −0.09 mm and 0.28 mm and the vertical standard deviations of error were within the ranges ±1.64 mm to ±1.98 mm (Table 1).
Visual inspection of the DoD suggested that the combination of the design of image acquisition (use of a high specification camera) and the use of control points had eliminated systematic error such as doming, reflected in standard deviations of error that are commensurate with the resolution of the acquired imagery (James et al., 2019).

| Morphological change
For each time period between DEMs, we subtracted the second DEM from the first DEM such that positive change represented erosion and therefore a contribution to sediment transport; and negative change deposition and a decrease in sediment transport. Prior to this, we removed the mean (i.e., systematic) error (Table 1) (Table 1) under the assumption that errors were point by point independent, random and Gaussian (Brasington et al., 2003;Lane et al., 2003): where δ DoD labels the propagated error of DOD; δ a and δ b are the standard error of DEM a and DEM b, respectively. We applied a Student's t test at 95% confidence to determine a threshold of the level of detection (Lane et al., 2003;Wheaton et al., 2009). Only the cells in the DoD with values outside the thresholds (1.96×δ DoD ) were regarded as real topographic changes and used in the calculations. There is some debate as to whether such thresholding should be applied (e.g., Anderson, 2019). When undertaking sediment routing calculations with erosion and deposition data, noise in the latter automatically leads to variance in the sediment transport rates that is spatially non-linear (e.g., the effects will be greater at gully heads where there is close to zero transport than further down gully where there may be higher transport rates). For this reason, we think it is important to address this problem using a threshold. It may be a problematic assumption for situations where there is spatially extensive, coherent but very small magnitude change, a problem that can occur in braided rivers . In gullies, especially with the high quality of data used here, change tends to be spatially concentrated and high magnitude rather than extensive and low magnitude. Note a crucial assumption in our analysis is that the gullies erode by surface erosion and not subsurface erosion; the latter can occur in gully systems (Nouwakpo & Huang, 2012), and may lead both to surface lowering that is not belonging to surface erosion and deposition and sediment transport that does not follow the surface flow paths. This is a weakness of the method.
We also added in one other calculation. If we sum positive and negative DODs across space to get a total volume eroded and a total volume deposited, respectively, we have an estimate of the net volume of sediment that should have been lost from the basin. After correction for bulk density, we can then determine a sediment export mass which we convert into kilograms per minute (kg/min) by dividing by rainfall duration; we call this the sediment export rate. We also measured the sediment export rate at the basin outlet directly (Table 1), and so we could compare the DEM-estimated and the measured sediment export rates, which is ultimately a validation of the global mass conservation of sediment under the assumptions we make (porosity, surface lowering only by surface erosion, DEM reliability, etc.).

| One-dimensional application of the morphological method
The morphological method is an inverse-solution to estimating sediment transport rate based upon the Exner equation (Exner, 1925) for sediment transport: T A B L E 1 Summary of experiments: rainfall, discharge and digital elevation model (DEM) quality measured against check points (Cui, 2002) Rainfall event (photography event) where Q b is the sediment transport in the x and y downstream and cross-stream directions respectively; p is the sediment porosity; Z is elevation in position (x, y); t is time and C b is the concentration per unit bed area of sediment in motion. By assuming negligible lateral sediment transport, and assuming that C b is constant in space and time, Equation 2 becomes in one dimension (Vericat et al., 2017): where Q b is the cross-sectionally averaged sediment transport rate; and Z is the cross-sectionally averaged morphological change.
For the one-dimensional (1D) application of the method, net volume changes per cross-section are measured by DEM differencing.
The net change (ΔV j ) in any one time period between adjacent sections j is then determined . A discretization of Equation 3 is then used to calculate a sediment transport rate in kilograms per second: where S j is the transport for a given cross-section j; S j-1 is the transport rate coming from the upstream cross-section; ρ is the material density; p is the porosity; ρ(1 − p) is equal to the bulk density in practice; ΔV j is the net volume change by cross-section measured by DEM differencing and t is the transport duration.
The main difficulty in a gully system is specification of cross-sections. A priori, this needs the definition of a gully network and then cross-sections at right angles to the gullies. As gullies evolve, the orientation of the cross-sections should change, making it very difficult to compare consecutive cross-sections in time. For this reason, we make a much simpler approximation which involves fixed crosssections oriented parallel to the x coordinate direction, and hence the orientation of the primary gully 2.5 | Two-dimensional application of the morphological method: Principle Two-dimensional application of the morphological method was first proposed and tested by Lane et al. (1995) for a small section of a braided river. This, and a much larger scale application to a braided river , used the output of depth-averaged 2D hydraulic models to route sediment spatially. For the case of gullies, a similar routing treatment is needed. The path followed by sediment will be a function of the resolution of two components of shear stress; that due to flow; and that due to gravity or topography . Following Nelson and Smith (1989), the resultant stresses are defined by:  where the first term of the right-hand side of Equations 5a and 5b describes the flow-related shear stress; j u j = ffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 x + u 2 y q is the magnitude of the velocity vectors in planform (u x ,u y ); the ρ w is the water density (1,000 kg/m 3 ); g is the gravitational acceleration (9.81 m/s 2 ); d is the water depth; and n is the value of Manning's n; the second term on the right-hand side represents the topography-related shear stress (Nelson & Smith, 1989); τ c is the critical shear stress defined from a Shields condition; θ is the arctan of s, the slope derived from the DEM; ϕ is the bulk angle of repose of the sediment. The two terms on the right-hand side are not independent as local topographic slope will influence the water surface slope and the momentum equations that lead to the velocity terms in Equations (5a) and (5b). When the topographic slope is low, the magnitudes of the flow-related and topographic terms will be similar; but as the slope steepens, so the topographic terms will become dominant. Even in flatter gully bottoms, gully bottom slopes may be substantially steeper than in a river, meaning that the topographic term dominates. To clarify definitions from here onwards, we refer to the first term on the right-hand side as flow forcing and the second term as topographic forcing. When just topographic forcing is included, and the flow itself does not influence routing, then we have a problem equivalent to hydrological routing.
When both flow forcing and topographic forcing are included, then we call it hydraulic routing.

| Two-dimensional application of the morphological method: Hydrological routing
If topographic forcing is dominant then we assume that sediment follows topographically driven flow paths with only s x and s y controlling routing; that is it can be treated as a classical hydrological routing problem. The right-hand term in Equations 5a and 5b is effectively a generalization of a multiple flow path routing, with two flow paths possible. Here we allow for multiple possible flow paths using the Holmgren (1994) function. Each cell receives a certain volume of sediment from surrounding cells that are higher in altitude, Q s . It then supplies sediment to cells that are lower in altitude in proportion to their slope. Hence, according to the Equation 2 and the multiple flow direction model (Holmgren, 1994), the sediment transport rate in any direction k, Q k b , is defined as: where k is the index of each of the eight surrounding cells that could deliver sediment to cell ij if they have positive slopes; s k is the slope of the surrounding cells k; ΔV ij is the net volume change measured by DEM differencing; t is the duration time; ρ(1 − p) is equal to the bulk density in practice; and α is the parameter that controls the degree of flow diffusion (or concentration) which tends to a classic D8 algorithm as α tends to infinity, and tends to multiple flow direction routing with an α of 1. Flow routing was applied in using the TopoToolBox system (Schwanghart & Kuhn, 2010).
Use of Equation 6 includes the parameter α which is not known a priori. Given the difficulty of measuring spatially distributed sediment transport rates, there are no independent calibration data. However, it is possible to use a no negative transport condition (e.g., Antoniazza et al., 2019) which states that there cannot be more deposition in any one cell in the model than there is sediment supplied to that cell. If the routing is not correct, then negative transport situations may arise. In the gully system case, and unlike in rivers, we can assume that sediment transport above the gully heads is zero, we do not need to specify upstream sediment supply conditions and so negative transport can be attributed either to errors in the morphological change estimation (i.e., the photogrammetrically-acquired DEMs) or the routing process itself. Thus, the extent of negative transport is an index of the quality of the sediment routing and we use it here to assess dependence on α. We ran the hydrological routing 200 times with α randomly sampled between 1 and 10 for the second period. and a Manning's n of 0.03 (after Arcement & Schneider, 1989), for the initial default simulations. It should be emphasized that it has been shown that in 2D hydraulic models, perturbations of Manning's n have only a very small impact on flow routing (Lane et al., 1994b;Yu & Lane, 2006).
To couple the hydraulic routing to the hydrological routing, we treated each branch of the main gully as an inlet boundary condition ( Figure 3). This is supported by the observation of a marked reduction in slope as branches enter the main gully. At each inlet, we needed to supply a flow rate. This was done by taking the discharge measured at the experimental outlet in the main gully and distributing it to each branch in proportion to the latter's upstream contributing area: where Q i is discharge of inlet i; A i is the upstream contributing area of inlet i; n is the total number of inlets; Q out is the recorded outlet discharge.
At the outlet, we specified a normal depth condition with model extension downstream such that this condition had no upstream influence upon the simulated zones of interest.
To avoid semi-wetted cells and numerical instabilities, a minimum water depth was required for a cell to be defined as inundated (Ingham & Ma, 2005). Here, we used a minimum water depth of 10 mm. Hydraulic simulations were performed with a fixed bed of the initial or final DEM (first or second DEM of the survey period) in every period. The DEM choice followed the hydrological routing; initial for situations where the hydrological routing used the initial DEM, and final where it used the final DEM. In all cases, BASE-MENT was run until a steady state had been reached with a mass balance between input and output instantaneous discharge error lower than 0.1%.
When a flow-related shear stress is combined with a topographyrelated shear stress as in Equations 5a and 5b, parameters are needed.
For the critical shear stress, we used classical formula (Parker et al., 2003) to calculate the particle Reynolds number (Re p ): where Re p denotes the explicit particle Reynolds number; R = 1.65 is the submerged specific gravity of the sediment; ν = 10 −6 m 2 /s denotes the kinematic viscosity of water and is evaluated at 20 C; and D is the median grain size. Then, according to a Shields curve for cohesive silt (Miedema, 2013), the Shields parameter was determined.
Finally, the critical shear stress was calculated from (Ni et al., 2020;Parker et al., 2003): where τ Ã c is the critical Shields number. Equations 5a and 5b were then applied following Antoniazza et al. (2019): where t is the duration of competent flow; k is the index of each of the eight surrounding cells that could deliver sediment to cell ij; and ΔV ij is the net volume change measured by DEM differencing.
Whilst the net volume change ΔV ij and critical shear stress τ c can be determined by the DoD and median grain size, respectively, the angle of repose and Manning's n are still uncertain in Equation 5a and 5b due to a lack of information on sediment composition. There is also no a priori reason why the Manning's n values needed in Equations (5a) and (5b)  Unlike when the hydrological routing method is applied across the entire gully system, the blended method requires the output from the hydrological treatment to be coupled to the hydraulic treatment.
Although the calculation of the hydraulic routing takes discharges from only the gully branches, assuming direct hillslope runoff to the F I G U R E 3 The inlets and outlet shown superposed on the second relief shaded DEM [Colour figure can be viewed at wileyonlinelibrary.com] main gully to be negligible, the sediment transport is fully coupled, with sediment routed hydrologically on the hillslopes and in the gully branches delivered to any point in the main gully zone where the hydraulic routing was applied.

| Effects of spatial resolution and temporal frequency on estimated transport rates
Finally, as research suggests that the spatial and temporal resolution of data can impact the results obtained using the morphological method (Ashmore & Church, 1998;Lane et al., 1994aLane et al., , 1995, we tested their effects. The second and third DEMs were resampled to 50, 100, 200 and 500 mm for the second survey period for investigating the effect of DEM resolution. Alternating erosion and deposition during a survey period could compensate morphological estimates.
Recorded morphological changes may vary with survey frequency within a certain period, which affects sediment transport estimations.
Thus, we also looked at how topographic survey frequency (i.e., temporal resolution) impacted results. The transport rates for the whole period studied (Table 1) were calculated from all nine DEMs, five DEMs (i.e., the first, third, fifth, seventh and ninth DEMs were used), three DEMs (i.e., the first, fifth and ninth DEMs were used), and two DEMs (i.e., the first and ninth DEMs were used), respectively.

| Erosion, deposition and net changes at the basin scale
The eight DoDs calculated from the nine available DEMs are shown in Figure 4, after application of the level of detection. The DoDs show erosion (red) at the gully heads and deposition (blue) in the gully bottoms. The gullying process involves erosion both laterally (gully widening) and headward (gully extension). As the gullies widen and lengthen, deposition can occur downstream, notably in periods 3, 5 and 8; but this is not a systematic evolution. Thus, the sediment yield (erosion minus deposition in any one period) does not systematically decline through time (Table 2). Erosion and deposition (  Figure 4). This will lead to erosion and deposition within a cross-section cancelling itself out, when applying the 1D morphological method (Lane et al., 1995).
For the whole catchment, the sediment export rates calculated from the DoDs ranged from 3.72 to 20.66 kg/min (Table 2). There was a good general agreement with the observed sediment export rate (Table 1) (r = 0.97, P < 0.01, Figure 5), suggesting that inference of sediment transport for morphological changes is reliable at the basin scale. This was particularly the case for periods 1, 2, 3, 5, 6 and 8 ( Figure 5) and suggests that surface erosion is likely to have been dominant for these periods. For periods 4 and 7, more sediment was observed as exported than was inferred from the DoDs (Figure 5), and this could be due to subsurface erosion (i.e., erosion that has no surface expression) during these two periods. It could also be due to erosion that has a depth that is smaller than the limits of detection of the photogrammetry. It should be emphasized that the calculation does not confirm whether or not the sediment routing implicit in 1D and 2D approaches is correct, and that the associated spatial distribution of transport is correct.
3.2 | One-dimensional sediment transport rate Figure 6 shows the results from applying the 1D morphological method. The positive slope of the transport rate versus distance means that the expanding gully network results in net sediment erosion from upstream to downstream even through there may be both deposition and erosion in any one cross-section ( Figure 4). Except for the first time period, the sediment transport rates tend to become stable at a certain distance downstream (i.e., cross-section 700-900 cm in Figure 6), especially in the latter periods (e.g., periods 7 and 8). As 3.3 | Two-dimensional sediment transport rate: Hydrological routing but it is also in the flat zone where there is substantial negative transport ( Figure 6) and this suggests that in the gully bottom hydrological routing moves sediment incorrectly because it is too sensitive to small topographic variations. These observations aside, the complete hydraulic routing (with both the flow and the topographic forcing) does not reduce the proportion of negative transport cells significantly in the main gully (  Figure 10). Table 3 shows the slight improvement that comes with this blended routing for all periods. Compared to the 1D results, the 2D results show substantial variation of sediment transport rates for a given value of the y-direction (Figure 10), as well as zones of convergence (e.g., 700-800 cm in the downstream direction in the fifth F I G U R E 1 0 Two-dimensional sediment transport rate estimates using blended hydrological-hydraulic routing (note that the negative transport rate is set to zero for visualization) [Colour figure can be viewed at wileyonlinelibrary.com] period in Figure 10) and divergence (e.g., 500-600 cm in the downstream direction in the fifth period in Figure 10) of transport. High transport rates (values bigger than 1 kg/min) are concentrated in the main gully due to the effects of accumulation in the dendritic drainage network. Low transport rates (values between 0 and 0.5 kg/min) are commonly found in the gully branches, but also at the margins of the main gully.
Although the transport rates in the main gully were always the highest, we found that branch transport rates increased with the gully system development ( Figure 10). Then, we calculated the proportion of the branch transport in every period, that is the ratio of the sum of transport rates at branch outlets to the sediment export rate of the catchment. Figure 11 shows that the proportion of branch transport and the number of branches. As the gully developed, the proportion of branch transport increased at first and then became stable. The number of branches showed the same pattern; and a positive correlation was found between the two (r = 0.86, P < 0.01). This finding means that although the main gully maintains a high transport, it becomes less important in terms of proportion as the gully system develops.

| Comparison of 1D and 2D approaches
To compare with the 1D results, the 2D results produced by blended hydrological-hydraulic routing were transformed into their 1D form by integrating across a given cross-section. The 2D results are more variable as a function of distance downstream ( Figure 12) and also plot higher in terms of total transport rate. This finding is not surprising. The 1D method accumulated the DoD within cross-sections at first and then routed them downstream, which inevitably compensates some zones of erosion (positive in a DoD) and some zones of deposition (negative in a DoD) within cross-sections ( Figure 4). The 2D method routes each cell in the DoD and so reduces the possibility for local compensation. When we integrated the 2D results back onto the 1D treatment (Figure 12), we obtained higher transport rates than the 1D method in many cross-sections ( Figure 12). The intensity of lateral transport within a cross-section in the presence of net erosion drives this effect; if any cell of eroded sediment is transported laterally (x coordinate direction) then the sediment contributes to more than one cell in the x direction; such that calculation of the integration across x direction causes higher increases in transport rate than with a 1D calculation. The reverse happens with net deposition in a section (the mean reduction in transport rate is higher). Note that the two solutions are still mass conservative. That said, the broad patterns of the 1D method and the integration of 2D method are the same: the transport rate increase with distance downstream.

| The effect of DEM resolution
Here, we used the hydrological routing to analyse the effect of DEM resolution because the hydrological routing is more straightforward and is easier to apply than the hydraulic routing. We use the second period of data for the test. With the DEM resolution becoming coarser, the spatial extent of sediment transport becomes wider ( Figure 13a-e). Some hillslope areas without erosion were mistakenly regarded as eroded areas, leading to over-estimation of transport extent. To clearly see how the value of transport rates changes, we transformed the 2D result into their 1D form by integrating them across the x (lateral) direction. As Figure 13(f) shows, the coarser resolution DEM tends to result in lower transport rates. This finding suggests that the coarser resolution DEM causes higher transport in low transport areas and leads to lower transport in high transport areas. Besides this local effect, the global sediment transport rates (i.e., sediment export rate) decline with increasing coarseness of the topographic data. Figure 13(f) also shows how with coarser DEMs, the reduction in the quality of routing estimation can lead to rapid decreases in transport rate (e.g., transport rates at 680-700 cm downstream), largely because sediment becomes stuck in parts of the gully system due to local averaging in the presence of steep topographic gradients. 3.8 | The effect of survey frequency Figure 14 shows sediment transport estimations with various survey frequencies for the whole period studied (Table 1). The sediment transport is much concentrated in main gully when only two DEMs were used (Figure 14a). As the survey frequency increases, there are more details of transport paths; and the paths in main gully become more diffusive (Figure 14a-d). We also integrated the 2D results into their 1D form (Figure 14e). The 1D form reveals a global effect of survey frequency: transport rates generally tend to increase with survey frequency, especially in downstream areas. This is likely to reflect the effects of compensating erosion and deposition, confirmed by a greater impact further downstream.

| DISCUSSION
The 2D morphological method was proposed more than 20 years ago for rivers (Lane et al., 1995). This article develops the morphological method for sediment transport estimates in a gully system. The method seems to work as the negative transport condition is generally met and the proportions of cells with negative transport are very low.
Progressive gully development creates a dendritic drainage network ( Figure 4) and this serves to focus sediment into the main gully, where deposition occurs. Although the rainfall intensity varies between periods (Table 1), the results show that there is enough flow accumulation to maintain significant sediment transport to the catchment outlet (Figures 6 and 10). This finding is different from braided rivers where continuous flow convergence and flow divergence leads to a strong autogenically-driven control on the braiding process .
The results suggest that the 1D approach should be applied with care to gully systems as it is hard to specify cross-sections for an evolving dendritic network. In this article, we used a simple approximation with cross-sections oriented orthogonally to the y coordinate direction, itself parallel to the main gully direction. This approximation method allowed visualization of a general trend of variation in sediment transport from upstream to downstream for the whole gully system and through time ( Figure 6). In all cases, transport rate increased with distance downstream through the gully system, eventually reaching a plateau. This plateau occurs approximately where the last branch joins the main gully and implies balanced erosion and deposition within the downstream main gully. The fact that from this point sediment delivery no longer increases but that sediment transport continues, and that in this zone, morphological changes are substantially lower than towards the gully heads ( Figure 4) suggests that the main gully has evolved to be able to transport all sediment supplied to it. This observation is also reflected in sediment export rates ( Table 2) that are of a similar magnitude to the plateau ( Figure 6).
Thus, as the gullies develop the lower parts evolve to being sediment transport conveyors whereas the gully heads continue to grow. This is reflected in the 1D sediment transport rates which begin to rise pro- Compared with applications of the 2D morphological method to rivers Bakker et al., 2019;Lane, 1997;Lane et al., 1995), the application to gully systems has fewer demands in terms of hydrological information; and it has a natural boundary condition, the gully heads, across which sediment transport rate must be zero. The datasets that result may provide a better means of calibrating and validating models of sediment transport in gullies (Fu et al., 2005;Gordon et al., 2007). It is possible that we can use the derived sediment transport fields to understand the spatial structure of sediment transport and to relate this to distinct geomorphic process signatures (e.g., Llena et al., 2020) and the space-time evolution in within gully sediment connection (Heckmann & Vericat, 2018) and ultimately to test and to refine sediment connectivity indices (e.g., Cavalli et al., 2013). The data can also be used to produce spatially-explicit sediment delivery ratios (Heckmann & Vericat, 2018) and to describe how these evolve through time.
This positive conclusion aside, there are some recommendations and limits that follow. First, the choice of DEM may affect sediment routing. Previous application of the morphological method used initial DEMs to drive sediment routing Bakker et al., 2019;Heckmann & Vericat, 2018). In our case, we compared the initial DEM, the average DEM, and the final DEM and found that the initial DEM gave marginally better results in terms of the negative transport condition (Figure 7b). The problem arises because the method uses a single initial topography for routing. However, the erosion-deposition process will modify topographic conditions continually (Mosselman, 2005;Nicholas et al., 1995) so that flow paths change through time during an event. The extent of the problem will depend upon both the time between surveys and the rates of morphological change (Heckmann & Vericat, 2018). There is likely to be a point at which the time between surveys is so long or rates of morphological change so great that the actual sediment transport pathways diverge significantly from those defined by the initial DEM. The best survey frequency needs to be judged with reference to rates of morphological change. Hence, it is difficult to produce a single recommended survey frequency other than to state that, for the most extreme events, the morphological method is likely to break down completely.
Second, this study tested a range of different approaches to the routing problem: hydrological; and hydrological blended with hydraulic. This reflected the fact that if there is divergence between the slope of the water surface and the bed surface then routing will be hydraulically controlled, and both flow forcing and topographic forcing will drive routing. This is the assumption made in 2D applications to rivers . On hillslopes, as the topographic forcing is likely to be significantly greater than the flow forcing, hydrological routing may be sufficient. In quantitative terms, using the negative transport condition, hydrological routing on the hillslopes and hydraulic routing in the bottom of the main gully resulted in a marginal improvement (Table 3). However, it was not great. Visually, there was a more marked difference in the results, with the hydraulic routing leading to greater diffusion of sediment transport paths, notably at the downstream end of the gully ( Figure 10). We have no further evidence to suggest which of these routing treatments is most likely to be correct except that as a gully becomes wider it is more likely to need a hydraulic treatment. For most of the gullies in this system, the hydrological treatment seems to be able to respect the negative transport condition as well as the blended hydrological-hydraulic treatment and it is therefore acceptable. It is also substantially easier to apply.
Third, coarser-resolution DEMs result in lower sediment transport rates in general, but over-estimate the extent of sediment transport as compared with the finer-case (Figure 13a-e). There is also some evidence that DEM coarsening may lead to local topographic artefacts that can substantially impact transport rates (Figure 13f). Research (Dai et al., 2019;Gómez-Gutiérrez et al., 2015;Saksena & Merwade, 2015) has addressed the effects of DEM resolution on gully characterization and shown that it depends on the scale and morphology of gully itself. It is quite probable that the resolution of topographic data will be defined by the scale of the gully; but with developments in SfM-MVS photogrammetry, such data are increasingly likely to be available. Thus, spatial resolution tests should be a routine element of understanding how the morphological method is performing in any given application, but new technologies are likely to be able to provide the necessary high-resolution data.
Fourth, erosion and deposition can alternate through time (Balaguer-Puig et al., 2017;Lindsay & Ashmore, 2002;Milan et al., 2007). High frequent surveys can capture more morphological changes, reducing the probability of compensation by alternating erosion-deposition through time, resulting in higher transport rates and providing more details of transport paths ( Figure 14). However, in practice, such high frequent surveys are unfeasible. Lane et al. (1994a) showed significant bias in the estimation of volumes of morphological change if the survey frequency does not match the rate of process change in time in a river. In this study, the time periods contained more than one rainfall event. It would probably be wise to increase survey frequency such that only one rainfall event occurs between surveys. However, given the ease of data acquisition and DEM generation using current technologies, as well as the ease of its application using hydrological routing, the morphological method could become a routine manner for experimental or field-based evaluation of controls on gully erosion processes. Even with a good survey frequency, sufficient to avoid compensating erosion and deposition, the actual estimated transport rates assume that the associated rates of Fifth, the 2D application is dependent upon the parameters that may be uncertain a priori (i.e., angle of repose and Manning's n, or exponent of α). On the one hand, although the parameters can be determined a posteriori via the kind of Monte Carlo parameter sensitivity analysis used here, with the negative transport condition as target, such methods do depend on the reliability of this condition.
Indeed, this parameterization approach runs the risk of identifying empirically-adequate but physically implausible parameters. Given that the negative transport condition did not change that much in response to parameter changes in this study, placing too much emphasis on statistically-identified more probable parameter values should be avoided. On the other hand, whilst the negative transport condition is a means of parameterization, it is a condition that needs careful application, not least because negative transport may come from many different sources. Given that we assume that error in the DEM follows a Gaussian distribution, we would expect a certain number of errors to remain in a DoD even after thresholding (the 5% of observations that fall outside of the detection limits). Such data points may lead to the localized deposition of sediment that in turn leads to insufficient sediment being delivered downstream. Despite being caused by a local error, it could explain the intermediate sized zones of negative deposition that we see for example in Figure 7(a). This is unlikely to explain the large zone on the bottom right of Figure 7 (a) which seems to be associated with an entire gully branch. This latter suggests a more systematic error in the calculation. The DEM error aside, if subsurface erosion occurs, the routed sediment may be delivered to a zone different to that where it is delivered by the subsurface and then produce negative transport.
Finally, compared to the laboratory case study, field application of the morphological method should address some additional issues that may become problematic. On the one hand, the experimental set up used here allowed us to determine the total transport rate because for the most part ( Figure 5) our global mass conservation was respected, that is; our observed exported sediment included both bedload and suspended load; and except for two periods, this matched our DoD estimates. However, in a field application, it is possible that there is sediment loss both below surface or leading to surface changes that are too small to detect feasibly with survey methods. In such cases, the method will under-estimate sediment export rates because it under-estimates actual erosion. On the other hand, hydrological data (e.g., average discharge) are necessary for hydraulic simulation. But for field sites, such data are frequently unavailable. We herein recommend the hydrological routing approach for field applications because there is only a small difference between it and the blended hydrological-hydraulic one, and this reduces substantially data requirements; no discharge data are needed. For field application, vegetation may need to be filtered or removed from the point clouds before constructing DEMs. Vegetation is dynamic and its density (amount and size of the leaves and branches) changes through the year, which ultimately will cause DEM errors. A promising alternative would be provided by multi-echo laser scan and interferometric synthetic aperture radar (InSAR). These technologies can penetrate vegetation and produce high accuracy DEMs.

| CONCLUSIONS
In this article, we applied the morphological method in both one-and two-dimensions for a laboratory gully system. Application of the morphological method in one dimension provided a general trend of longitudinal (main-gully parallel) variation of sediment transport for the whole gully system, that is transport rates always increase from the upstream to downstream areas and become progressively stable in downstream areas in the gully system as the gullies extend headwards.
Application of the morphological method to gullies in two dimensions is facilitated by the growing ease of repeat measurement of gully morphology and the fact that sediment transport upslope of gully heads has to be zero. The latter means that the sediment supply condition that is a challenge for rivers does not apply in this case. It is also facilitated by being able to apply a no negative transport condition to calibrate and to evaluate the approach: eroded sediment, after routing, must be sufficient to sustain deposition measured downstream.
The 2D approach in gullies needs a solution for routing sediment and here we considered blending hydrological routing for the slopes and hydraulic routing for the gully bottoms. The latter makes sense physically as with the reduced slopes of a gully bottom hydraulically-driven flow forcing as well as topographic forcing may be important. We found that degree of negative transport condition was relatively insensitive to whether or not a blended hydrological-hydraulic routing was used and that results from applying the hydrological routing throughout were not significantly degraded. This is an important finding because application of hydrological routing is significantly more straightforward than the blended solution.
The 2D application of the morphological method still requires critical parameters. In the hydrological-routing based application, we used a Holmgren (1994) routing and this has the parameter α that controls the degree of diffusion. Equally, knowledge is required of bulk density. We showed that consideration should be given to DEM spatial resolution, even if this is becoming less of an issue with the growth of new DEM acquisition methodologies, and also temporal resolution. However, the framework that we have developed might be an additional tool for both experimental and field quantification of the spatial patterns of sediment transport in gullies; and quantification of how these patterns change under different forcing factors.