Sediment Modelling of a Catchment to Determine Medium-term Erosional Trends

1. Introduction

Soil erosion in tropical countries is triggered due to geologically weathered soils, intensive rainfall, inappropriate soil and land management practices, and mining activities. Where mining activities cause disturbance in catchment stream ecosystems, they are impacted by heavy loads of eroded material entering the streams [1] if not managed correctly. Transport of hydrophobic pollutants and heavy metals is often facilitated by suspended sediment particles eroded and taken down-stream during heavy discharge events [2]. Thus, suspended sediments correspond to particle-related pollutant transport from the catchment across the river reaches.

Stream-suspended mud also acts as an indicator of post-mining landform stability. During landform disturbance such as mining, fine suspended sediment (the silt and clay components of suspended sediment - FSS) spikes have been observed in the receiving waters following the disturbance. As the landform is rehabilitated and reaches equilibrium with the surrounding catchment, FSS spikes following rainfall events return to pre-mining levels [3]. The magnitude of FSS spikes following a known rainfall event can be a tool to measure landform stability [3]. Since FSS is highly mobile, many contaminants released by the catchment disturbance attach themselves to FSS and are transported downstream [4]. Thus, erosional trends of a catchment imply associated sediment and pollutant transfer within the river reaches.

Quantification of stream sediments at the catchment outlet using gauging stations is a time-consuming process. It involves collecting water samples at high-discharge events which are then dried, sieved, and weighed to quantify sediments of different particle sizes. However, there is a strong correlation between turbidity (NTU) and mud concentration (silt and clay component of suspended sediment) [5]. Thus, continuous turbidity (NTU) data collected with turbidimeters at the catchment outlet of a stream can be used to measure changes in FSS load exiting the catchment.

Continuous discharge and FSS data measured at the catchment outlet facilitates catchment hydrology and erosion modeling. Modelling also requires knowledge of model input parameters that depict the catchment characteristics, such as transpiration and surface roughness. This provides an avenue to comment on dynamics of landform stability through simulated erosional trends and catchment erosion quantification based on rainfall events. Hydrology and sediment modelling of the catchment enables simulations by the model for projected future rainfall events. This provides medium-term and long-term erosional trends of the catchment under different rainfall scenarios.

The broader objective of the project is to validate an approach using an event-based stream FSS discharge relationship in combination with a landform evolution model (LEM) to determine when erosion of a rehabilitated landform is at equilibrium with the surrounding catchment [6]. A rainfall event-based stream FSS/discharge relationship was developed using continuous receiving stream monitoring data in which a catchment disturbance caused a change in this relationship [3]. To evaluate the change in event stream FSS on account of a probable disturbance due to mining in the catchment, CAESAR-Lisflood LEM is calibrated and validated for the mining catchment. It is then run for future rainfall simulations to show changes over time and changes in response to differences in those rainfall scenarios [7]. It was found that CAESAR Lisflood accurately predicts event FSS for a specific discharge for early wet season but underpredicts for the later wet season. To run the LEM for a 1000-year period and to investigate how the event FSS/discharge relationship changes across the years, upstream continuous stream discharge and sediment input to the catchment is required. Therefore, the upstream catchment needs to be calibrated and modelled with the help of HEC-HMS to simulate future continuous discharge and sediment input to the lower catchment where the mine resides. The aim of this study was to use HEC-HMS to determine sediment and discharge quantity upstream of a catchment with mining to determine how those upstream inputs affect the downstream development of the catchment landform, in terms of sediment changes and geomorphology.

This study focusses on hydrology and sediment modeling of a catchment with HEC-HMS (Hydrologic Engineering Center-Hydrologic Modeling System). The Hydrologic Modeling System (HEC-HMS) is designed to simulate the complete hydrologic processes and soil erosion of dendritic watershed systems [8]. The application is used to model the catchment by adjusting the parameters of the model such that an input like an observed rainfall event will give discharge and sediment output at the catchment outlet similar to that measured at the site. In this study, HEC-HMS was calibrated for two individual water years and the same model parameters were used to validate the model for a third water year. Observed discharge, fine suspended quantities and rainfall used for calibration and validation of the model were obtained from gauging stations. Rainfall for future simulations was computed with the help of a weather generator incorporating climate change factors. The model was then run for a time period of 50 years to explore the medium-term erosional trends of the catchment. The application also facilitates quantification of sediment mobilized and eroded from the catchment for different wetter and drier rainfall scenarios.

Previous studies indicate that the parameters that influence the erosion rate likely to occur on the landform will change in time. A study that focused on temporal trends in hydrology and erosion for the post-mining landform at Ranger mine, Alligator Rivers Region, Northern Territory, Australia, indicates that the amount of sediment transported from the rehabilitated landform at Ranger will decrease with time, particularly during the first 50 years after rehabilitation [9]. The SIBERIA Landform Evolution Model (LEM) was used to predict medium-term erosional characteristics for another site in the same region, the Scinto 6 site, and results shows that erosion in gullies reduces over time [10]. The catchment under consideration in this study for more detailed modeling of discharge and sediment quantification is the upper Gulungul Catchment in the Alligator Rivers Region.

2. Materials and Methods

2.1. Study Site

Ranger Mine (12°41’ S, 132°55’ E) is an open cut uranium mine operated by Energy Resources of Australia Ltd. (ERA) in Northern Territory, Australia. The mine is located 8 km east of the township of Jabiru in the Alligator River Region (ARR) and within the 78 km² Ranger Project Area which is surrounded by, but separate from, the World Heritage-listed Kakadu National Park. It has been producing uranium oxide (U₃O₈) since 1981 and is presently undergoing rehabilitation and closure. Mining ceased in 2012, milling of the mined ore ceased in early 2021 and work is now focused on rehabilitation [11]. The regional geology, climate and the location significance are detailed in [7]. ERA is responsible for the rehabilitation of Ranger Mine based on laid out principles called Environmental Requirements (ERs). The Australian Commonwealth Government Supervising Scientist has a supervisory role. ERs pertaining to erosion equilibrium of the landform require erosion characteristics which, as far as can reasonably be achieved, do not vary significantly from those of comparable landforms in surrounding undisturbed areas. It is expected that the erosion rates will initially be high then tend slowly towards the natural rates. These timeframes are expected to be quite long so the outcome is to use the best available modelling to demonstrate that the erosion characteristics of the final landform will eventually be comparable to natural landscapes.

The Ranger Mine is adjacent to Magela Creek, a left-bank tributary of the East Alligator River [11]. Gulungul Creek is a small tributary of Magela Creek that is adjacent to the tailings dam. It is one of the tributaries that would be the first to receive sediment generated from the mine site during and after rehabilitation[12]. Gulungul and Magela Creeks are ephemeral, braided, sand-bed streams which carry very large sand loads (bed and suspended bed) and small FSS (fine suspended sediment) loads. The Environmental Research Institute of the Supervising Scientist (eriss) monitored sites in Gulungul Creek downstream (GCDS), Gulungul Creek upstream (GCUS), Magela Creek downstream (MCDS), and Magela Creek upstream (MCUS) of the mine (Figure 1). The data obtained from eriss for this study were continuous discharge (m³s^-1), and turbidity (NTU) and FSS (<63 µm fraction of sediment samples collected in the auto-samplers; mg). Data were from August 2004 to August 2015 measured at a frequency of 6 min for Gulungul Creek GCDS and GCUS. Rainfall data for GCDS and GCUS at 10 min intervals were also obtained for modelling Gulungul catchment in CAESAR-Lisflood.

2.2. HEC-HMS Model

HEC-HMS is designed to simulate complete hydrological processes of dendritic watershed systems. The software includes many traditional hydrologic analysis procedures such as event infiltration, unit hydrographs, and hydrologic routing. It also includes procedures necessary for continuous simulation including evapo-transpiration, snowmelt, and soil moisture accounting. The physical representation of a watershed is accomplished with a basin model. Hydrologic elements are connected in a dendritic network to simulate runoff processes. Available elements include sub basins, reaches, junctions, reservoirs, diversions, sources, and sinks. Computation proceeds from upstream elements in a downstream direction [8]. Modeling of the Upper Gulungul Catchment involved the following components:

• Precipitation methods which can describe an observed (historical) precipitation event or a frequency-based hypothetical precipitation event.
• Evapotranspiration methods
• Loss methods which can estimate the amount of precipitation that infiltrates from the land surface into the soil. By implication, the precipitation that does not infiltrate becomes surface runoff.
• Direct runoff methods that describe overland flow, storage, and energy losses as water runs off a watershed and into the stream channels generally called transform methods.
• Baseflow methods that estimate the amount of infiltrated water returning to the channel.
• Hydrologic routing methods that account for storage and energy flux as water moves through stream channels.
• Models of naturally occurring confluences (junctions) and bifurcations (diversions)
• Basin and stream erosion methods [14].

The hydrology, soil erosion and sediment routing capabilities of HEC-HMS were applied to the Upper Gulungul catchment (UGC) (Figure 2) an area of 39.1 km² which drains to an imaginary outlet at GCUS. Thus, the catchment modelling is done such that the continuous discharge and sediment output of UGS is similar to the observed values at gauging station GCUS (Figure 2). For this study, a 10m resolution DEM of UGC (Figure 3) was used in HEC-HMS and the watershed was divided into 77 sub-watersheds/basins and 38 reaches for modelling. The sub-basins and reaches were delineated by HEC-HMS based on DEM data.

2.2.1. Parameters

The parameters used for modelling are tabulated in Table 1 below and a detailed description of the same follows.

Rainwater is infiltrated until rainfall exceeds the infiltration capacity of the soil after which excess precipitation formed on the surface that contributes to the run-off. Deficit and constant loss method was used to compute precipitation losses in this study since this allows for continuous simulation when used with a canopy method that will extract water from soil in response to potential evapotranspiration computed in meteorological model [18]. Since the watershed in the study has distinct wet and dry seasons and model calibration is done for water years from the start of wet season, initial storage was 0% and initial deficit equals maximum storage in deficit and constant loss method. The initial deficit, constant loss and imperviousness for basins were attained from ARR data hub and calibrated accordingly [19].The calibrated value of initial deficit was 175mm and constant loss rate varied from 0.175 to 0.225 mm/hr. The evapotranspiration value for the catchment was taken as 5mm/day [20].

The actual surface runoff calculations for a sub basin element are performed by a Transform method. SCS (Soil Conservation Service) unit hydrograph was used as a transform method in this study. The unit hydrograph is a technique for modeling the transformation of excess precipitation to runoff at the watershed scale. The SCS unit hydrograph method defines a curvilinear unit hydrograph by first setting the percentage of the unit runoff that occurs before the peak flow. The percentage of runoff occurring before the peak is reflected in the Peak Rate Factor (PRF). The default unit hydrograph has a PRF of 484 which is used in this study [21]. The catchment lag time is calculated from equation 1[15].

Lag Time(min) = Time of Concentration/1.67,
Time of Concentration(min) = (107nL^0.333)/S^0.2

(1)

where, n= roughness value (n=0.035 for bushland)

L=overland sheet flow path length (m)
S=slope of surface (%) = [m/km]/10

L and S are attained from sub-basin characteristics in HEC-HMS

Recession Baseflow method is used to compute baseflow contributions to sub-basin outflow. The two constants, namely Recession constant and Ration to peak constant used in the method can vary from 0 to 1. The calibrated constant values for this study were 0.8 and 0.3 respectively.

Surface erosion can be computed at sub-basin elements using the MUSLE approach. It computes the total sediment load transported out of the sub-basin. Sediment yield using MUSLE is calculated using equation 2 [22].

Sediment yield (tons) = 95(Q_surf×q_peak)0.56×K×LS×C×P

(2)

where Q_surf is the surface runoff volume (acre feet), Q_peak is the peak runoff rate (cubic feet per second), K is the soil erodibility factor, LS is the topographic factor, C is the cover and management factor, and P is the support practice factor. The erodibility factor K describes the difficulty of eroding the soil. The topographic factor LS describes the surface's susceptibility to erosion due to length and slope. The practice factor P describes the effect of specific soil conservation practices. The cover factor describes the influence of plant cover on surface erosion. The values for K, LS, C and P for the watershed were determined from CSIRO (Commonwealth Scientific and Industrial Research Organisation) data raster maps [16,23] . The threshold factor is the peak flow less than which there will be no erosion or sediment yield. The calibrated value of threshold factor in this study was 1 m³/s. The exponent is used to distribute the sediment load into a time-series sedigraph which was 0.56 in this study. Gradation Curve is the sub-basin sediment particle size distribution of the watershed area where the catchment belongs [24].

The transport potential of the watershed specifies which method to use to calculate the stream flow sediment carrying capacity for non-cohesive sediments. The transport potential method used in this study was the Wilcock and Crow method [25] and the fall velocity method which determines the time required for the deposition of sediment from water column to the stream was the Van Rijn method [26].

While a reach element in HEC-HMS conceptually represents a segment of stream or river, the actual calculations are performed by a Routing Method contained within the reach. The routing method used in this study was Muskingum Cunge. The initial type of method was selected as assuming that the initial inflow is the same as the initial outflow to the reach from upstream elements. The length and slope of the reach were taken from reach characteristics in HEC-HMS. The Manning’s n roughness coefficient value was a calibrated value of 0.07. For the space-time method, Auto DX Auto DT Method was selected so that the program will automatically select space and time intervals that maintain numeric stability. The Index flow was taken as the average value of hydrograph under consideration. The reach cross section shape was given as trapezoid and the corresponding parameters of width and side slope were attained from Digital Elevation Model (DEM) reach profiles using GIS (Geographic Information System).

Sediment processes within a reach are directly linked to the capacity of the stream flow to carry eroded soil. The volume ratio method links the transport of sediment to the transport of flow in the reach using a conceptual approach. The bed width was the width of the sediment bed which is used in computing the volume of the upper and lower layers of the bed model. Bed depth is the total depth of the upper and lower layers of the bed. The active bed factor is used to calculate the depth of the upper layer of the bed model. The values for bed width, bed depth and active bed factor for the catchment were taken as 10m, 0.2m and 1 respectively [27]. The initial bed curve defines the distribution of the bed sediment by grain size at the beginning of the simulation. The computed particle size distribution of different grain sizes at GCUS was used as input [7].

2.3. Rainfall

The rainfall gauging station at GCUS (figure 1) gives continuous rainfall data at 10 min intervals. Since GCUS lies at the outlet of the Upper Gulungul catchment, using the same data to account for daily rainfall in the whole catchment was not accurate for calibration and validation of the HEC-HMS model. Thus, daily rainfall data were attained for water years 2012-13, 2013-14 and 2014-15 for modeling purposes from SILO which is a database of Australian climate data from 1889 to the present [28]. SILO gives rainfall point data at 5 km intervals. SILO data from points available across and near the catchment were considered and Voronoi polygons were drawn in QGIS [29] to determine the area of the catchment influenced by each rainfall point data [7]. Daily rainfall data from six such data points were sorted to determine the rainfall that influences respective sub-basins of the catchment.

Area Reduction Factors (ARF) were then applied to the rainfall data. The point data gives rainfall intensity at a single point. For large catchments, the rainfall intensity is always not constant across the whole of the catchment during a storm. To account for this, an ARF is included to estimate the average rainfall over the whole of a large catchment [30]. The equation for ARF is given in equation 3. Computed ARF value of 0.89 for the catchment was then multiplied with the point data rainfall intensity available at six SILO points to account for the whole catchment. This data was used as rainfall input for respective sub-basins in the model for calibration and validation.

ARF=[1−a(Area^b−clog₁₀Duration)Duration^−d+eArea^fDuration^g(0.3+log₁₀AEP)
+h10^{iAreaDuration/1440}(0.3+log₁₀AEP)]

(3)

where a, b, c, d, e, f, g, h and i are ARF parameters specific to the region under study [19].

For future simulations of over 50 years, simulated rainfall data needs to be used as input with adequate climate change factors incorporated. The Intergovernmental Panel on Climate Change (IPCC) is the United Nations body for assessing the science related to climate change [31]. CSIRO and the Australian Bureau of Meteorology (BOM) prepares tailored climate change projection reports for Australia which provides guidance on climate changes that needs to be considered in planning and projections [32]. The projections are based on the outputs from the ensemble of model simulations brought together for the Coupled Model Inter-comparison Project phase 5(CMIP6) [33]. Coupled Modeling has facilitated using eight different models for climate projections in four greenhouse gas emission criteria [34]. It provides future rainfall dataset available for discrete 30-year future periods and not continuous time series data over 100 years. Each future time series data is regarded as representative of the mean state of future climate for specified time period. Multiple future time series cannot be joined to make longer continuous time series.

Future rainfall simulations for this study shown in Figure 4 are taken from a project on Tom’s Gully which is 100 km away from Ranger mine [35]. The study uses WGEN weather simulator [36] to generate 1000 year continuous daily rainfall series from year 2025. WGEN is a stochastic weather generator which uses monthly and annual statistics to generate daily time series of precipitation, minimum temperature, maximum temperature, and solar radiation. In the absence of on-site observation data, long-term daily climate data from 1957 to 2019 was sourced from the SILO Australian climate online database and used as reference in WGEN.

3. Results

3.3.1. Model Calibration & Validation

The HEC-HMS model was calibrated with parameters given in Table 1 and described in the previous section. Calibration was done for two water years (31^st August to 31^st August) of 2011-12 and 2012-13. Daily rainfall point data from SILO improvised with area reduction factor was used as input for respective sub-basins. The model was run for each water year to obtain daily output of discharge (m³s^-1) and sediment (tonnes) at model output which is GCUS (Figure 3). The continuous discharge and turbidity data at 6 min intervals obtained from eriss from their monitoring station at GCUS was used to compare the simulated output with the observed site data (Figure 6a & 6b). Observed discharge at 6 min intervals was averaged to hourly discharge and turbidity data(NTU) at GCUS was converted to continuous FSS hourly data in tons (Table A1 in Appendix A) using equation 3 [3].

FSS=0.52∗NTU

(2)

Figures 5 and 6 thus compare the calibrated model output to the observed data for discharge and FSS respectively. The efficiency of the calibrated model to match the observed values is determined by Nash-Sutcliffe efficiency (NSE) number generated at the end of each simulation. The NSE value is computed as in equation 4.

$$\mathrm{N}\mathrm{S}\mathrm{E}=1.0-{\sum }_{t=1}^T\frac{({yt-ft)}^2}{\sum _{t=1}^T{(yt-y)}^2}$$

(2)

where y_t is the observed data values for time period t, f_t is the simulated data values for the same period, y is the mean observed data values per time period and T is the number of time periods [37]. The maximum NSE value possible is 1.0 and occurs if simulated values perfectly match observed values. Larger NSE values denote better model performance. The NSE values for 2011-12 and 2012-13 hydrology calibration in this study are 0.5 and 0.6 respectively. Since the turbidity values give the quantity of FSS using equation 3 which is the combined clay and silt component, the clay and silt quantities from the model output are summed up and compared to the observed values in Figure 6.

Since the model performs well for each calibration simulation, it was run for another water year of 2013-14 with the same model parameters used for the calibration study, for model validation. Comparison between observed values and the simulated output for hydrology and FSS are shown in Figure 5a and 5b respectively. The NSE value for hydrology simulation was found to be 0.7 which was a good fit.

3.3.2. Medium Term Erosional Trends

Long-term rainfall simulations obtained from WGEN were used to run the calibrated and validated model for 50 years. The rainfall generated had a uniform average value across the years as shown in Figure 4. The corresponding simulated discharge at GCUS for 50 years is shown in Figure 6a which also exhibits a uniform trend. Unlike the rainfall and discharge, the simulated FSS quantities across the years decrease gradually (Figure 6b). Thus, the FSS output of the Upper Gulungul Catchment decreases over the years for a corresponding uniform rainfall event in UGC and corresponding discharge at GCUS.

3.3.3. Sediment Quantification and Total Denudation Rate

The FSS sediment output for each water year from HEC-HMS simulations were 296 tonnes, 172 tonnes and 287 tonnes at GCUS for 2011-12, 2012-13 and 2013-14 respectively. The total denudation rate of the Upper Gulungul catchment after 50 years of simulations was found to total 1.2269 mm. This accounts to 0.0245 mm year^-1 of denudation rate in the catchment as predicted by HEC-HMS. The long term average total denudation rate calculated from the available observed data at lowlands near Ranger mine is 0.075 ± 0.013 mm year^-1[38]. As per CSIRO data, the denudation rates in the upper Gulungul catchment varies from 0.576 tonnes ha-1 year^-1 to 19.923 tonnes ha^-1year^-1 for different basins [16,23]. This equates to 0.0217mm year^-1 to 0.727mm year^-1 for the denudation rate of basins across the catchment. Annual denudation rate predicted by HEC-HMS from 50-year simulations is 0.0245mm which is thus within the bounds of these estimates.

4. Conclusion

This study showed that HEC-HMS was effectively calibrated and validated as a catchment model for observed rainfall and it quantified discharge and sediment output from the catchment over 50 years. The model output proved to be a good fit to the observed data. The output from the HEC-HMS model showed that the erosional trends tended to decrease over time which is in concordance with the previous studies which states that parameters that influence the erosion characteristics of a post-mining landform decreases over time, especially in the first 50 years after rehabilitation. The annual denudation rate of the study catchment predicted by the model after the 50-year simulation is in the same order of magnitude of the natural denudation rate of the area calculated in previous studies. It also falls in the range of denudation rates determined by CSIRO for the catchment area. Thus, the methodology and selection of parameters for calibration and validation of the model in this study seems adequate to predict discharge and sediment output from the catchment in concordance with the field observed data and thus enabling it to predict future erosional trends.

Previous studies on the Ranger Mine site have evaluated landform stability and temporal changes in erosional trends based on gully erosion and incision depths. This study enables modelling of continuous discharge and FSS output from the catchment and enables the evaluation of event FSS output for a specific discharge event. Thus, it enables the use of the output for evaluating FSS following specific rainfall events as an indicator of disturbance and achieving equilibrium. While this was previously carried out within a catchment [3] the use of HEC-HMS enables upstream inputs to be included. HEC-HMS also allows hourly rainfall inputs, and easily enables differing wet and dry season rainfall scenarios to be included. This was important for its use in tropical regions with highly seasonal rainfall conditions. In contrast the CAESAR-Lisflood LEM was modelled [7] but while the annual simulations gave event FSS output matching the observed site values for the early wet season, they underpredicted the FSS output for the later part of the wet season. HEC-HMS modelling simulations of FSS gave uniform comparison to observed values throughout the wet season. HEC-HMS also runs faster than CAESAR-Lisflood to predict event discharge and FSS quantities. Moreover, HEC-HMS predicted output can be used to run CAESAR-Lisflood LEM for the mine site to determine how the mine site affects the downstream development of the catchment landform, in terms of sediment changes and geomorphology.