Taking theory to the field: streamflow generation mechanisms in an intermittent, Mediterranean catchment

. Understanding streamflow dynamics at the catchment scale remains an arduous task; this is especially true for non-perennial networks. While modelling tools offer important advantages to study streamflow dynamics, the highly nonlinear, unsaturated dynamics associated with the transitions between wetting and drying in non-perennial systems make modelling cumbersome. This has stifled previous modelling attempts and alludes to why there is still a knowledge gap. In this study, we first construct a conceptual model of the physical processes of streamflow generation in an intermittent river system in South 10 Australia, based on the hypothesis that the vertical and longitudinal soil heterogeneity and topography in a basin control short-term (fast flows), seasonal (slow flow), and a mixture of these two. We then construct and parameterise a fully integrated surface-subsurface hydrologic model to examine patterns and mechanisms of streamflow generation within the catchment. A set of scenarios are explored to understand the influences of topography and soil heterogeneity across the catchment. The results showed distinct flow generation mechanisms develop in the three conceptualized areas with marked soil and 15 topographic characteristics, and suggested that capturing the overall distribution of soil types across the catchment was more important than capturing the wide variability of soil hydraulic properties. This study augments our understanding of catchment scale streamflow generation processes, while also providing insight on the challenges of implementing physically-based, integrated surface-subsurface hydrological models in non-perennial stream catchments. hundreds of similar non-perennial river systems in the semi-arid, coastal Mediterranean climate regions of South and Western Australia, California, South Africa, and around the Mediterranean itself. This study provides an initial step towards understanding non-perennial streamflow generation processes at the catchment scale and provides a template for using ISSHM’s for process understanding in these stream systems. The development of a conceptual model of the most important factors impacting flow generation processes within Pedler catchment presented a hypothesis that combined our understanding of field data with lessons learned from previous studies. This conceptualization was key to inform the model setup and to capture dynamics of streamflow generation in this non-perennial stream system. Difficulties in setting up and running the model reaffirmed the numerical difficulties experienced in large scale unsaturated models, as accurately reproducing the topography and observed initial conditions was a challenge and model runtimes prohibited extensive exploration of multiple scenarios. In particular, the importance of preserving channel representation to model streamflow generation on non-perennial


Introduction
With water scarcity increasing globallyIn a world of increasing water scarcity, understanding the hydrology of rivers in arid and semi-arid regions has become increasingly especially important. In these regions, most streams and rivers are nonperennial, meaning surface flow ceases for some or most of the year. The understanding of processes that lead to streamflow 25 generation in non-perennial rivers is incompleteIt is widely recognised that our hydrologic understanding of these nonperennial streams is still insufficient (Costigan et al., 2017). In particular, our understanding of the processes that lead to streamflow generation in non-perennial rivers is incomplete (Costigan et al., 2017;Gutierrez-Jurado et al., 2019;Shanafield et al., in press2021). This is partly due to lack of appropriate data; indeed, it has been shown that streamflow gauges, which provide the most fundamental data to understanding the hydrology of rivers, are preferentially located on perennial rivers 30 (Fekete and Vörösmarty, 2002;Poff et al., 2006). In addition, the particular challenges associated with characterising unsaturated flow and highly transient streamflow, the natural wetting and drying cycles, and other highly non-linear processes complicate the efforts to model non-perennial systems (Beven, 2002;Ye et al., 1997).
Among semi-arid regions, Mediterranean climate regions are relatively well-represented in the literature, because they are widespread globally, as unlike inland regions, they have long been population centers and therefore have experienced severe 35 anthropogenic alteration, and are experiencing ever increasing anthropogenic water demands (Merheb et al., 2016). This type of subtropical, dry summer, wet winter climate exists not only in southern Europe, but also represents much of California, South and Western coastal Australia, South Africa, and central coastal Chile (Geiger, 1961). Intermittent rivers (i.e. those that flow seasonally) and ephemeral rivers (i.e. those that flow only after rain events) are already prevalent in these regions.; however, Understanding streamflow generation mechanisms in these rivers is particularly needed Mediterranean climate 40 regionsbecause Mediterranean climate regions are sensitive to climate change (Cudennec et al., 2007) and expected to experience significant drying due to shifting climate patterns, which will greatly impact streamflow regimes (Milly et al., 2005). Therefore, understanding of streamflow generation mechanisms in these regions is particularly neededNumerical models offer a methodare used to explore the complex drivers that lead to streamflow production in a catchment. Several modelling studies of non-perennial river catchments have provided insight on the role of vegetation, soil cover, topography, 45 antecedent wetness, and soil heterogeneity on runoff generation (Carr et al., 2014;Pierni et al., 2014;Ebel et al., 2016;Maxwell & Kollet, 2008;Vivoni et al., 2007) and on the evolution of saturated area patterns (Weill et al., 2013), as well as the importance of unsaturated storage dynamics as major controls on the processes of runoff generation (Vanderkwaak & Loague 2001;Gutierrez-Jurado et al. 2019;Mirus & Loague 2013 [Table 1]). Nevertheless, the required level of information to adequately parameterize boundary value problems have restricted the use of fully integrated surface-subsurface hydrologic 50 models (ISSHM) in non-perennial river catchments to mostly small-scale hillslope or headwater catchments (0.001-0.9 km 2 ).
Moreover, most ofa majority of these the numerical models studiesfor runoff and streamflow generation in non-perennial rivers captured investigated relatively short time periods ranging from 2h to 330 days and . Under these approaches, the aim is to investigate reported the integrated system response to a set of distinctive scenarios and/or conditions (Carr et al., 2014;Di Giammarco et al., 1996;Heppner et al., 2007;Kollet et al., 2017;Mirus et al., 2009;Panday and Huyakorn, 2004). In contrast, 55 there is still a knowledge gap regarding rather than to understand longer-term hydrological processes that controls on the drywet transition.
Given these limitations, tThe overarching, spatiotemporal processes that control key catchment-level dynamics in nonperennial rivers remains a knowledge gap. For instance, in small-scale systems, the hydrological processes occurring at a given time and place (i.e. 'active' processes [Ambroise, 2004]) might be the same as those contributing to flow generation at that 60 same time. However, in large-scale systems the hydrologic response is influenced by different surface water-groundwater travel times, initial losses (e.g. evapotranspiration or infiltration) and the connectivity of the areas where hydrological processes are occurring. Consequently, the active processes occurring at a time and place do not necessarily contribute to the integrated catchment response at another given point at that or a later time. This is particularly important for non-perennial rivers, where In contrast, a defining characteristic of the dry-wet transition in non-perennial rivers is that the dry initial conditions exacerbate 65 initial runoff and streamflow losses due to high rates of streambed infiltration, causing the development of saturated areas and the generation of runoff and streamflow to occur discontinuously throughout the catchment.
Within the context of Mediterranean climate catchments, Gutierrez-Jurado et al. (2019) used an ISSHM in an idealized concept development study to begin exploreing the processes leading to the transition from dry streambed to flowing stream. This theoretical study concluded that soil hydraulic properties and unsaturated storage dynamics exhibit strong control over 70 streamflow generation and determine the spatiotemporal development of runoff generating areas and dominant flow generation mechanisms. It also highlighted the importance of understanding the development and progression of active areas (i.e. where processes are active) and their dominant flow generation mechanisms to understand the pathways and threshold of streamflow generation. But how applicable are the findings of this small-scale, idealised and simplified model to real, complex, and largerscale catchments? 75 The goal of this study is to investigate the extension of an ISSHM to a mid-sized, real world Mediterranean climate catchment in South Australia with a stream network composed of ephemeral and intermittent reaches. We hypothesise that topography, groundwater level, and soil properties are also the dominant controls over streamflow generation mechanisms in mid-sized, Mediterranean climate coastal catchments such as those typical of South Australia. We first present a conceptual model of the potential flow mechanisms for a representative catchment, based on the theoretical findings of Gutierrez-Jurado et al. (2019), 80 information of shallow soil profiles, and hydrological, meteorological, and geologic data. This conceptual model is then tested using an ISSHM, which is needed to fully capture the physical properties of interest. Given the inherent difficulties of modelling a large, unsaturated domain with contrasting soil layers during state changes (dry to wet), the goal of this study was not to reproduce the field observations, which in non-perennial rivers are strongly a function of antecedent moisture conditions, but to more broadly understand the interplay between shallow soil propertiesprofiles, groundwater levels found throughout the 85 catchment. Finally, through simulations across the range of likely soil properties and with varying degrees of model discretisation, we explore whether the general lack of detailed soil data at catchment scale and computational difficulties in capturing localised variation in stream geometry impact characterisation of streamflow generation mechanisms. WWe first present a conceptual model to develop hypotheses of the potential physical flow processes based on the theoretical findings of Gutierrez-Jurado et al. (2019), information of shallow soil profiles, and hydrological, meteorological, and geologic data. This 90 conceptual model combines field data with our understanding of the roles of the shallow soils, geology, topography, and water levels in the catchment. Then, we use available shallow soil distribution and hydrological, meteorological, and geologic data indevelop a three-dimensional, HydroGeoSphere (Therrien et al., 2010) model of the catchment, coupled with the Hydraulic Mixing-Cell method (HMC) developed by Partington et al. (2011). The HMC method tracks rainfall as it moves through the catchment, allowing the identification of active areas and the quantification of the contributing flow generation mechanisms 95 on those areas. Given the inherent difficulties of modelling a large, unsaturated domain with contrasting soil layers to capture the physical processes at a rather sudden state change (dry to wet), the goal of this study is not to reproduce the field observations, which in non-perennial rivers are strongly a function of antecedent moisture conditions, but to test and re-evaluate our conceptual understanding of the physical processes occurring in the catchment. Both the insights gained from this modelling process and the challenges encountered in attempting to accurately reproduce the shallow soil distribution, climate, 100 geology, and hydrology of such a dynamic system are discussed.

Study Area
The catchment used for this model is Pedler Creek which is part of the larger Willunga basin located roughly 30 km south of Adelaide, South Australia ( Figure Fig. 1 an area of high agricultural value mainly for the viticulture industry. The total catchment area is approximately 107 km 2 and discharges into the sea on the Gulf of St Vincent to the west of the catchment. A wastewater treatment plant located in the town of McLaren Vale discharges water into the creek; therefore, we only considered the 69 km 2 area of the catchment upstream of the treatment plant (this represents 80% of the stream network). area upstream of the creek before it passes through the town for this study. The Pedler sub-basin has an area of 69 km 2 , encompassing over 60% of the catchment and more 110 importantly, it contains roughly 80% of the total length of the stream network. From the 4 major land uses in the catchment, aAgriculture (45%-30.4 km 2 ) and grazing (46%-31.8km 2 ) dominate the over 90% of the landscapeland use, with small urban and plantation forestry areas. Urban (residential and commercial) account only for 8% (5.4 km 2 ), and plantation forestry covers less than 1% (0.38 km 2 ) also dotting of the catchment ( Figure Fig. 1). The creek flow regime is intermittent generally flows continuously from July to September in response to the winter rains, with isolated and ephemeral flows during the rest of the 115 year, flowing only after extreme rainfall events. For the (daily) period of record 2000-2018 at (gauge ID A5030543, Department of Environment and Water, Government of South Australia), the creek flows on average 120 days per year, ranging from 33 to 199 d yr -1 . Mean annual flow volume is 3.88 x 10 6 m 3 with the higher flows occurring between July and September A wastewater treatment plant located in the town of McLaren Vale discharges water into the creek; therefore, we only considered the area upstream of the creek before it passes through the town for this study. The Pedler sub-basin has an area of 120 69 km 2 , encompassing over 60% of the catchment and more importantly, it contains roughly 80% of the total length of the stream network. The creek flow regime is intermittent flows an average of 120 days continuously during winter (33 to 199 d yr -1 over the 2000-2018 period, gauge ID A5030543, Department of Environment and Water, Government of South Australia), and only after extreme rainfall in other seasonsgenerally from July to September in response to the winter rains and ephemeral during the rest of the year, flowing only after extreme rainfall events. For the period of record 2000-2018 at gauge ID 125 A5030543, the creek flows on average 120 days per year, ranging from 33 to 199 d yr -1 . This intermittent flow represtents a Mmean annual streamflow volume ofdischarge is 3.88 x 10 6 m 3 , with most of the flow with the higher flows occurring between July and September (Water Data Services, 2019). The Mmean annual precipitation for the basin is 550 mm, ranging and ranges from 289 to 812 mm for over the period of record   in January and lower daily temperatures registered during June and July (SILO -Australian Climate Data, 2019).
As part of the Willunga basin, Pedler catchment presents a complex multi-aquifer system. The groundwater system consists of four main aquifers: Quaternary sediments, Port Willunga Formation, the Maslin Sands and the basement fractured rocks (Aldam, 1989). The Maslin Sands and the Port Willunga aquifer are separated in some locations by the Blanch Point Formation which acts as an aquitard. All hydrogeological units outcrop at the surface and regional groundwater flows towards the coast 135 from northeast to southwest.
The catchment topography consists of a low-lying coastal plain with mild undulating hills towards the north of the catchment and separated by the Willunga Fault to the steep hills located on the east of the catchment ( Figure Fig. 2). The eElevation ranges from ~400 m on the northeast of the catchment (steep hills) to ~50 m at the catchment's outlet. The hills area on the east of the fault is characterized by having a shallow sediment profile (0.5 -2 m) which is underlain by the basement rocks 140 while west of the fault the sediments thicken seaward. Surficial soil types (the upper 1.5 m) in the catchment can be clustered into three major soil groups: loam, sand and clay. Covering roughly 62% (42.36 km 2 ) of the catchment, the loam soils are distributed on the middle-eastern area where the majority (80%) of the stream network is located. ; sSandy soils cover around 32% (21.7 km 2 ) and are located mainly on the north part of the catchment with some patches present in the middle section (valley), while; the clay soils account for only 6% (3.94 km 2 ) of the catchment area and are located on the further downstream 145 section towards the west of the catchment. Most of the stream network (over 80%) is located within the loam soil. Detailed sSoil profiles obtained from the Department of Environment Water and Natural Resources (DEWNR) consistently show a distinctive clay layer starting from 1.1 to 1.5 meters depth in the sandy soil areas and at around 0.5 m within the loam areas (Department of Environment Water and Natural Resources, Government of South Australia). . Despite this general understanding of soil types, there are no available estimates of hydraulic conductivity for these shallow soils. Below these 150 shallow soils, regional groundwater flows from northeast to southwest towards the coast. The groundwater system consists of four main aquifers: Quaternary sediments, Port Willunga Formation, the Maslin Sands and the basement fractured rocks (Aldam, 1989). The Maslin Sands and the Port Willunga aquifer are separated in some locations by the Blanch Point Formation, which acts as an aquitard.
Surface water-groundwater interactions within Pedler catchment have been documented in previous studies to play a critical 155 role in the creek flow's regime. Sereda and Martin (2000) observed rapid groundwater level rises in response to large precipitation events in some shallow monitoring wells adjacent to creeks. They also noted while noting that GW level declines in the Quaternary aquifer during 1995-1999 could be attributed to a decrease of yearly precipitation during that period.
Moreover, Harrington (2002) observed that groundwater levels seemed to mirror streamflow records for the Creek. While these observations confirm that GW recharge occurs from precipitation and creek seepage, these and other studies have also 160 indicated that GW discharge occurs in some areas of the creek (Harrington 2002, Anders 2012. Further studies have indicated that the creek presents both gaining and losing stream sections, which are not only spatially but temporally variable and which are dependent on rainfall and shallow groundwater levels (Harrington 2002;Brown 2004;Irvine 2016;Anders 2012).

Conceptual Model of Streamflow Generation Process in Pedler Creek
For medium to large size catchments such as Pedler Creek, the interactions between topographic features such as slope and 165 mean soil thickness, with surficial soil heterogeneity, various aquifers properties, and a spatially variable depth to GW are likely to result in variability in streamflow generation processes developing at different spatiotemporal scales within the catchment. To understand the integrated catchment response as well as the stream network dynamics (development, expansion, and contraction) it is paramount to capture the spatiotemporal occurrence of the different streamflow generation mechanisms. Based on field data, available soil and aquifer information and the topographyUsing available soil and topography 175 information for the catchment (DWLBC, 2004;Hall et al., 2009;DEWNR, 2016), we first developed a conceptual model to outline the most likely processes leading to streamflow generation and the resulting dominant streamflow generation for Pedler Creek (Table 21). We identified three major areas with distinctive characteristics: 1) the steep hills on in the east, 2) the undulating hills on the north, and 3) the flat valley on the south-west area of the catchment ( Figure Fig. 3a). These three main areas in the map aim to provide a spatialal understanding of the most likely streamflow generating processes, with 180 different processes and areas of the catchment contributing to the short ephemeral flows during summer and late fall, the early winter build up to intermittent flow, and throughout the rainy season continuous flow. . , the ephemeral flows in Pedler Creek occur in direct response to extreme precipitation events (characterized by low-duration, high-intensity precipitation) where the temporal component can be considered negligible. Instead, ephemeral flows are mostly linked to the spatial characteristics of the catchment that favour infiltration excess overland flow. Conversely, it is paramount to include the 185 temporal component to understand how the threshold of flow is reached for intermittent flow during the rainy season. We hypothesize that the dry conditions at the beginning of the rainy season will result in most of the rainfall to infiltrate due to the high infiltration capacity of the soil (Blasch et al., 2006;Batlle-Aguilar & Cook, 2012;Mihevc et al., 2001). However, differences in topography and shallow soil profiles will promote different processes to develop as the rainy season progresses. As the ephemeral flows in Pedler Creek occur in direct response to extreme precipitation events (characterized 190 by low-duration, high-intensity precipitation), the temporal component can be considered negligible. Instead, flow is mostly linked to the spatial characteristics of the catchment where we hypothesize that the dominant flow generation mechanism is infiltration excess overland flow originating from the loam and clay soils throughout the catchment. Conversely, it is paramount to include the temporal component to understand how these processes developed and the threshold of flow is reached for intermittent flow during the rainy season. We hypothesize that the dry conditions at the beginning of the rainy 195 season will result in most of the rainfall to infiltrate due to the high infiltration capacity of the soil. However, differences in topography and soil characteristics will promote different processes to develop as the rainy season progresses. A A A detailed description of the processes for each area and the spatiotemporal development of the most likely dominant streamflow generation component for intermittent flow is provided below ( Figure Fig. 3b-d).

Steep Hills; Fast Flow 200
The steep hills are characterized by a permeable and shallow (top 0.5 m) loam soil underlain by a heavy clay profile with steep slopes (Figure3b). The combination of the shallow loam soil permeability, the high infiltration capacity, and the steep slopes are likely to allow the water to infiltrate and to flow relatively fast as unsaturated-IF towards the stream (FigureFig. 3b.1-2).
We hypothesize that the shallow loam soil profile and the water holding capacity of the loam will promote a perched GW mounding along the riverine area which will result in SE-OF from the riverine area and the adjacent hillslope to develop as the 205 dominant streamflow generation mechanism (FigureFig. 3b.3). We hypothesize that this area contributes heavily to the temporally isolated ephemeral flows that occur in direct response to extreme precipitation events (characterized by lowduration, high-intensity precipitation). Theand the dominant flow generation mechanism for these events would be infiltration excess overland flow.

Undulating Hills; Slow Flow 210
The undulating hills consist of a highly permeable deep (top 1.1 m) sandy soil profile underlain by a heavy clay layer with mild slopes (FigureFig. 3.b). We expect that Tthe high soil infiltration capacity and permeability of the sand will result in a large infiltration rate allowing most of the early winter precipitation to infiltrate in this area (FigureFig. 3.3.c.4) (Blasch et al., 2006;Batlle-Aguilar & Cook, 2012;Mihevc et al., 2001). As the infiltrated water reaches the low permeable clay layer it will move in the subsurface as IF towards the low-gradient areas (FigureFig. 3c.5). We hypothesize that the high infiltration rates 215 in combination with the mild slopes (or low-gradient areas in the valley) will favour the development of a perched GW that will rise uniformly, allowing the river to develop into a gaining condition. As the infiltrated water moves as IF it will discharge into the downstream areas (FigureFig. 3c.6). Due to the larger unsaturated storage and the mild slopes, this area will likely take longer to contribute to flow (i.e. more water will be needed and therefore more time to reach the threshold of flow generation). We hypothesize that these areas will therefore provide the "slow flows" necessary to sustain intermittent flow for 220 the days without rainfall during the intermittent season and conversely, they are not likely to contribute to flow during ephemeral events (FigureFig. 3e).

Flat Valley; Mixed Flows
The flat valley comprises a mix of the previous two soil profiles (deep sand and shallow loam underlain by heavy clay) and a heavy clay area, all located in a low-gradient topography (FigureFig. 3d). The GW becomes shallower near the riverine areas 225 in the valley and depth to GW decreases towards the outlet area (the bore with the shallowest GW is located near the outlet where the GW is ~2 m below surface elevation). This zone has the largest draining area with both the steep and undulating hills draining towards it. The diversity of conditions in this area is likely to result in a combination of the processes previously discussed as well as additional ones. We hypothesize that the processes originating on the sandy soil areas on the valley will be similar to those on the undulating hills with the difference that the unsaturated-IF that might originate early during the 230 season might only contribute with a small amount of flow that might reflect further downstream. We also expect to see some saturated/unsaturated-IF originating early in the season in the loam areas on the valley (FigureFig. 3d.8-9). However, we hypothesize that the low-gradient terrain along with the water holding capacity of the loam soil will slow down water moving as interflow and rather promote the soil saturation to build up in the shallow soil profile. As the saturation increases, we expect the dominant streamflow generation mechanism will switch to saturation excess overland flow from both the hillslope and the 235 river area.
The clay's low permeability will limit infiltration and favour water to pond on the surface on the clay areas, which will eventually result in infiltration excess overland flow (Table 43d.10-11). The large draining area of the valley combined with the low-gradient topography is likely to promote the development of a perched GW along the riverine area which will result in SE-OF along some sections of the river (FigureFig. 3d.9). During wet years, sections of the creek near the outlet where GW 240 is shallow, are likely to develop into a gaining state with old-GW contributing to streamflow (FigureFig. 3d.10). Once the saturation threshold has been met along the riverine area in the steep hills and throughout the loam areas in the valley, SE-OF from those areas and the IE-OF from the clay are likely to contribute with the "fast flows" as travel times for overland flow are generally smaller than those for subsurface processes (FigureFig. 3.b).

Modelling Platform and HMC Method 245
To explore the potential for streamflow within the structure of the conceptual model, we built a fully integrated, numerical model of the catchment with the hillslope fraction divided into the three dominant soil types (FigureFig. 4). We used HydroGeoSphere (HGS), a 3-D fully integrated surface-subsurface hydrological model (ISSHM) that allows physically-based simulations of hydrological processes by using the control-volume finite element method to simultaneously solve the surface and subsurface flow equations. The numerical code uses the diffusion wave approximation to the Saint-Venant equations for 250 2D surface flow and a modified form of the Richards' equation to solve the variably saturated subsurface flow. Further details on the physical and mathematical conceptualisation and the implementation of the HGS code can be found in Aquanty Inc.
The decomposition of flow into the different flow generation mechanisms is provided by coupling HGS with the HMC method which is based on the modified mixing-cell method (Campana and Simpson, 1984). Using the standard hydrological output 255 from a numerical model, the HMC method allows the partition of flow in any node within the catchment. To do this partition the HMC method tags the existing water at the beginning of the simulation and any new water as it enters the model domain by area of origin (i.e. stream, hillslope, and the porous media), by and by boundary condition (i.e. the source of water). and for the case of the nNew water is tagged, by the internal model state of saturation of the area of origin (i.e. saturated or unsaturated soil profile). Using these tags, the water is tracked as it moves through the model domain and after each time step 260 of the flow simulation the method calculates the fraction of water in each cell that derives from the different flow components (Table 21)

Model Setup
A three-dimensional, HydroGeoSphere (HGS, Therrien et al., 2010) model of the catchment, coupled with the Hydraulic 265 Mixing-Cell method (HMC) developed by Partington et al. (2011) was used to test the conceptual model. The HMC method tracks rainfall as it moves through the catchment, allowing the identification of active areas and the quantification of the contributing flow generation mechanisms on those areas.

Model Discretization
The topography for the surface elevation was implemented by using a digital elevation model ( nodes and 159,192 total triangular elements. As the streams in the study area were only a few metres wide the DEM did not capture the incision of the streams into the landscape. In order to overcome this limitation, the digital elevation model was post-processed using a Python routine script to depress the elevation of the stream nodes.

Porous Media Properties 280
Hydraulic conductivity, porosity and specific storage based on literature estimates for the Quaternary Sediments, Port Willunga Formation, Blanche Point formation and the Maslin Sands were assigned using rasters of the top and bottom elevations for each unit (Aldam 1990;Anders 2012;Irvine 2016;Martin 1998Martin , 2006 (Table 32). Unsaturated hydraulic parameters were then estimated from Carsel and Parrish (1988) and Mirus et al. (2011), which had the closest hydraulic conductivity values to the estimates for each hydrogeological unit. The basement fracture rock formation is believed to act as an impervious layer 285 throughout most of their its length (Knowles et al., 2007;Martin, 1998); therefore, the elements for this layer were assumed to be inactive on the time-scale of these simulations.
The shallow soils were considered as the top 1.5 m of the subsurface domain which was the average depth reported in the soil characterisation datasheets. The information of the horizontal and vertical distribution of soils was assigned into the model using 2D overlays (horizontal) for the three main soil areas and the mesh layers generated during the grid discretization 290 (vertical). We used a digital soil-landscape map obtained from the Department for Environment, Water and Natural Resources of South Australia (DWLBC, 2004;Hall et al., 2009) to differentiate the spatial distribution of the three dominant shallow soil types. The vertical heterogeneity was determined by analysing soil characterisation datasheets from detailed soil profiles available within the Pedler sub-catchment ( [DEWNR, 2016]).

295
As quantitative soil hydraulic properties were not available, we tested a range hydraulic parameters representative of the three main soil types (sand, loam and clay) obtained from Carsel and Parrish (1988) as shown in Table 32. Data of particle size distribution at different depths (soil layers) from six soil profiles (two within the Pedler sub-catchment and four nearby) were used to estimate soil hydraulic parameters for the three dominant soil types in the catchment. We To validate the selected range, we estimated estimate soil hydraulic values from six soil profiles (two within the Pedler sub-catchment and four nearby) 300 used the ROSETTA model H2 (Schaap et al., 2002t)hat included data of particle size distribution at different depths (soil layers) usinged the ROSETTA model H2 (Schaap et al., 2002). As explained below, the influence of the hydraulic conductivity values for each shallow to select a range of values typical of each soil type to be tested. As explained below, the influence of the hydraulic conductivity values for each shallow soil, which is typically highly variable and not well known at catchment scale, was then explored using scenarios. 305

Overland Flow Properties
Manning's roughness coefficients (n) derived from Chow (1959) were implemented for the three prevalent land uses (i.e. agricultural, pasture and urban) which account for over 99.5% of the catchment area. We used the values for cultivated areas with mature row crops (4.05E-7 d/m 1/3 ) for the agricultural areas, for pasture with no bush and short grass (3.47E-7 d/m 1/3 ), and the value of asphalt (1.85E-7 d/m 1/3 ) was applied to the urban area (Table 32). For the stream network, we used values for 310 a clean and straight natural channel for the headwater sections (4.05E-7 d/m 1/3 ) and of weedy reaches for the middle-lower sections (1.15E-6 d/m 1/3 ). Rill storage height was set to 0.01 m uniformly across the domain and no obstruction storage height was implemented. The coupling length was set at 0.001 to warrant a good coupling of the surface-subsurface domains which is paramount to capture streamflow generation processes (Liggett et al., 2014).

Simulation Period and Initial Conditions 315
We selected a 4-year simulation period from January 2015 to December 2018 to ensure a representative set of years with average (2017 ~500 mm/y), below average (2015 and 2018 ~400 mm/y) and above-average (2016 ~800 mm/y) annual rainfall amounts. Precipitation records (15-minute data) from the McLaren Vale (located in the valley) and the McLaren Flat (located near the steep hills) stations (MEA, 2019) were averaged and applied as a fluid flux to the surface of the model domain. To determine the optimal time resolution for the precipitation forcing we tested preliminary models with quarterly-hour, 1-hour 320 and 24-hour inputs. Results from the preliminary models show better convergence and smaller errors in the water balance for the hourly precipitation inputs. Estimates of potential evapotranspiration (ET 0 ) were only available at a daily time step; therefore, we used values of solar radiation to approximate ET 0 at hourly intervals to match the precipitation inputs. Values of ET 0 that were less than 0.0001 m/h were considered numerical noise and were excluded from the input dataset. The resulting ET 0 dataset was applied to the surface domain. Actual evapotranspiration (ET) and interception are simulated as mechanistic 325 processes within HGS using the concepts by Kristensen and Jensen (1975) and Wigmosta et al. (1994) which require plant and soil conditions (Aquanty, 2016). Vegetation characteristics cited in the literature for Eucalyptus were used on the riverine area and values typical of grass (Banks et al., 2011;Geeroms 2009;Hingston et al., 1997) were used for the rest of the catchment (Table 32). Although a large area on the catchment consists of vineyards, during the winter months the vines are dormant without leaves and grass is commonly used as an inter-row soil cover. We did not include the effects of irrigation, ET, and 330 interception during the vines growing season as we considered the overall effects for streamflow generation would be negligible since they occur during the driest and hottest months of the year when the stream network is dry. With the simulation starting in January (the hottest month) we assumed completely dry initial conditions for the surface domain (i.e. no presence of surface water).
Initial groundwater levels were iteratively achieved by draining a fully saturated model and comparing the resulting water 335 table to field data; this process of draining required model simulations to run for days to weeks. The goal was to obtain realistic soil moisture profiles characteristic of dry, summer conditions and a smoothly varying water table surface. The resulting water levels were compared with the long-term average of data obtained from the Government of South Australia (https://www.waterconnect.sa.gov.au/) for the McLaren Vale prescribed wells area. Preference was given to matching wells with shallow GW heads (< 10 m depth) located on the flat valley area and close to streams, which are known to transition from 340 losing to gaining conditions in response to increases in the groundwater level. Finally, an output time where the simulated GW heads at four of the five wells in this area were within 1 m of average recorded levels was selected as the initial conditions for the porous media in the subsequent simulations. Initial groundwater levels in the steep hills area matched average recorded levels within 1 -28 m, where the depth to water is quite deep and varies by tens of meters across the fault line (FigureFig. 4).

Boundary Conditions 345
Boundary conditions for the outflow in the surface domain were set as a critical-depth boundary at the catchment's outlet and as a no-flow boundary condition for the rest of the domain. We applied a fluid transfer boundary condition around the catchment outlet which allows for the discharge of groundwater through the subsurface. The hydraulic gradient for the fluid transfer was given by setting a hydraulic head ~4 m below the surface elevation (the known deepest GW head near the outlet) at 10 m from the outlet faces. 350

Simulation Implementation and Data Post-processing
The simulations were performed in HGS using the control-volume finite-element mode and the dual-node approach for surfacesubsurface coupling. We used an adaptive time step with a computed under-relaxation factor scheme to aid the computational efforts. Adaptive time-stepping was applied by an initial step size of 0.001 days, a maximum step multiplier factor of 2.0, and a maximum time-step of 5 days. The simulations were run in parallel mode using 6 CPUs cores from an AMD EPYC 7551 355 @2.55Ghz (with 32 cores / 64 threads) compute node to partition the model domain. The HMC method was set up to track the flow generation mechanisms originating from the different soils in the overland areas (clay, loam, and sand), directly in the river, and from the porous media (Table 21).
We ran over 52 preliminary models testing different mesh discretisations and time resolutions for the model forcings (i.e. precipitation and ET), simulation control values, and draining simulations to try to select the optimal model setup. From the 360 final setup, we developed a final set consisting of 8 scenarios to be tested; four corresponding to sets with different combinations for the shallow soils hydraulic properties (K sat and their corresponding unsaturated storage parameters α, β, and θr); and four scenarios with different values for incising the river nodes (Table 43). Due to the computational constraints, only one set of soil hydraulic properties was used to test the scenarios with the incised stream. Results from these two sets of scenarios were used to evaluate the need to modify and test further scenarios. 365 Model output for the surface domain (2D) was post-processed to identify and quantify the activation of areas (flow onset) and the flow generation mechanisms. The output files were processed in Python to determine the HMC fraction (flow generation mechanism) that contributed most of the flow (dominant fraction) at every single node and for each output time. Results of the dominant fraction were then included as a new variable to the overland output file. A water depth threshold equal to the rill storage height (0.01 m) was used to determine when an area was considered active (i.e. values of 0.01 m or less were considered 370 as rill storage and not flow). Output for the porous media (3D) was used to support the HMC dominant fractions findings.

Results
The computational demands of modelling a large and variably saturated domain subject to sudden state changes from dry to wet conditions led to in extremely slow model convergence. From the set of scenarios testing different values for incising the stream (Table 43), only scenario 8 (incision = 10m) finished under a reasonable computation timeframewithin two months of 375 run-time. Scenarios 5-7 showed less than 10% of progress after 20 days of computation time. Nevertheless, the comparable results between scenario 8 and 4 which shared the same soil hydraulic properties but had the two ends of the spectrum in respect to the river incision (the most vs. none) suggest that results from scenarios 5-7 would have likely shown similar results.
Therefore, from here on we will only focus on the results from scenarios 1-4 (different soil hydraulic properties) and 8 (river incised). 380

Active Areas and Dominant Flow Generation Processes Determined with HMC
The development of active areas (initiation of flow) in terms of the timing and extent was similar among scenarios 1-4, while for scenario 8 (scenario with the incised stream nodes) the aerial extent was consistently smaller (FigureFig. 5). Across all scenarios, flow was generated first on areas from the steep hills and these areas expanded and contracted throughout the simulation. Fragmented active areas developed along the stream network for scenario 8 while for scenarios 1-4, the active 385 areas along the stream developed as result of the flow from the hills connecting to stream network and expanding from there.
Although the overall active areas along the river were larger for scenarios 1-4, a larger length of the stream network showed flow for scenario 8. The goals of this study were to provide insight into streamflow generation processes at the catchment scale for an intermittent, Mediterranean climate catchment to better understand the importance of controlling characteristics identified previously with small-scale theoretical modelsduring previous modelling efforts in non-perennial stream catchments (Table 1). In particular, the role of soil heterogeneity on different streamflow generation processes has been documented in previous studies. This includes soil heterogeneity, as For instance, studies in both perennial and non-perennial catchments have shown that vertical 410 soil heterogeneity can result in the development of perched saturated zones that contribute to flow generation (Hathaway et al., 2002;Maxwell and Kollet, 2008). Other studies indicated that horizontal heterogeneity contributed to the spatio-temporal variability of flow generation under different mechanisms which overall resulted in longer flow durations due to delays in runoff occurring from areas with high infiltration capacities (Ebel et al., 2016;Luce & Cundy, 1994;Smith & Hebbert, 1979).
Therefore, iIn this study we developed a conceptual framework of the hydrological processes identified in three distinct hypothesized that these distinct topographical conditions and soil types have a definite influence on streamflow generation mechanisms. However, although each catchment has its own set of conditions, most Mediterranean climate catchments would have similar topography to what was modelled here, with hills graduating to a coastal plain. Moreover, because the modelled 425 catchment included a range of soil types, we were able to explore the variation in streamflow generation processes across several soil types. Moreover, it is likely that many other seasonally-flowing catchments would have similar variation in soil, as the periodic and often flashy nature of streamflows carries fine material from the steeper headwaters and deposits it on the plains (Jaeger et al., 2017).A The results of a final set of 8 5 models comprising 4 scenarios testing different sets of soil hydraulic properties for the shallow soils and 4 1 scenarios testing the effects of testing different values for incising the stream 430 nodes were testedare discussed.
Model results overall supported our conceptual understanding that distinctive topographic and soil characteristics explain of the flow generating processes in Pedler Creek. Results from the active areas showed distinct mechanism developing in the three major areas (FigureFig. 5), supporting the idea that there is a spatial and temporal variation of flow generation processes in Pedler Creek. In the model, flow developed first at the steep hills areas (fast-flow) and the dominant mechanism was SE-435 OF with a few areas showing unsat-IF as hypothesized in our conceptual model ( FigureFig.s 3b and 5). An unexpected development was the contribution of pre-event GW during flow recessions. In this area the pre-event GW was likely to be pre-event soil water since evaluation of the model showed that the groundwater level did not rise to intersect land surface in this area. The flows generated in the valley near the outlet were similarly simulated via the conceptualized mechanisms (Fig.ures   3d and 5). We saw small areas with flow originating from IE-OF from the clay areas, and a combination of unsat-IF and pre-440 event GW for the rest of the active areas in this region. The GW did rise above the surface elevation in this area, supporting the GW contribution to flow in this area. Finally, in the few small areas close to the sandy undulating hills region, flow was simulated through the unsat-IF mechanism as predicted in our conceptual framework ( FigureFig.s 3c and 5). These results support the findings by Gutierrez-Jurado et al. (2019), who suggested that soil properties largely dictate the dominant flow generation mechanisms and that unsaturated storage dynamics control the thresholds and pathways of flow. For real catchments 445 such as Pedler Creek, soil properties and topography evolve in tandem and it is impossible to fully disentangle their relative influence on streamflow generation. Such undulating and variable topography is not captured by theoretical models. The development and extent of active areas and the dominant flow generation mechanisms estimated by the HMC method and the water balance results were almost identical for scenarios 1-4, indicating that knowledge of the exact hydraulic conductivity value of a given soil type is less importance than capturing the general vertical and longitudinal soil heterogeneity across the 450 catchment. The small differences simulated between scenarios 2 & 3 and 1 & 4 show that the models were more sensitive to variations on the hydraulic properties of the loam soil as the scenarios with identical responses (2 & 3 and 1 & 4) shared the same loam but different sand hydraulic properties. The loam's smaller hydraulic conductivity for scenarios 2 & 3 (0.0624 m 3 /d) limited infiltration, which translated to more OLF. At the same time, the higher water holding capacity in the loam areas might have resulted in slower subsurface flows to either exfiltrate to the surface or to contribute to the fluid transfer (FT) 455 through the subsurface boundary. In contrast, both tested values for the sand resulted in a low water holding capacity that allowed the incoming precipitation to drain past the root zone and move in the subsurface contributing to the FT. This is supported by the higher FT values shown for scenarios 1 & 4 (FigureFig. 6).
The greatest differences for both the water balance and HMC results among the tested scenarios were caused by differences in stream definition between scenarios 1-4 (no incised stream) and scenario 8 (incised stream). The effects of the incised stream 460 in scenario 8 resulted in a larger OLF flux, a smaller ET flux, and in an overall smaller extent of active areas (FigureFig. 5 and supplemental information). The roles of a good channel representation in ISSHMs is extensively discussed by Käser et al. (2014). While their discussion of channel representation revolves on the ability of ISSHMs to quantify GW-stream interactions (which was not a major component for this model), the hydrological principles are relevant and transferable to explain the importance of channel representation to capture streamflow generation processes. Without a defined channel, the connection 465 to the Pedler Creek floodplain was lost and the shallow sheet-flow and most of the precipitation infiltrated and stayed within the porous media. The importance of the connection of the floodplain to the channel is also discussed by Käser et al. (2014), as they argued that not only is the channel topography important but also its connection to the floodplain given that riverbank geometry is key for bank storage and overbank flooding. Although overbank flooding is not considered important for this study (flows in Pedler only rarely will experience overbank flooding), the stream-floodplain connectivity and bank storage 470 were key aspects under our model conceptualisation. Namely, the predicted dominant mechanisms relied upon either the saturation to build up along the riverine zone in the loamy areas which would lead to saturation excess overland flow; and we expected a perched groundwater developing on the sandy hillslopes which would, after intersecting the stream, contribute with interflow into the stream. While we observed these processes developed, they only occurred briefly as very shallow runoff.
Another important consideration for the channel representation in streamflow generation studies for non-perennial streams is 475 the relationship between flow and the wetted area. The larger the channel (both vertically and horizontally), the larger the area of exchange to the unsaturated zones during a flow event (Doble et al., 2012), which would be exacerbated under low flows (Käser et al., 2014). This is particularly significant when evaluating streamflow generation for non-perennial streams where high streambed infiltration and transmission losses are common (Gutierrez-Jurado et al., 2019;Levick et al., 2008;Shanafield & Cook, 2014;Snelder et al., 2013) and often prevent flows from even reaching the catchment outlet (Keppel & Renard, 1962;480 Aldridge, 1970). In this study, we observed that without a defined stream to 'channel' the water, the little overland flow that was simulated in scenarios 1-4 (no incised stream) spread over a larger area than in scenario 8 where the stream was incised (FigureFig. 5). The same was true for the patterns of increased saturation of the porous media across the catchment (FigureFig. 6). Results from the water balance reflected the effects of having both flows and porous media saturation spread over larger areas by exacerbating ET and decreasing the overall amount of overland flow for scenarios 1-4 (Supplemental Information). 485 This is consistent with the remarks by Käser et al. (2014) regarding the likely impacts to the water flow budget by the spatiotemporal aspects linked to channel representation due to spatial exchange patterns.
Finally, this study highlighted both the need for further studies examining streamflow generation processes in additional nonperennial catchments. For instance, our results underlined the importance of channel representation and future studies should investigate the effects of channel morphology in streamflow generation in non-perennial catchments. Moreover, while for 490 Pedler creek the GW-stream interactions were conceptualized to occur only near the catchment outlet (and likely only during certain 'wet' years), flow intermittence in many rivers can be attributed to a water table fluctuation relative to the stream channel elevation (Snelder et al., 2013). Future work on flow intermittence as a result of GW-stream interactions would be valuable. Similarly, and the inherent challenges associated with capturing unsaturated zone dynamics at catchment scale were underpinned in this study. Indeed, modelling this non-perennial river system confirmed the inherent difficulties of using 495 ISSHMs in medium size non-perennial river catchments and reiterated why so few studies have been done on this topic.
Extensive work was needed to set up this model, and the large computational time to run the simulations was a major constraint to both establishing initial conditions and exploring scenarios. For example, when draining the model, simulations running for over 10 days only progressed to day 100 of the simulation. Relaxing the mesh allowed us to develop reasonable initial conditions after testing over 37 scenarios. For the scenarios running for the full 4-year simulation, simulation convergence 500 consistently slowed down when the simulation encounter a precipitation input and particularly during prolonged precipitation events (i.e. consecutive precipitation inputs), which is expected given the highly non-linear equations for unsaturated flow in the unsaturated surface domain. Despite the conceptual advantages of using a fully physically-based model to explicitly capture all surface and groundwater processes, future studies may try to identify a suitable simplified surrogate model to speed up simulations and focus on specific areas where particular streamflow generation processes are thought to be dominant. 505

Conclusion
There are hundreds of similar non-perennial river systems in the semi-arid, coastal Mediterranean climate regions of South and Western Australia, California, South Africa, and around the Mediterranean itself (Davies et al., 1993;Tzoraki and Nikolaidis 2007;Skoulikidis et al., 2017). This study provides an initial step towards understanding non-perennial streamflow generation processes at the catchment scale and provides a template for using ISSHM's for process understanding in these 510 stream systems. The development of a conceptual model of the most important factors impacting flow generation processes within Pedler catchment presented a hypothesis that combined our understanding of field data with lessons learned from previous studies. This conceptualization was key to informed the model setup and to captured dynamics of streamflow generation in this non-perennial stream system. Difficulties in setting up and running the model reaffirmed the numerical difficulties experienced in large scale unsaturated models, such as accurately reproducing the topography and observed initial 515 conditions, was a challenge. and modelModel runtimes prohibited extensive exploration of multiple scenarios. In particular, the importance of preserving channel representation to model streamflow generation on non-perennial systems became apparent in the scenarios. Yet overall, the model results confirmed our conceptual understanding that soil type, unsaturated storage dynamics and topography are major controls for streamflow generation processes in non-perennial streams. The similarity in the results from scenarios comparing soil hydraulic properties across the literature range for each soil type showed 520 that exact knowledge of these values for a given soil type is not critical for identifying streamflow generation processes if the conceptual model is accurate and the vertical and longitudinal soil heterogeneity is captured. Given that soil properties are often highly heterogeneous within a catchment and rarely well-known, this result is will be important to future modelling studies.

Code/Data Availability 525
The data for this model were sourced through publicly available resources as cited in the methods. Commercially available AlgoMesh (HydroAlgorithmics) software and HydroGeoSphere (Aquanty) software was used to prepare and run model simulations. Post-processing routines and model files are currently in the process of being archived by Flinders University. The fraction of grass/bare soil is the main determining factor explaining the runoff response to different rainfall events