Implementation of the TOPKAPI model in South Africa : Initial results from the Liebenbergsvlei catchment

Flash floods and droughts are of major concern in Southern Africa. Hydrologists and engineers have to assist decision makers to address the issue of forecasting and monitoring extreme events. For these purposes, hydrological models are useful tools to: • Identify the dominant hydrological processes which influence the water balance and result in conditions of extreme water excess and/or deficit • Assist in generating both shortand long-term hydrological forecasts for use by water resource managers. In this study the physically-based and fully distributed hydrological TOPKAPI model (Liu and Todini, 2002),which has already been successfully applied in several countries in the world (Liu and Todini, 2002; Bartholomes and Todini, 2005; Liu et al., 2005; Martina et al., 2006), is applied in Africa for the first time. This paper contains the main theoretical and numerical components that have been integrated by the authors to model code and presents details of the application of the model in the Liebenbergsvlei catchment (4 625 km2) in South Africa. The physical basis of the equations, the fine-scale representation of the spatial catchment features, the parsimonious parameterisation linked to field/catchment information, the good computation time performance, the modularity of the processes, the ease of use and finally the good results obtained in modelling the river discharges of Liebenbergsvlei catchment, make the TOPKAPI model a promising tool for hydrological modelling of catchments in South Africa.


Introduction
Planning and management of water resources are of major concern in Southern Africa.Droughts and floods are extreme hydrological hazards that regularly face the population.Food security and public health are directly dependent on these hydrological phenomena (Kleinschmidt et al., 2001;Jury et al., 2002).In this context, environmental scientists and engineers are involved in two important challenges: • Providing short-term forecasts of the extreme hydrological events • Assessing how these events could possibly evolve in the context of ongoing global climate change.
For this purpose, the numerical modelling of the water cycle is essential in order to provide tools for operational hydrology in order to assist decision makers to develop suitable environmental strategies, and for improved understanding of the hydrological processes involved in the variability of the water cycle.
As part of a current Water Research Commission funded research project (K5/1683: Soil Moisture from Satellites: Daily Maps over RSA), as referred to in Pegram et al. (2006), the TOPKAPI model was adapted and coded using algorithms from the literature.The early development performed by Parak (2007) was reported in Pegram et al. (2007).One of the objectives within the current project is to use the distributed TOPKAPI catchment model to estimate the soil moisture from hydrological data and then to compare this estimate with remote-sensing estimates using two different types of satellite data.Thus the research-based model had to be further developed to run with practical hydrological applications in an operational manner.The initial results of the application of the model to compare simulated and remotely sensed soil-moisture estimates have been reported on in a companion paper (Vischel et al., 2008), which does not include much of the detail of model setup and calibration, issues of more interest to practising hydrologists.By contrast, the purpose of this paper is to present the development and application of the TOPKAPI model to a South African catchment, including the details of the model.There is naturally some material appearing in this paper which is common to Vischel et al. (2008), but the emphases are different.
In this paper the model TOPKAPI, initially developed by Liu and Todini (2002), is presented and applied to a catchment in South Africa (Liebenbergsvlei, 4 625 km 2 ).TOPKAPI meets most of the requirements expected by both operational and research hydrology: • It is a fully distributed model; the spatial range of the grid cell discretisation within which the model is valid up to 1 km (Martina, 2004).• It is physically-based, meaning that it explicitly represents the hydrological processes on the basis of the fluid mechanics and soil physics, while the input parameters can be directly obtained from existing spatial datasets.• There are relatively few (15) input parameters for a distributed model of which only three or four typically require calibration (Liu and Todini, 2002;Liu et al., 2005).
The TOPKAPI model has already been successfully implemented as a research and operational hydrological model in several catchments in the world (Italy, Spain, France, Ukraine, China) (e.g.see Liu and Todini, 2002;Bartholomes and Todini, 2005;Liu et al., 2005;Martina et al., 2006).The existing literature describing the TOPKAPI model was the point of departure of the work presented here in order to develop the code and devise an efficient hydrological modelling tool for South Africa.In this way, this work is also the first attempt to implement TOPKAPI in Africa.This paper presents the main aspects of the model: • Its principle and physical concepts • An example of an application of TOPKAPI to the modelling of a regional-size catchment in South Africa.
The TOPKAPI model background TOPKAPI is an acronym which stands for TOPographic Kinematic APproximation and Integration and is a physically-based distributed rainfall-runoff model.In the original version proposed by Liu and Todini (2002), TOPKAPI consists of 5 main modules comprising soil, overland, channel, evapotranspiration and snow modules.The first 3 modules take the form of non-linear reservoirs controlling the horizontal flows.The evapotranspiration module has been slightly modified in this application, compared to the original module presented in Liu and Todini (2002).The snow module component is ignored in the present study as the influence of snow-falls on runoff is negligible in the Liebenbergsvlei catchment.

Model assumptions
• The TOPKAPI model is based on 6 fundamental assumptions (Liu and Todini, 2002): Precipitation is constant in space and time over the integration domain (namely the single grid cell or pixel and the basic time interval, usually a few hours).This assumption simply means that the model is lumped at the grid scale.• All precipitation falling on the soil infiltrates, unless the soil is already saturated (Dunne, 1978) • The slope of the groundwater table coincides with the slope of the ground • Local transmissivity, like horizontal subsurface flow in a cell, depends on the integral of the total water content of the soil in the vertical plane • In the soil surface layer, the saturated hydraulic conductivity is constant with depth and, due to macro-porosity, is much larger than in deeper layers • During the transition phase, the variation of water content in time is constant in space.

Model equations
The equations controlling the level of the three main reservoirs that comprise a cell (soil, overland and channel reservoirs) are obtained by combining the physically-based continuity and mass equations under the approximation of the kinematic wave model.The well-known point-scale differential equations obtained are then analytically integrated in space to the finite dimension of a grid cell, which is taken to be the pixel of the digital elevation model (DEM) that describes the topography of the catchment.
For a detailed description of the theoretical framework used to formulate the TOPKAPI equations, the reader should refer to Liu and Todini (2002) or to Pegram et al. (2007).An overview of the relationship between the equations is given below.
The equation of mass continuity of each of the three reservoirs that compose a cell i can be written as a classical differential equation of continuity: (1) where: all the variables are observed at time t V i is the total volume stored in the reservoir is the rate of change of water storage is the total inflow rate to the reservoir is the total outflow rate from the reservoir The kinematic wave approach used to resolve the continuity and mass balance in TOPKAPI (by neglecting the acceleration terms in the St Venant energy equation) leads to a nonlinear relationship between and V i , transforming Eq. ( 1) into an ordinary differential equation (ODE) of the form: (2) where: is a combination of the forcing variables which are dependent on the reservoir type (soil, overland or channel), interconnecting flows between the elemental storage reservoirs within the cell and from upstream connected cells, including rainfall and evapotranspiration.For instance for the soil reservoir, is the sum of rainfall P (in mm) and flows coming from overland and soil reservoirs of the upper cells that have not been drained by the channel reservoir of the upper cells.Details are given in Fig. 1 in association with Table 1.b i is constant in time (it frequently varies spatially) and is a function of the geometrical and physical characteristics of the reservoir α (constant in space and time) and b i are dependent on each type of reservoir, as described below.

Soil reservoir
For the soil reservoir, the coefficient b i is expressed as:

with
(3) where: X is the lateral dimension of the grid-cell (in m) L i is the soil depth (in m) Ks i is the saturated hydraulic conductivity (in m/s) tan (b i ) is the tangent of the ground slope angle b i (with b i in degrees) θ s is the saturated soil moisture content (in cm 3 /cm 3 ) θ r is the residual soil moisture content (in cm 3 /cm 3 ) α s is a dimensionless pore-size distribution parameter (Brooks and Corey, 1964) originating from the infiltration equations (Liu and Todini, 2002).It usually takes on values between 2 and 4.

Overland reservoir
For the overland flow reservoir b i is expressed as: where: n oi is Manning's roughness coefficient for overland flow (in m -1/3 s -1 ) tan (β i ) is the tangent of the ground slope angle β i (with β i in degrees) β o is the well-known dimensionless power coefficient equal to 5/3 originating from Manning's equation Channel reservoir where: X c is the channel length (in m) (X c =X or X c = 2 X) W i is the width of the channel (in m) n ci is Manning's roughness coefficient for channel flow (in m -1/3 s -1 ) tan (βc i ) is the tangent of the channel slope (βc i ) (with (βc i ) in degrees) α c is as for the overland reservoir the dimensionless power coefficient equal to 5/3 again originating from Manning's equation.

Evapotranspiration
In this model coding, the evapotranspiration module has been slightly modified from the original version of Liu and Todini (2002).In the channel, the evaporation is extracted at the rate of the potential evaporation from a free water surface.On each cell i, the actual evapotranspiration Et a is computed as a proportion of the reference crop evapotranspiration Et r using, as a first approximation, a constant crop factor k c and the current saturation of the reservoir computed at each time-step t as the ratio between the effective and maximum soil water content (respectively V s (t) and V smax ), as shown in Eq. ( 6). (6)

Numerical framework of the TOPKAPI model
The automation of the TOPKAPI model has been achieved by coding the model in the Python programming language (Python, 2008).This free software was chosen mainly as it is suited for the management of large data arrays and because it can be extended and interfaced with many other existing languages (FORTRAN, C, C++, etc.).This feature is especially useful for resolving speed bottlenecks in computationally demanding routines.

TOPKAPI algorithm
For each cell, at each simulation time-step t, the inflow rate is computed, assumed to be a constant over the interval Δt, then the ODE equation is solved by numerical integration (see section 'ODE solution algorithm').
In Table 1 all the variables that are computed for each reservoir from the ODE finite difference solution showing the reservoir and cell connectivity are reported.Table 1 is associated with Fig. 1 which illustrates the fluxes and connections for a typical modelled cell.
At t+Δt, the evapotranspiration losses are then either extracted from the channel as well as from the overland flows if the cell is saturated or from the soil store alone if the demand is not satisfied by the overland storage.
with nb_cell_up the number of upper cells for the cell i.

Figure 1
Water balance in the TOPKAPI model (note that for clarity, the evapotranspiration losses are not represented on the figure)

ODE solution algorithm
The key point of the algorithm summarised above is the resolution of the ODE that controls the three reservoirs levels (Eq.( 2)).Particular attention was given to this point, since the ODE solution governs the ability of the model to: • Accurately compute each component of the water balance and satisfy local and global mass continuity • Provide a small enough computation time to use the model as a suitable tool for operational purposes.
Three methods have been carefully tested and compared to solve the ODE in the TOPKAPI model (with input term different from 0 in Eq. ( 2), since a straightforward analytical solution is known for this particular case): • The Runge-Kutta-Fehlberg (RKF) method (see e.g.Gerald and Wheatley, 1994) an extension of the Runge-Kutta algorithm, which yields the bonus of error estimates at each step is a well-known adaptive time-step algorithm that is widely recognised as one of the fastest and most numerically stable algorithms to solve ODEs in forward difference mode.• A quasi-analytical solution (QAS) based on a local linear approximation of the ODE.This method was proposed by Liu and Todini (2002) in order to reduce the computation time compared to the Runge-Kutta performances.• An analytical solution derived from integral tables (e.g. Gradshteyn and Ryzhik, 1994).
The last method was rapidly dismissed, since it performed poorly (in terms of computation time) when compared to the RKF method and it has the disadvantage of being restricted to a limited number of α values.The QAS method is faster than the RKF method but was shown to be unstable for some conditions in practice, since unexpected numerical divergences may occur, leading to an inability to satisfy the continuity of mass in the model.
Eventually a numerically stable and accurate hybrid scheme was devised to integrate the appropriate variations of the TOPKAPI equations.This hybrid method is based on the QAS procedure that is used by default, and switches to the RKF algorithm in cases where the mass continuity equations are not satisfied.This hybrid method was shown to be as stable and accurate as the RKF method when used on its own, but reduces the computation time by more than 50%.

Methodology to derive parameters for the TOPKAPI model from catchment information
The application of the TOPKAPI model to a catchment is relatively simple and requires data that are available from field measurements or remotely sensed observations.The methodology is based on the: • Determination of the characteristics of the geometrical entities of the catchment used in the TOP-KAPI model: the grid cell size, the cells composing the river network and how the cells are connected • Assignment of adequate parameter values in order to obtain realistic modelling of the catchment hydrology.
The first aspect is linked to the geomorphological features of the catchment which are extracted exclusively from the elevation data.The second aspect refers to the link between the data describing the elevation, soil and surface characteristics of the catchment (measured in situ or remotely sensed) and the physical parameters displayed in the TOPKAPI model equations.These two aspects are detailed below after presenting the features of the study catchment.

General description
The Liebenbergsvlei catchment (4 625 km 2 ) is a sub-catchment of the Vaal River catchment and is situated near Bethlehem in the Free State Province (Fig. 2).The climate is semi-arid, characterised by a mean annual rainfall ranging from 600 to 700 mm and a mean annual reference evaporation of between 1 400 and 1 500 mm (Midgley et al., 1994).

Validity of the TOPKAPI model assumptions for the Liebenbergvlei
Most of the assumptions listed in the Section 'Model assumptions' are independent of the catchment location except the second assumption which considers that the runoff production is exclusively due to saturation excess mechanisms.One might question the validity of this assumption for the Liebenbergsvlei catchment since infiltration-excess runoff processes (Hortonian processes) usually predominate in semi-arid environment.However, the predominance of saturation excess runoff production

335
seems to be realistic in the area for the major part of the season.This has been confirmed by some field experiments conducted at the hill slope scale in the region (Colin Everson, 2007).

Catchment data
The landscape is characterised by: • Hill slopes and steep relief in the southern part of the catchment which corresponds to the border with Lesotho and the Maluti Mountains • Grassland and cropland over the bulk of the catchment since farming is the main activity in the region.
These features are shown in the two first digital maps of Fig. 3a depicting a digital elevation model (DLSI, 1996) and land-cover/ use (GLCC, 1997).Information about soil types (SIRI, 1987) and soil texture (Midgley et al., 1994) is also available.

Hydrometeorological data
Hydrometeorological data on the catchment are available (Fig. 4, next page).A unique rain-gauge network consisting of 45 tipping bucket rain-gauges provided 5 min time-step ground rainfall measurement for the period 1993-2002.This network was not maintained after 2002.Since then, a network consisting of only 19 rain-gauges has remained.Ground measurements of rainfall are supported by the MRL5 C-Band weather radar operational since 1996.
Of the 9 flow gauges that can be identified on the catchment from the DWAF database, only three were shown to provide sufficiently long-term and good quality measurement to be used as reference flow data for testing the TOPKAPI model.Two flow gauges (labelled 1 and 2 in Fig. 4 -originally labelled C8H020 and C8H026 in the DWAF database) are available at the outlet of the catchment and further upstream, with unequal data availability and quality between 1993 and 2001.External flows arrive from Lesotho into the Ash River (tributary of the main river) via an inter-basin transfer, beginning in September 1997.These inter-basin transfer flows are recorded at a station located at the outlet of the transfer tunnel (labelled 3 in Fig. 4 -originally labelled C8H036 in the DWAF database).The quality of the flow data at Stations 1 and 2 has improved since 2002 but the recent flow data are not considered since, for the present study, it is not in the period covered by the dense rain-gauge network.

From catchment DEM to cell connection
The automation of the TOPKAPI model requires the definition of a numerical grid that divides the catchment space into squared cells that must be connected in order to model the transfer of flow (surface and subsurface) within the catchment.The Digital Elevation Model (DEM) is used as the base map to: • Define the grid and thus the spatial resolution of the model • Delineate the stream network.
These two steps can be achieved by using GIS software (here the ArcGIS TM software package has been used) and are thoroughly detailed in Parak (2007) and Pegram et al. (2007).
The spatial resolution of the model (parameter X in model equations, see the Section 'Model equations') is imposed by the grid characteristic of the DEM.In this application it was desirable to use the standard 1 km resolution usually provided by the freely available USGS DEM (USGS DEM, 2006).Thus, the initially fine resolution of 200m provided by the DEM available The 1 km DEM was then treated in order to delineate the streamflow directions.A common problem in DEMs is the occurrence of sinks.A sink, also referred to as a depression or pit, is a cell or area that is surrounded on all sides by higher elevation values.As such, a sink is an area of internal drainage and prevents the down-slope routing of water and, unless it is an actual case such as a lake or swamp, is an error.These errors often arise due to the sampling techniques used in processing a DEM or due to the rounding off of elevation values to integers (Mark, 1988).In order to create an accurate representation of the flow direction, it is best to use a DEM that is free of sinks, unless they are naturally occurring.In order to simplify the DEM treatment, Liu and Todini (2002) suggest the consideration of only 4 drainage directions (D4).However the limitation of the drainage to 4 directions can lead to an unrealistic representation of the relief variability.Indeed, the filling of the sinks using D4 in the Digital Elevation Model treatment results in a strong smoothing

337
of the relief variability because of the limitation of the drainage in only 4 directions (D4).For this reason, the TOPKAPI model was adapted to be compatible with 8 direction drainage (D8), which includes the 4 extra pixels beyond the diagonals.To fill a sink in D8, the methodology proposed by Mark (1988) was used exploiting the function 'sink' already available in the ArcGIS TM software.
From the depressionless DEM so created, the outflow direction of each cell is determined, i.e. the direction of the steepest outflow path from an active cell to the neighbouring downstream cells.A drainage direction code is assigned to each cell depending on the direction of maximum slope, which is calculated as the maximum difference in elevation divided by the horizontal distance from the centre of the active cell to the centres of the four surrounding cells.If the maximum slopes to several cells are the same, then the neighbourhood around the active cell is enlarged until the direction of steepest slope is found.The D8 flow direction raster map for the Liebenbergslvei catchment is shown on Fig. 5a.
The next step in the delineation of stream networks is to determine the number of upslope cells that contribute flow into each cell, i.e. the flow (in terms of contributing cells) accumulated in each cell.This is also achieved in ArcGIS™ using the standard tool 'Flow Accumulation'. Figure 5b shows the flow accumulation raster for the Liebenbergsvlei catchment.The colour palette indicates, for each cell, the number of upslope cells that feed it.The main trunk of the stream network is easily visible from this image, since it has the largest number of upstream cells.
The final step in delineating the stream network from a DEM is to assign a threshold value to the flow accumulation raster, for the minimum number of upslope cells that are required to initiate a channel in an active cell.The determination of this threshold value depends, according to Tarboton et al. (1991), on climate, slope and soil characteristics.They present procedures in order to 'rationally select the scale at which to extract channel networks' which correspond to networks obtained through more traditional methods, such as from topographic maps or fieldwork.Figure 5c shows a comparison between the stream network delineated from the DEM in the manner described above and a stream network digitised from a topographic map of the catchment (digitised by Midgley et al., 1994) at a spatial scale of 1:250 000.The comparison shows good correspondence between the two sources of networks, which is an important check.The threshold value chosen for the area over which a cell is considered to initiate a river channel (a parameter referred to as A threshold hereafter) was eventually fixed at 25 km 2 .In reality, the river network extension evolves within the seasons and with flow intensities but this value was considered to be representative of the average river network.Setting A threshold equal to 25 km 2 is in accordance with Todini's (1996) recommendation that the ratio between the number of cells containing channels and the total number of cells takes on a value ranging between 5% and 15% of the total catchment area; this threshold defined a total channel length of 522 km by 1 km pixel width, which is 11.3% of the total Liebenbergsvlei area of 4 625 km 2 .

From catchment data to physical model parameters
One of the most noticeable advantages of the TOPKAPI model, underlined by all the studies dealing with the model, is its physical basis that allows the user to link the model parameters with the catchment characteristics.The model parameters can be estimated a priori from the elevation data, soil and surface properties such as those available on the Liebenbergsvlei catchment but also available on the major part of South Africa.
A total of 15 parameters have to be assigned in the TOP-KAPI model.Among them, 11 are cell specific, meaning that they are potentially spatially variable (depending on the detail of the information available), and they mainly refer to physical characteristics (tan(β), tan(β c ), L, K s , θ r , θ s , n o , n c , α s , k c , W see Eqs. ( 2) to ( 6)).The 4 others are constant and refer to geometrical characteristics of the channel or grid cell (X, A threshold , the minimum width of channel W min , the maximum width of channel W max ).All the parameters and their associated values or range of values are reported in Table 2 (next page), as well as the references used to infer the parameter values.
As already noted in section 'From catchment DEM to cell connection', the constant parameters X and A threshold ,, have been assigned according to the DEM treatment to the respective values of 1 km and 25 km 2 .W min and W max were fixed at respectively 5 m and 35 m (estimations made from photographs taken at the flow stations).
For the cell-specific parameters, different degrees of treatment have been used to assign the values from the available data.Some of the parameters were directly provided by available datasets, some of them required references to the existing literature.Figure 2b illustrates how the catchment data maps were used to derive the maps of spatially variable model parameters.Maps of soil depths L and saturated soil moisture θ s were already available over the catchment in the data set of soil properties (SIRI, 1987).The slopes of the ground tan(β) (used for the soil and overland reservoirs) were computed from the DEM according to a neighbourhood function more representative of the mean slope within the cell and thus more representative of the transfers inside the cell (in and over the soil).This was achieved by using the function 'slope' in ArcGIS TM .The slopes used to transfer the flows in the channel drainage network tan(β c ) were computed from cell to cell in a down-stream direction using differences in altitude.Following Liu and Todini (2002), the poresize distribution parameter α s was uniformly set to the value 2.5.A sensitivity analysis (not presented here) showed that varying the value of α s in the realistic range of its values (between 2 and 4 according to Liu and Todini, 2002) had only a small influence on the simulations.As a first approximation, and because of the relatively homogeneous cropland/grassland land cover, the crop factor k c was assumed to be spatially uniform over the catchment and equal to 1.The width of channel in each channel cell was assigned by applying a linear relationship between the drained area and the channel width as proposed by Liu and Todini (2002): (7) where: A total is the area drained in the catchment A drained i is the area drained by a given cell i Particular attention was given to the 4 remaining parameters (K s , θ s , n o , n c ) that are not easily measurable in situ and have thus to be indirectly inferred by using literature references and tables.
The ordering method of Strahler (1957) was used to infer the values of the channel roughness Manning's coefficients n c .In Liu and Todini (2002), channel orders of 1, 2, 3 and 4 were assigned values of 0.045, 0.04, 0.035 and 0.035 for the Upper Reno catchment in Italy.These values were assumed to be suitable as starting values for the Liebenbergsvlei catchment for which no more precise information was available.The values of the overland roughness Manning's coefficient n o were derived from the land use/cover map (GLCC, 1997), using the widely used tables proposed by Chow et al. (1988).The residual soil moisture θ r and the hydraulic conductivity at saturation K s were derived from the soil texture map (Midgley et al., 1994) according to parameter tables for the Green-Ampt infiltration model (Maidment, 1993).

Model configuration
Once the model parameters are estimated a priori from the catchment data (as detailed in section 'Methodology to derive parameters for the TOPKAPI model from catchment information'), the application of the model requires some specific features defined in four main steps: • The definition of a modelling period • The choice of a simulation time-step • Pretreatment (aggregation/disaggregation) of the forcing data to conform with the spatial scale and time-step of the simulation • If necessary, an adjustment of some parameters in order to improve the model performance and account for the uncertainties in the a priori estimates.

Selected period
From the data set presented above, 2 seasons of 8 months each were selected during which the rainfall and flow data were both continuous and of good quality.The first season (Season 1) between November 1993 and June 1994 was used to adjust the parameters of the TOPKAPI model.The second season (Season 2) between November 1999 and June 2000 is used further as a model validation period.As it can be seen in Fig. 6, compared to the mean annual rainfall value 636 mm computed from the 1 st of November to the 31 st of October over the period 1993-2002, Season 1 is extracted from a relatively dry year (580 mm) and Season 2 from a relatively wet year (780 mm).

Model resolution
As explained in the Section 'From catchment DEM to cell connection', the spatial resolution of the model is imposed by the grid of the DEM.However, the simulation time is chosen by the user.Here, a 6 h time-step was chosen which is small enough to model the main discharge variations, since the catchment  (Midgley et al., 1994 ;Maidment,1993) fac Ks 60.

339
response time is estimated to be between 1 and 2 d.
It is worth noting that, at this time-step, only 20 min on a current PC are needed to run an 8 month season over the entire catchment (4 625 cells of 1 km).

Forcing variables
For the 2 seasons considered in this study, 6 h time-step rainfields were computed by Kriging the rain-gauge measurements at 1 km resolution using a climatological variogram with range of 30 km and a zero nugget (Wesson and Pegram, 2006).Because of the very good density of the network (10 km spacing) for the chosen period, the rainfall estimates from the weather radar were not used to create the rain-fields.
As no evapotranspiration data are available for the simulated periods on the catchment, the mean annual evapotranspiration over the region was used and disaggregated at a 6 h time-step, according to a mean seasonal signal determined from McKenzie and Craig (1999).

Calibration of the parameters
The TOPKAPI model is considered to be a physically based model, with all its parameters having physical meaning and which can be measured directly through fieldwork.But as in every physically-distributed model, the TOPKAPI model is subject to several uncertainties associated with the data on: • The information on topography, soil characteristics and land cover • The approximate methods and tables used to infer physical parameters from the data • The approximations introduced by the scale of the parameter representations.
For these reasons, Liu and Todini (2002) suggest that the calibration of parameters is still necessary but all studies dealing with the TOPKAPI model maintain that the calibration of the model is 'more of an adjustment' that can be achieved through simple trial-and-error methods.In the present application of the model on the Liebenbergsvlei catchment, the parameters estimated a priori from the catchment data were judged not accurate enough to be used as default parameters for the modelling experiment.
A calibration procedure was thus required.
In order to be convinced that only a small adjustment of the parameters is effectively sufficient for an efficient modelling of the Liebenbergsvlei catchment, an automated method of calibration was implemented.Inspired by the so-called Ordered Physicsbased Parameter Adjustment Method (OPPA) proposed by Vieux et al. ( 2004), the method consists of a dissociated calibration of the parameter responsible for the production of the runoff, from those responsible for the routing of runoff.According to a sensitivity analysis (not shown here, see also Liu et al., 2005) the most sensitive parameters controlling the runoff production are the soil depth L and the soil conductivity K s , while the Manning roughness of channel n c and overland n o are the main routing parameters.In the absence of any quantitative information, the initial soil moisture V s_initial , which was shown to have a strong influence on the simulations, was also calibrated.Ten values of mean catchment saturation between 1% and 90% were tested.Note that the a priori parameters are not calibrated cell by cell.This would lead to an extreme over-parameterisation of the model and to multiple and inconsistent combinations of parameter values.
The spatial pattern of the parameter maps is relevant information that was chosen to be conserved in the calibration procedure by using a multiplicative factor applied uniformly in space to the maps of the a priori parameters.For our application the 4 multiplicative factors to be applied were fac L (for the soil depth), fac Ks (for the hydraulic conductivity), fac no (for the overland roughness) and fac nc (for the channel roughness).
The triplet ( fac L , fac Ks , V s_initial ) was then adjusted in order to minimise the root mean square error (RMSE) objective function comparing modelled and observed discharge volumes aggregated at a monthly time-step.Once the optimum triplet was found, the pair ( fac no , fac nc ) was adjusted using the R 2 coefficient (as the objective measure) to match the timing of observed and modelled discharges at a 6 h time-step.
The calibration was first carried out using the flows estimated at station 2 (see Fig. 4).At this station, the drainage area is 3 563 km 2 , which effectively preserves the main soil heterogeneity of the entire catchment.

Calibrated period
Calibration was done on the flows of Station 2 in Season 1, shown in Fig. 7 (a).All other simulated flows in (b), (c) & (d) were verifications based on the one calibration.There is a good visual correspondence between observed and modelled hydrographs that is supported by a high value of the coefficient of efficiency (Nash and Suttcliff, 1970) of 0.78 and a good correlation illustrated in Figure 8a.
In Table 2 the values of the 4 calibrated multiplying factors are reported.It is worth noting that all the values of the parameters estimated a priori were quite appropriate except for the soil conductivity which has been increased by a factor of 60.The initial soil moisture was adjusted by calibration to a mean value of 40% over the catchment.
This strong increase of the soil conductivity made to satisfy the minimisation of the objective functions at first seemed puzzling, but it is actually rationally explainable and such calibration treatments of conductivity are fully reported in the literature (see e.g.Beven, 1997).First of all, it must be remembered that the tabled values used to assign the permeability values refer to the Green-Ampt infiltration model and thus are indicative of values used for vertical permeability of a column of soil experimented on in the laboratory.In this application, the concern is the horizontal transfer of water in soil that is explicitly  1994 1994-1995 1995-1996 1996-1997 1997-1998 1998-1999 1999-2000 2000-2001 2001-2002 Rainfall (mm) Avge = 636 mm Seasons (from the 1 st of November to the 31 st of October)

Figure 6
Mean (wet) seasonal rainfall computed between 1 st November and 31 st October from 1993 to 2002 over the Liebenbergsvlei catchment.These seasons match the duration of the calibration and validation periods in the study as indicated in Fig. 7.
modelled by TOPKAPI, whereas the vertical infiltration is implicitly lumped as an activity within the soil store.Lateral transfers mainly occur in preferential paths due to the presence of macropores, connected natural pipes and cracks.This means that the values of vertical permeability of the Green-Ampt model are not suitable for representing horizontal permeability and have to be increased to adapt to the representation of transfer rates that effectively occur laterally in the first metre of soil.

Evaluation of the calibration relevance on the entire catchment
As a verification of the relevance of the calibration procedure and its effect on other discharge time series, the calibrated model was applied to the entire catchment.For the same season (Season 1) the observed and modelled discharges at the outlet of the catchment (Station 1) are plotted in Fig. 7b.Globally, there is once again a good correspondence between observed and modelled flows, however at some points, the observed data seem to be unreliable since some peaks recorded at Station 2 do not appear as they should at the outlet and the recession shape of the main peak discharge seems somewhat unrealistic.This is clearly indicated on Figure 8b where the shape of the scatter plot of simulated versus observed discharge values presents an unexpected accumulation of observed discharge around the value 10 m 3 •s -1 .

Validation of the calibration procedure using data from another season
In order to validate the model calibration on more reliable data, the model was then applied to a different season (Season 2) some years later.During this season, the discharges are influenced by the inter-basin transfer flows now arriving from Lesotho.In order to reliably compare the modelled and observed discharges, the external flows observed at Station 3 had to be taken into account in the simulations.The modular structure of the TOPKAPI model allows the user to easily integrate additional fluxes that could come from new modelled processes, external flows or reservoirs.This is facilitated by the fact that a cell can be considered as an independent entity controlled by its own characteristics and its input.Here the external tunnel flows were simply added to the 'natural' input term of the cell positioned at the location the closest to Station 3. Again in the absence of any information about the initial soil moisture, the value of 40% calibrated for Season 1 was assumed to be applicable for Season 2. Hydrographs are plotted in Figs.7c and d, and scatter plots of observed versus simulated discharge are plotted in Figs.8c and d.Again, good simulations of the hydrographs were obtained even if the main peak discharges are underestimated.One can also note that the timing of the flows is remarkably good, especially at the beginning of the season, when in the absence of rainfall, the flows are mainly explained by the external flows that are routed from the upper part of the catchment; these appear pulsed due to hydropower generation.

Conclusion
TOPKAPI is a physically-based distributed hydrological model that couples the continuity equation with the soil and surface dynamic equations, resulting in a set of three kinematic wave equations applied to soil, overland and channel non-linear reservoirs within each cell, to which an evapotranspiration module is added.Its implementation is mainly based on elevation data (provided by a Digital Elevation Model) and also requires information about catchment surface properties and land use.
The present study is the first application of the TOPKAPI model on an African catchment, while the performance of the model has already been demonstrated in several catchments throughout the world.Applied on the Liebenbergsvlei catchment (4 625km 2 , South Africa), the model showed good ability in modelling the river discharges at a small (6 h) time-step with a limited adjustment of the parameters, and low computation times.341 Through the exercise of numerical programming of the TOPKAPI model from the existing literature, its implementation and the present application on the Liebenbergsvlei, significant expertise and understanding of the model have been gained.Several improvements to the original formulation have been made, for example an improved evapotranspiration scheme and a competitive algorithm for solving the differential equations, as outlined in this paper.From this experience, the following points indicate that TOPKAPI is a very promising hydrological model in the South African context: • Physical realism.A comprehensive and robust physical framework has been developed by Liu and Todini (2002) to give an explicit representation of the most significant hydrologic processes explaining the hydrological response of catchments.This response is the runoff production due to the subsurface flows and contributing saturated areas, surface transfers on hill slopes and river channels and evapotranspiration losses.• Spatial representation of the catchment properties.The spatial discretisation of the catchment at fine resolution (1 km) displayed by the DEM gives a distributed knowledge of the water fluxes that occur throughout the catchment, both surface and subsurface.This fine resolution is also coherent with the scales at which the processes effectively occur.

• Sensible parameterisation based on observed values from
the field.The parsimonious parameterisation of the TOP-KAPI models limits the number of parameters to 15, many of which are preset by the data available.Because of its physical basis, the parameters can be assigned according to field data obtained from in-situ measurements or remote sensed data, the parameters have a physical meaning attached to them, only a limited number of the parameters have to be adjusted due to a priori uncertainty, the parameters are relatively scale independent up to the limit of 1 km resolution as shown by Martina (2004).• Accuracy and computation speed.In this application of TOPKAPI a combination of a quasi-analytical solution (Liu and Todini, 2002) with a numerical integration procedure based on the Runge-Kutta-Fehlberg method was used.This hybrid scheme, which differs from the scheme initially suggested by Liu and Todini (2002), was chosen after a comparison with other methods, by carefully analysing the performances in terms of computation speed, numerical stability and accuracy of the simulations.• Modularity.The structure of the TOPKAPI model is based on the connection of independent cell entities.This allows the easy addition of reservoirs or dams, external flows (as was done for the Liebenbergsvlei case study) or supplementary hydrological processes that might modulate the hydrological response of the catchment.• Ease of use.Because of its simple and highly comprehensive reservoir structure, its parsimonious parameterisation and the direct physical link between the catchment data and the model parameters, the TOPKAPI model is very easy to use.
Combining the representation of the rapid flows associated with flash flood events and the fast computation time, TOPKAPI is a suitable tool for operational hydrology.It is also a useful tool for scientific and research hydrology.The TOPKAPI model has already been used as a research modelling tool in the project SHARE (a TIGER Innovator contract with the European Space Agency) to evaluate the soil moisture products derived from the satellites ERS 1-2.Beside its relevance of modelling the flow discharges of the Liebenbergsvlei catchment, the model also shows very good performances in modelling spatial soil moisture within the catchment as shown by the strong correspondence found between remotely sensed and modelled soil moisture (Vischel et al., 2008).Obviously one can discuss the physical basis of the model regarding the high calibration factor applied to the hydraulic conductivity K s during the calibration process (multiplying factor of 60).As already discussed in the Section 'Calibrated period', the problem mainly comes from the fact the values of K s estimated a priori were derived from Green-Ampt infiltration model tables that are associated with the local scale of a column of soil and for vertical infiltration fluxes.However these values underestimate flow due to the effect of macropores and preferential paths that play a crucial role in catchment response (Beven and Germann, 1982;McIntosh et al., 1999).The difficulty of estimating hydraulic conductivity of soils that control the subsurface lateral transfers at hillsope scale remains an issue.Of course it raises the question of the applicability of the model on

342
ungaged catchments, but recent studies propose several methods to regionalise the model parameters even if they have been calibrated (see e.g.Merz and Blöschl, 2004).
As a future extension of this work, improvements will be made in assessing the evapotranspiration rates which force the hydrological model.In the absence of data over the Liebenbergsvlei, the potential evapotranspiration has been coarsely estimated from mean annual values disaggregated into 6 hourly values by simply applying a seasonal signal extracted from the literature.This aspect will be improved by using the evapotranspiration estimates provided by Sinclair et al. (2007) who recently developed a method to compute reference evapotranspiration, on an hourly and daily basis over South Africa by combining remotely sensed and weather model data.
Then, by exploiting the easy modularity of TOPKAPI, it is planned to implement supplementary processes in the model in order to enlarge its applicability to a wide variety of hydrological systems.The response of some hydrological systems is more strongly controlled by runoff production mechanisms associated with rainfall rates exceeding the infiltration capacity of soil than by rain falling on already saturated areas.For these systems which mainly exist in very dry semi-arid to arid areas, the infiltration excess mechanism, also called the Hortonian mechanism, needs to be explicitly represented in the model.
Finally, the infiltration of water from subsurface to deep aquifer, whose contribution to rapid surface flows was assumed to be negligible in the present study, needs to be taken into account in order to have a complete insight of the water balance at regional scale.
Figure 2The Liebenbergsvlei catchment Figure 3 (a) Catchment data maps.(b) From catchment data to TOPKAPI model parameters

Figure 5
Figure 5From DEM to stream delineation.The red line superimposed on Figure5 c, is the river network digitalised from a topographic map at a spatial scale of 1:250 000(Midgley et al., 1994).

Figure 4
Figure 4Hydrometeorological data availability on the Liebenbergsvlei catchment

Figure 7
Figure 7 Comparison between modelled and observed discharges on the Liebenbergsvlei catchment for calibration and validation procedures.Calibration was done on the flows of Station 2 in Season 1, shown in (a).All other simulated flows in (b), (c) & (d) were verifications based on the one calibration.

Figure 8
Figure 8 Scatter plot of observed vs. simulated discharges on the Liebenbergsvlei catchment for calibration and validation procedures.Calibration was done on the flows of Station 2 in Season 1, shown in (a).All other simulated flows in (b), (c) & (d) were verifications based on the one calibration.

TABLE 2 Values of the TOPKAPI model parameters estimated a priori from data and literature, and values of multiplying factors obtained from the calibration procedure
3 ) Available on website http://www.wrc.org.zaISSN 0378-4738 = Water SA Vol.34 No. 3 July 2008 ISSN 1816-7950 = Water SA (on-line)