A channel transmission losses model for different dryland rivers

Channel transmission losses in drylands take place normally in extensive alluvial channels or streambeds underlain by fractured rocks. They can play an important role in streamflow rates, groundwater recharge, freshwater supply and channel-associated ecosystems. We aim to develop a process-oriented, semi-distributed channel transmission losses model, using process formulations which are suitable for data-scarce dryland environments and applicable to both hydraulically disconnected losing streams and hydraulically connected losing(/gaining) streams. This approach should be able to cover a large variation in climate and hydro-geologic controls, which are typically found in dryland regions of the Earth. Our model was first evaluated for a losing/gaining, hydraulically connected 30 km reach of the Middle Jaguaribe River (MJR), Cear á, Brazil, which drains a catchment area of 20 000 km 2. Secondly, we applied it to a small losing, hydraulically disconnected 1.5 km channel reach in the Walnut Gulch Experimental Watershed (WGEW), Arizona, USA. The model was able to predict reliably the streamflow volume and peak for both case studies without using any parameter calibration procedure. We have shown that the evaluation of the hypotheses on the dominant hydrological processes was fundamental for reducing structural model uncertainties and improving the streamflow prediction. For instance, in the case of the large river reach (MJR), it was shown that both lateral stream-aquifer water fluxes and groundwater flow in the underlying alluvium parallel to the river course are necessary to predict streamflow volume and channel transmission losses, the former process being more relevant than the latter. Regarding model uncertainty, it was shown that the approaches, which were applied for the unsaturated zone processes (highly nonlinear with elaborate numerical solutions), are much more sensitive to parameter variability than those approaches which were used for the saturated zone (mathematically simple water budgeting in aquifer columns, including backwater effects). In case of the MJR-application, we have seen that structural uncertainties due to the limited knowledge of the subsurface saturated system interactions (i.e. groundwater coupling with channel water; possible groundwater flow parallel to the river) were more relevant than those related to the subsurface parameter variability. In case of the WEGW application we have seen that the non-linearity involved in the unsaturated flow processes in disconnected dryland river systems (controlled by the unsaturated zone) generally contain far more model uncertainties than do connected systems controlled by the saturated flow. Therefore, the degree of aridity of a dryland river may be an indicator of potential model uncertainty and subsequent attainable predictability of the system.

a coupling of formulations which are more suitable for data-scarce dryland environments, applicable for both hydraulically disconnected losing streams and hydraulically connected losing(/gaining) streams. Hence, this approach should be able to cover a large variation in climate and hydro-geologic controls, which are typically found in dryland regions of the world. Traditionally, channel transmission losses models have 10 been developed for site specific conditions. Our model was firstly evaluated for a losing/gaining, hydraulically connected 30 km reach of the Jaguaribe River, Ceará, Brazil, which controls a catchment area of 20 000 km 2 . Secondly, we applied it to a small losing, hydraulically disconnected 1.5 km channel reach in the Walnut Gulch Experimental Watershed (WGEW), Arizona, USA. The model based on the perceptual hydrological 15 models of the reaches was able to predict reliably the stream flow for the both case studies. For the larger river reach, the evaluation of the hypotheses on the dominant hydrological processes was fundamental for reducing structural model uncertainties and improving the stream flow prediction, showing that both lateral stream-aquifer water fluxes and groundwater flow in the underlying alluvium parallel to the river course for disconnected streams, because in-channel ponding depth and gravitational terms are time-dependent (Freyberg et al., 1980). To overcome this difficulty, Freyberg (1983) proposed a numerical solution (trapezoidal quadrature) of the Green-and-Ampt equation for a uniform alluvium. His algorithm was initiated by the analytic solution to a nongravity approximation due to the singularity in infiltration rate at time equal to zero and 10 the inadequacy of the trapezoidal quadrature for rapid rate of change in infiltration rate at small time steps (Freyberg, 1983). Therefore, unsaturated flow through the alluvium, together with in-channel variable ponding depth, hampers a transmission losses model for disconnected losing streams. An extra difficulty might be the existence of an underlying stratified alluvium, which can often be found in dryland riverscapes (Parissopoulos 15 and Wheater, 1992;El-Hames and Richards, 1998).
Another approach for disconnected losing streams is the Smith-Parlange infiltration equation used in KINEROS2 model, which is based on an approximate solution of the basic equation of unsaturated flow (Smith et al., 1995;Semmens et al., 2008). The model requires basically three parameters (the integral capillary drive, the field effective 20 saturated hydraulic conductivity and soil water content) to describe the infiltration behavior, but the underlying soil profile can only be represented by two-layers with each layer allowed to have different infiltration parameters (Smith et al., 1995;Semmens et al., 2008).
Pressure-head-based Richards' equation enables us to model unsaturated flow 25 through the alluvium considering both in-channel variable ponding depth and stratified alluvium as done by El-Hames and Richards (1998). This might be the most physically comprehensive approach to model channel transmission losses for disconnected losing streams. However, its application can require a long processing time to simulate 8907 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | large-and meso-scale catchments (El-Hames and Richards, 1998) and large sets of alluvium data, which are usually not available, especially in dryland environments. Alternatively, some authors have used constant infiltration rates in the channels (Lange et al., 1999;Morin et al., 2009) neglecting both in-channel variable ponding depth and unsaturated flow. 5 In this paper, we aim to develop a semi-distributed channel transmission losses model, a coupling of formulations which are more suitable for data-scarce dryland environments, applicable for both hydraulically disconnected losing streams and hydraulically connected losing(/gaining) streams in dryland environments, considering a possible transition from the first behaviour to the former one and vice-versa, too. Hence, 10 this approach should be able to cover a large variation in climate and hydro-geologic controls, which are typically found in dryland regions of the world. Traditionally, channel transmission losses models have been developed for site specific conditions.
Our channel transmission losses model is firstly evaluated for a losing/gaining, hydraulically connected 30 km reach of the Jaguaribe River, Ceará, Brazil, which controls 15 a catchment area of 20 000 km 2 . Secondly, we apply it to a small losing, hydraulically disconnected 1.5 km channel reach in the Walnut Gulch Experimental Watershed (WGEW), Arizona, USA, which is well-known for its long-term database of semiarid hydrology and studies on channel transmission losses (e.g. Renard, 1970;Renard et al., 2008;Stone et al., 2008). 20 The model application to those channel reaches will be undertaken in order to evaluate the model capabilities in two very different dryland environments. Moreover, we will test hypotheses on the dominant hydrological processes for the larger river reach with a view to generating insights into reach functioning through comparisons of model performance (Savenije, 2009;Buytaert and Beven, 2011;McMillan et al., 2011;Clark 25 et al., 2011;Li et al., 2010). 8908 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Modelling of channel transmission losses
Conceptually, we considered the following processes, which have been shown experimentally to be the most influential on channel transmission losses (see discussion in introduction): 1. stream flow in natural rivers; 5 2. unsaturated seepage under in-channel variable ponding depth through a stratified alluvium; 3. vertical unsaturated subsurface water redistribution beneath the stream; 4. lateral (stream-)aquifer interaction, which includes the development of a groundwater mound and; 10 5. groundwater flow, parallel to the river course, in unconfined aquifers.
We also established possible in/outflow through the model boundaries of the channel transmission losses model, such as surface and subsurface hillslope runoff, evapotranspiration in the streambed and groundwater transpiration. These additional variables can be provided by process-oriented and (semi-)distributed hydrological mod-15 els. The model structure is composed of five components, which link spatially the sub-models of the above mentioned processes. Figure 1 shows an overview of the processes' spatial evolution: The channel transmission losses model couples six approaches: flood wave routing and in-basin stream flow distribution (Sect. 2.1), unsaturated stream infiltration model 20 (Sect. 2.2), vertical soil water redistribution model (Sect. 2.3), lateral (stream-)aquifer dynamics model (Sect. 2.4) and groundwater flow model (Sect. 2.5). These approaches interplay and proceed temporally as showed in Fig. 2.
The calculation begins with the flood wave routing without stream-aquifer interaction, i.e. we predict firstly stream flow and stream water stage excluding stream-aquifer in- 25 teraction flux. Then, we use these predicted "intermediate" values of stream flow and 8909 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | water stage to run the other sub-models (2, 3, 4 and 5), which estimate (a) the streamaquifer interaction flux and (b) the moisture in the underlying aquifer. Afterwards we apply the flood wave routing again, but now with the estimated stream-aquifer interaction flux, to predict finally the stream flow and water stage at the end of the time step. This kind of solution of stream flow and water stage is called a two-step procedure, 5 which was used e.g. by Mudd (2006) and Bronstert et al. (2005). As long as the stream-aquifer column is not saturated, the stream (1, 6 and 7 sub-models) and groundwater (4 and 5 sub-models) flows hydraulically disconnected, while channel transmission losses are dominated by the unsaturated zone beneath the stream (2 and 3 sub-models). Once the stream-aquifer column has been saturated, the 10 stream and groundwater turn into a hydraulically connected system, wherein channel transmission losses are driven by the saturated zone (4 and 5 sub-models), which can either discharge to (no losses) or recharge from the stream.
The following sub-sections from 2.1 to 2.5 describe the physical assumptions and the main mathematical formulations for the sub-models of our channel transmission 15 losses model. We detail the stream-aquifer interaction calculation in the last Sect. 2.6.

Flood wave routing
Normally the full Saint-Venant equations and its simplified diffusive-based form are applied to simulate stream flow in a natural drainage network, when the up-and downstream boundary conditions are available. However, in fact, most dryland streams have 20 no "fixed" downstream boundary conditions because many hydrographs end somewhere between initial stream flow and an assumed outlet. Moreover, poor monitoring can make the entire drainage network from the initial stream flow completely ungauged. Therefore, we proposed here an alternative flood wave routing to be applied to dryland rivers. 25 First, we use a form of conservation of mass equation (based on Fread, 1988): 8910 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | where t is the time (T), x is the length along the channel axis (L), Q is the stream discharge (L 3 T −1 ), A is the wetted cross-sectional area (L 2 ), s is the sinuosity coefficient (dimensionless), q is the lateral inflow per unit of length of channel (L 3 T −1 L −1 ) and I RA is the stream-aquifer interaction term per unit of length of channel (L 3 T −1 L −1 ), which can be stream infiltration (negative) or groundwater discharge (positive). 5 Applying the four-implicit numerical scheme (see Fread, 1993) to (1), we have where j and i are indexes of time and stream section, respectively. Equation (2) has two unknown variables: the stream discharge and the wetted cross-sectional area (related to the stream water stage) at the future time and at the next stream section: 10 Q(j + 1,i + 1) and A(j + 1,i + 1), respectively. Since, in natural streams, the channel morphology is a response of the stream hydrology, then the channel cross-section is a function of the past and upstream flood events. It means that all the information for the future and further stream discharge Q(j + 1,i + 1) is already "printed over" the wetted cross-sectional area A(j + 1,i + 1). 15 Taking this hypothesis into account, we get all the states for A(j + 1,i + 1) and substitute into Eq. (2) to find possible states for Q(j + 1,i + 1) according to the conservation of mass equation. Then, we average over the possible states of Q(j + 1,i + 1) and A(j + 1,i + 1), which obey the following simple physical rules Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Equation (3) can seem to be a physical filter of the states for Q(j + 1,i + 1) and A(j + 1,i + 1). If the next stream section is the last section of a sub-basin, the predicted stream discharge at this section, i.e. the catchment runoff from a sub-basin, is added as lateral inflow into a stream reach (in-basin stream flow distribution model).
The Courant-Friedrichs-Lewy (CFL) condition may be used as a condition for numer- where ν max is the maximum expected stream velocity (L T −1 ), ∆x min is the minimum stream reach and ∆t sim is the time step (T) for simulation. 10 We adapted here the modified Green-and-Ampt model proposed by Chu and Mariño (2005), because it might be a suitable compromise between computation time, data requirement and simplifying assumptions (e.g. constant infiltration rates). The alluvium beneath the stream ( Fig. 1) consists of N layers with hydraulic conductivities K N (L T depths of cumulative infiltration z N (L). When the wetting front is in a layer y at location z (Z y−1 < z ≤ Z y ), the governing equations are

Unsaturated stream infiltration model
8912 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | where f is the infiltration rate (L T −1 ), F is the cumulative infiltration (L), t is the time for the wetting front to arrive at location z and H 0 is the hydraulic head at surface (L), which was admittedly negligible in Chu and Mariño's formulation because their focus was on hillslope hydrology. Substituting Eq. (5) into Eq. (7) yields Separating Eq. (8), since the hydraulic head at surface is constant in a certain time step, we have Solving Eq. (9): which is similar to the equation of the travel time of the wetting front from Chu and Mariño (2005), but with the hydraulic head at surface H 0 . We use Eqs. (5) and (10) to estimate the actual infiltration and the location of the wetting front z. 15

8913
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Before applying the above procedure to the next time step, when we have a new hydraulic head at the surface, the initial soil moisture has to be updated according to the location z (Z y−1 < z ≤ Z y ) of the wetting front using In this model, the hydraulic head at the surface (the upper boundary condition) is the 5 average "intermediate" predicted values of stream water stage obtained from the solution of the flood wave routing (Sect. 2.1). The lower boundary condition is a layer, which can either represent fractured bedrocks (time independent) or be the soil layer immediately above the groundwater level (time dependent). Once the wetting-front achieves the lowest layer, a hydraulically connected streamlowest layer should now be considered and a groundwater mound is to be developed (Sect. 2.4). In contrast, the wetting-front flows vertically downward to the lowest layer (Sect. 2.3). For the first case, the infiltration rate tends to be constant and the capillary head zero as in Chu and Mariño (2005). Equation (5) can be rewritten as where Z N is the depth of the considered alluvium profile above the groundwater level and f z N is the infiltration rate for a hydraulically connected surface-boundary condition system. The infiltration rate from unsaturated to saturated regime can be formulated as Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | where ∆t is the time step. Note that the second element of the right term of Eq. (13) represents the first recharge to groundwater, if it exists, before the development of a groundwater mound.

Vertical soil water redistribution model
In most unsaturated zone studies, the fluid motion is assumed to obey the classical 5 Richards' equation (Hillel, 1980) and its 1-D soil moisture-based form is shown in the first two terms of Eq. (14), which is applicable in homogeneous media only and requires soil head-conductivity-moisture curves. We use here a simplification of the classical equation, which allows application in unsaturated heterogeneous media and needs less fitting parameters than the original form. First, we neglect the pressure head term ψ(θ) in Eq. (14), but we assume that percolation from one soil layer to the next layer below occurs if and only if the actual soil moisture exceeds soil moisture at field capacity θ FC . This assumption was also used in other hydrological models (e.g. Arnold and Williams, 1995;, leading to 15 We include in Eq. (14) the actual soil evaporation Eva for the upper soil layers and the actual evapotranspiration Eta for the soil layers in the root zone, in the case of existence of in-channel associated or riparian vegetation, which may be important for eco-hydrological studies and may allow insights into the relationship between channel transmission losses, in-alluvium temporal water storage and ecological water demand. 20 Furthermore we apply an explicit finite difference scheme to it: where k and j are indexes of depth and time, respectively. The percolation terms of Eq. (15) are solved as follows where K (θ) is a given approximation of K and θ relations. Note that if the lower layer 5 (k + 1) is groundwater, there is a recharge to groundwater before the development of a groundwater mound because of the vertical movement of the soil water.
A separate hydrological catchment model can provide the potential soil evaporation and the potential evapotranspiration. Then, we assumed that evapotranspiration and soil evaporation occur if and only if the actual soil moisture exceeds soil moisture at 10 permanent wilting point θ pwp and at hygroscopic water θ ha , respectively. The computation begins with percolation, followed by an updating of θ j k and then the transpiration calculation.

Lateral (stream-)aquifer dynamics model
We consider that each aquifer unit is formed by M columns, which have saturated and 15 unsaturated zones (Fig. 1). All these columns can be stratified such as that below the stream (Sects. 2.2 and 2.3). The lateral flow between the columns is considered saturated, consequently, we do not account for lateral unsaturated flow. Our aim is to predict in-column groundwater level (stream and groundwater levels for stream-aquifer columns), comparing the hydraulic heads between two column neighbours. During a time step, the calculation begins from the centre of the stream-aquifer column to the right (or the left) lateral boundary conditions (Fig. 1).
First, we calculate the hydraulic head of two column neighbours at the equilibrium (h e ), i.e.
where A is the column index (-), h is the in-column hydraulic head (L) and Cw is the column width (L). Then, assuming a subsurface water flow velocity similar to the order of magnitude of the lateral saturated hydraulic conductivity, we estimate the necessary time (d t e ) to reach that equilibrium head using whereK A+1 is an average lateral saturated hydraulic conductivity from the actual head to the equilibrium ones (L T −1 ). If d t e is equal to or smaller than the simulation time step ∆t sim , then the heads of the column neighbours reach the equilibrium head, otherwise

15
where h * A+1 is the new hydraulic head of column A + 1 due to the exchanges with the column A. Afterwards, the column A + 1 with the new hydraulic head h * A+1 will interact with its next neighbour A + 2.

8917
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Groundwater flow model
We use a simple water balance-based approach (similar to Niu et al., 2007) in order to simulate groundwater flow between aquifer units parallel to the river course (see Fig. 1) where S is the groundwater storage in the aquifer unit (L 3 ), Q Up,GW and Q La,GW are the 5 upstream and the lateral groundwater flow from other aquifer units (L 3 T −1 ), respectively, which are known from a previous time, Q V,Inf is the vertical channel transmission losses (L 3 T −1 ), which come from unsaturated seepage (Sect. 2.2) or unsaturated soil water redistribution (Sect. 2.3), Q Do,GW is the downstream groundwater flow (L 3 T −1 ), Q S is a sink term (L 3 T −1 ), which can be groundwater pumping and/or transpiration, 10 and Q V,DP is the vertical deep percolation (L 3 T −1 ), which is considered a constant (in)outflow. The downstream groundwater flow between aquifer units is estimated as follows where u is an index of aquifer unit,h is the average groundwater head of the aquifer 15 unit (L), W is the aquifer unit width (L), dx u is the aquifer unit length (L),K is the average aquifer unit saturated hydraulic conductivity (L T −1 ). Note that the downstream groundwater flow is compensated by a time factor, which is adopted similarly as was done in the previous section. The upstream boundary conditions are a constant (in)outflow. The modeller can define the downstream boundary conditions (a) as Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | no-flow or (b) assuming that the gradient of the downmost aquifer unit is equal to its closest upstream one.
After the estimation of aquifer water balance components, the difference between aquifer inflow and outflow is distributed for each column of the aquifer unit as follows where Q in(out),A is the in-column inflow or outflow from the aquifer water balance (L 3 T −1 ). If Q in(out),A is inflow, than the updating of in-column groundwater level due to the aquifer water balance is modelled by the actual groundwater level and b is the new groundwater level. On the other hand, if Q in(out),A is outflow, than where θf c is the soil moisture at field capacity (L 3 L −3 ). Moreover, if the in-column groundwater head overcomes the topographical maximum, the excess does not return 15 to the river network. Instead, it might firstly accumulate in depressions in the floodplain and then evaporate.

8919
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Stream-aquifer interaction calculation
The stream-aquifer interaction term per unit of length of channel I RA (L 3 T −1 L −1 ) (Sect. 2.1) can be estimated by wherehs andP s are the average stream water stage and wetted perimeter, respec- 5 tively, and f * is the potential infiltration determined in Sect. 2.2 as long as the streamaquifer column is not saturated. Once the stream-aquifer layer is saturated, then f * is calculated as follows where ∆h * is the increase or the decrease difference of the hydraulic head in the 10 stream-aquifer column determined using Eq. (20) in Sect. 2.4. If all the available stream water is to be infiltrated, then we apply no flood wave routing and set the predicted stream discharge and wetted cross-sectional area (related to stream water stage) equal to zero, in order to avoid numerical fluctuations when we use the stream-aquifer interaction term in the flood wave routing. 15

Case studies of the channel transmission losses model
We evaluated our channel transmission losses model for two stream reaches with different scales and dominated processes: a large reach of the Jaguaribe River, Ceará, Brazil and a much smaller one in the Walnut Gulch Experimental Watershed (WGEW), Arizona, USA. The data description of these sites and their parametrization are pro-20 vided in the following sub-sections.
We assessed different model structure strategies for the reach of the Jaguaribe River, in order to find out which model structure should better fit the channel transmission 8920 losses processes for that study site. This procedure is important if one intends to predict channel transmission losses in an ungauged stream reach, which preserves similar scale, database, climate and hydro-geologic controls with that reach. Therefore, since we intended to achieve the best long-term prediction of channel transmission losses, the best model structure was that which minimizes the root mean squared error 5 (RMSE) and the mean absolute error (MAE) of stream flow peak and event volume series. For the reach studied in the WGEW, we only calculated the RMSE and the MAE of stream flow peak and event volume series, because the uncertainty in the dominant processes involved in its channel transmission losses is relatively low in our point of view, which was based on previous publications and reports (e.g. Renard, 10 1970;Goodrich et al., 2004;Renard et al., 2008;Stone et al., 2008;Emmerich, 2008;Osterkamp, 2008). After that, we carried out a simple parameter sensitivity analysis, in order to guide the efforts on data acquisition and parameter calibration in future applications. We used for the parameter sensitivity analysis the following standard formulation where φ is the sensitivity coefficient, y is here a simulated variable, stream flow peak or event volume, and P is a model parameter. To carry out the sensitivity analysis, we selected the driest and the wettest stream flow events, whose upstream flow reached the lowest stream section.

Data and parametrization
We simulated a losing/gaining, hydraulically connected 30 km reach of the Jaguaribe River, Ceará, NE-Brazil, which controls a catchment area of 20 000 km 2 . The Jaguaribe River Basin's hydrology is determined by an annual cycle of rainy and dry seasons, 25

8921
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | which are driven mainly by the position of the Intertropical Convergence Zone and secondarily by cold fronts from the South Atlantic (Xavier, 2001;Werner and Gerstengarbe, 2003) and which produce about 1 m yr −1 of rainfall and 2.2 m yr −1 of potential evaporation (class A pan). The rainy season lasts up to six months (December-May) on average.

5
The simulated reach is dominated by unconfined aquifers (Fig. 3) belonging to an alluvium with a 20 m average depth and composed of layers of fine and coarse sand, gravel and clay (IBGE, 2003). According to Costa et al. (2011), on the one hand, during the dry and at the beginning of the rainy seasons, no river flow is expected for preevents and stream flow events have predominantly vertical infiltration into the alluvium. 10 On the other hand, at the middle and end of the rainy seasons, river flow sustained by base flow occurs before and after stream flow events and lateral infiltration into the alluvium plays a major role during events. Moreover, most channel transmission losses certainly infiltrated only through streambed and banks and not through the flood plains (Costa et al., 2011). 15 Measurements on the initial moisture of the aquifer-system were not available. However, since at the middle of the rainy seasons river flow is expected to be sustained by base flow, we may assume the groundwater level to be close to the river bed at the middle of the rainy seasons. Therefore, we applied the model from the middle of the rainy seasons, in which there was a big enough time shift between its beginning 20 and its middle, since at beginning of the rainy seasons the aquifer-system moisture is unknown and rather difficult to assume.
We assumed from Costa et al. (2011) that the actual inflow into the simulated reach is a sum of the actual stream flow measured at the N2 stream gauge, close to the confluence of the Cariús River into the Jaguaribe River, and the one-day-before stream 25 flow measured at the N1 stream gauge in the Jaguaribe River (see Fig. 3). The simulated output stream flow was compared to the stream flow measured at the N3 stream gauge in the Jaguaribe River.

8922
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | We used alluvial stratigraphy data, 15 boreholes and one electrical resistivity survey (Carneiro, 1993), and alluvium extension information from a hydrogeological map (Fig. 3) to derive the aquifer units (see Fig. 1). We used remote sensing-based data available from Costa et al. (2011) to delineate the channel length and the maximum channel width, whereas field observation provided the maximum channel depth. Then, 5 we derived stream cross-sectional areas by assuming a triangular channel crosssectional area. We did not account for infiltration into floodplains; since for our example it is not considered relevant for channel transmission losses (see discussion above in this section).
The simulated reach of the Jaguaribe River was spatially modelled (see Fig. 1) as 10 one basin system, which has one river with 4 reaches and 5 sections. Its aquifer system was formed by 4 units containing, respectively 7, 17, 13 and 21 (stream-)aquifer columns from up-to downstream. The typical up-to-down stratigraphy of an aquifer column was: sandy loam (topsoil), fine to coarse sand (1st alluvial layer), coarse gravel and very coarse sand (2nd alluvial layer) and silty clay (boundary condition), being 15 the last three for the stream-aquifer columns. Moreover, the soil layer interval was set at 0.2 m for all (stream-)aquifer columns. The texture of the aquifer was used to derive its soil physical properties, such as saturated hydraulic conductivity and porosity, obtained from laboratorial-experiments-based tables published in Rawls et al. (1993) and Dingman (2002). 20 The time step of the calculation (lead time), which gave the best numerical stability of flood wave routing and, consequently, used for this simulation, was 12 h. Since the original input time series were not sampled every 12 h, but only daily, we had to disaggregate them. 25 We Using those rainy seasons, we evaluated which model structure would provide the best simulation, i.e. the minimum of both RMSE and MAE of peak and event volume time series. Using the same parameter set and the spatial discretization, which were derived without calibration as shown in the last sub-section, we defined three possible model structures: (a) flood wave routing only, i.e. no aquifer system, (FW); 5 (b) flood wave routing with lateral (stream-)aquifer dynamics, but without groundwater flow parallel to the river course, (FW + LD); and (c) the same as (b) but now with parallel groundwater flow (FW + LD + GW). Figure 5a-c show the simulated and observed output stream flow series.

Model structure strategies evaluation
The FW-based model overestimated both the stream flow peak and the volume as 10 expected. The (FW + LD)-and (FW + LD + GW)-based models predicted similar peaks, but the (FW + LD)-based simulated hydrograph decreased more sharply during the recession flow than the (FW + LD + GW)-based one. The models' performance is shown in Table 1. The (FW + LD)-and (FW + LD + GW)-based models had comparable performance 15 and both were better than the FW-based. Because the (FW + LD + GW)-based model had the most similar behaviour to the observed hydrographs than the (FW + LD)-based one, we consider the (FW + LD + GW)-based model structure as the best suited for this study site. 20 Once the (FW + LD + GW)-based model structure presented the best simulation performance, we chose the following set of parameters in order to carry out the sensitivity analysis: (a) porosity/ soil moisture at field capacity/ initial soil moisture, which are related to the updating of the in-column groundwater level (Eq. 24); (b) lateral saturated hydraulic conductivity, which is related to lateral (stream-)aquifer dynamics (Eq. 19); 25 (c) "parallel" saturated hydraulic conductivity, which is related to groundwater flow parallel to the river course (Eq. 22).

Parameter sensitivity analysis
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Stream flow volume and maximum peak simulated by the (FW + LD + GW)-based model for the years 2005 and 2009 were used as reference variables (see Eq. 27), because those years were the driest and the wettest. Then, we multiplied a variable factor with the original values of the parameter sets (a), (b) and (c) and ran the (FW + LD + GW)-based model again, in order to estimate the sensitivity coefficients In general, high parameter values did not change the reference simulation, because large fluxes between model units are restricted by their hydraulic gradient (see Eqs. 19 and 22) and the reference simulation was already driven by the hydraulic gradient be-10 tween the model units. The sensitivity coefficient of porosity, soil moisture at field capacity and initial soil moisture was negligible even for small values (10 % of the original values). Also, between 100 % and 50 % of the original values of lateral and parallel saturated hydraulic conductivities, the sensitivity coefficient can be considered negligible. On the other hand, from 50 % of the original values of lateral and parallel saturated 15 hydraulic conductivities, the sensitivity coefficient could no longer be considered negligible. However, from 50 % to 10 % of the original values of those parameters, the sensitivity coefficient was between the range [−0.20; 0.20]. 20 We simulated here a losing, hydraulically disconnected 1.5 km channel reach in the Walnut Gulch Experimental Watershed (WGEW), Arizona, USA, from the flume FL008 (input flow) to the FL006 (output flow) (Fig. 8). Based on previous publications (e.g. Renard, 1970;Goodrich et al., 2004;Renard et al., 2008;Stone et al., 2008;Emmerich, 2008;Osterkamp, 2008), we assumed that stream flow infiltrates into an sandy alluvium 25 with enough depth such that it never becomes completely saturated during a stream flow event, because depth to groundwater within the WGEW ranges from ∼ 50 m at the  Spangler, 1969). Hydrological data and geo-information were made available at http://www.tucson.ars.ag.gov/dap/.

Data and parametrization
We selected hydrographs from stream flow events in which: 1. the input flow was only registered by the selected upstream flume (FL008); 5 2. the event volume, duration and peak flow at the selected upstream flume (FL008) were greater than at the downstream flume (FL006); 3. the soil moisture content of the underlying alluvium could be assumed close to the residual moisture content, i.e. at the beginning of the rainy season or after a long time between runoff events during the rainy season, since no soil moisture 10 data of the underlying alluvium were made available.
Maximum channel cross-sectional area and channel width were derived by stream channel morphology relationships provided by Miller et al. (2003). The stream reach under study presented a significant variation (+50 %) in its cross-sectional area from that its largest affluent reaches it (see Fig. 8). However, this variation was picked up by 15 Miller's relationships. The stream cross-sectional areas were then derived assuming a triangular channel cross-sectional area. We did not account for floodplains, because no data about them were available and also because we considered them to be of minor importance for transmission losses in that stream reach. Consequently, the floods had to be assumed as in-bank flows. 20 The simulated reach in the WGEW was spatially modelled (see Fig. 1) as one basin system, which has one stream with 3 reaches and 4 sections. Its aquifer system was formed by 3 units, each containing only one stream-aquifer column. The aquifer system was assumed to be uniformly sandy. Moreover, its soil layer interval was set 0.1 m for all stream-aquifer columns. The texture of the aquifer was used to derive its soil physical properties, such as saturated hydraulic conductivity and porosity, obtained from laboratorial-experiments-based tables published in Rawls et al. (1993) and Dingman (2002). 8926 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | The time step of calculation (lead time), which gave the best numerical stability of flood wave routing and which was consequently used for this simulation, was 2 min. Since the original input time series were not sampled for every 2 min, we had to resample them. 5 We selected 6 stream flow events which met the conditions described in the previous sub-section, in order to simulate the channel transmission losses from flume FL008 to flume FL006 (Fig. 8) using the parameters set and the spatial discretization derived without calibration, as shown in the last sub-section. Table 2 gives the comparison between the observed and simulated volume and peak flow of those events. 10 The volume of the events was clearly underestimated, its MAE being equal to 0.4 × 10 3 m 3 and its RMSE equal to 0.5 × 10 3 m 3 . The peak flow of the events was better predicted than its volume, where its error did not show a clear trend, its MAE being equal to 0.2 m 3 s −1 and its RMSE equal to 0.3 m 3 s −1 . We show the best and the worst predicted output hydrograph, which occurred on 29 August 1972 and on 2 August 15 1968, respectively, as follows.

Parameter sensitivity analysis
We selected the following set of parameters to carry out the sensitivity analysis: soil moisture at field capacity, pore size distribution index, porosity, wetting front suction and saturated hydraulic conductivity. Stream flow volume and maximum peak simulated for 20 the events on 28 July 1972 and 2 August 1968 were used as reference variables (see Eq. 27), because those were the driest and the wettest events. Then, we multiplied a variable factor with the original values of those parameters and ran the channel transmission losses model again, in order to estimate the sensitivity coefficients (Eq. 27) for stream flow volumes and maximum peaks. The sensitivity practically did not vary with changes on soil moisture at field capacity and pore size distribution index. In contrast, the sensitivity varied significantly with changes on porosity, wetting front suction and saturated hydraulic conductivity. Figure 11a-c shows the results of sensitivity analysis of those parameters for 28 July 1972 and Fig. 12a-c for 2 August 1968. 5 As expected, the sensitivity showed the largest values with changes in saturated hydraulic conductivity followed by wetting front suction and porosity. The higher those parameters are, the smaller is their sensitivity, i.e. the higher the infiltration from the stream into the alluvium. However, fluctuations in that behaviour could be found for peak flow in relation to wetting front suction, which might be related to numerical insta-10 bilities.
The sensitivity analysis for both applications showed that the sensitivity coefficient of saturated hydraulic conductivity and wetting front suction, which drive the unsaturated part of the channel transmission losses model, showed much higher values than that of the parameters which drive the saturated part of the model.

Discussion and conclusions
We developed in this study a new and fairly comprehensive channel transmission losses model, which was designed to simulate the surface-subsurface water fluxes in data-scarce dryland environments. Channel transmission losses modelling is indispensable for simulation of arid and semi-arid watersheds hydrology, as long as the 20 underlying aquifer system has not been fully saturated, as is expected to occur in river reaches at the beginning and in the middle of the rainy seasons (Renard, 1970;Costa et al., 2011). Moreover, after drought periods or during extensive groundwater pumping, (sub-)humid river reaches may resemble dryland rivers in relation to channel transmission losses processes. 25 Our model was able to predict reliably the stream flow for a large losing/gaining, hydraulically connected river and a small losing, hydraulically disconnected stream, because of its broad coverage of the dominant processes involved in channel transmission losses. Channel transmission losses models have not been developed and applied to different climate and hydro-geologic controls and scales as undertaken in this work.
The model structure strategies evaluation used for the application to the larger river 5 was shown to be pivotal for reducing structural model uncertainties and improving the stream flow prediction. Moreover, this evaluation provided a hydrological concept for ungauged river reaches, which preserves similar scale, database, climate and hydrogeologic controls with those of the studied reach. This concept means that both lateral (stream-)aquifer water fluxes and groundwater flow in the underlying alluvium parallel to the river course are necessary to predict stream flow and channel transmission losses, the former process being more relevant than the latter. The sensitivity analysis showed that the model results were bound by the parameters relating to the saturated part of the model (lateral stream-aquifer dynamics and groundwater flow parallel to the river course) because the water fluxes between the 15 model units were explicitly driven by the hydraulic gradient. In other words, even if the parameters can "potentially" produce large flow exchanges between model units in the saturated part, large flow exchanges do not happen because they are restricted by the actual hydraulic gradient between the model units. Moreover, the saturated-part-based parameters produced much smaller variation in the sensitivity coefficient than those 20 which drive the unsaturated part of the channel transmission losses model (unsaturated stream infiltration and vertical soil water redistribution). Thus, future efforts on data acquisition and parameter calibration should taken into account these results of sensitivity analysis.
The model presented not only a reliable prediction of stream flow volume, but stream 25 flow peak as well. The application for the small stream showed that the model performed even better simulating peak rather than volume. In this way, further research will be carried out to investigate the applicability of our model for runoff peak prediction in small and medium-sized dryland catchments. Moreover, a natural improvement 8929 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | will be the coupling of our model with a landscape hydrological model, such as the WASA model adjusted particularly for semiarid hydrology , for application and research in data-scarce dryland environments. The increase in data availability, principally from the subsurface, can allow a finer 5 spatial discretization of model units of the channel transmission losses model or even the development of more complex model structures, i.e. moving from the actual semi-distributed to a distributed hydrological modelling and/or adding new processes. This was, for example, the case for hydrological modelling in the Okavango Delta in Botswana, where surface-subsurface fluxes were simulated initially by concep-10 tual models and then by fully-distributed ones over the past decades (e.g. Bauer et al., 2006;Milzow et al., 2009), also due to increasing society demands on studies of manmade and climate change impacts on the Okavango Delta hydrology (e.g. Bauer et al., 2006;Milzow et al., 2009). However, even if more data are available, before increasing model complexity, it should first be evaluated that the actual simpler model matches 15 the objectives of the study (model parsimony).
Acknowledgements. The first author thanks the Brazilian National Council for Scientific and Technological Development (CNPq) for the PhD-scholarship. We thank David Goodrich for his comments on the literature review section and the section applicable to Walnut Gulch Research Watershed. We thank the Brazilian Geological Service (CPRM), the Brazilian Water