{'fraction_non_alphanumeric': 0.025576337538060025, 'fraction_numerical': 0.011193272437291577, 'mean_word_length': 5.2014026254270815, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-linear Routing Scheme at Grid Cell Level for Large Scale Hydrologic Models: A Review

New tools and concepts in the form of mathematical models, remote sensing and Geographic Information System (GIS), communication and telemetering have been developed for the complex hydrologic systems that permit a different analysis of processes and allow watershed to be considered as an integrated planning and management unit. Hydrological characteristics can be generated through spatial analysis, and ready for input into a distributed hydrologic models to define adequately the hydrological response of a watershed that can be related back to the specific environmental, climatic and geomorphic conditions. In the present paper, some recent development in hydrologic modeling will be reviewed with recognition of the role of horizontal routing scheme in large scale hydrologic modeling. Among others these developments indicated the needs of alternative horizontal routing models at grid scale level that can be coupled to land surface parameterization schemes that presently still employed the linear routing model. Non-linear routing scheme will be presented and discussed in this paper as possible extension.


INTRODUCTION
The hydrological response of river basins is controlled by a complex function of environmental (Idris and Mahrup, 2017;Welde and Gebremariam, 2017), climatic (Cristiano et al., 2017;Meaurio et al., 2017), and geomorphic processes (Dhali and Biswas, 2017;García-Comendador et al., 2017;Sahoo and Jain, 2018), especially in the tropical environment (Shukla and Gedam, 2019;Zema et al., 2018). However, for such knowledge to have wide applications, the hydrological response function needs to be understood in terms of the processes driving it, such as climate (Joo et al., 2017), vegetation (Sun et al., 2017), drainage network (Sofia and Tarolli, 2017), soils and land uses (Zhang et al., 2016); and recognize the patterns governing nature of the system. The introduction of system concept into the complex hydrological processes of watershed system in the past decades has made it possible to develop hydrologic system models that incorporated all the driving factors in an integrating manner, partly realizing the dreams of theoretical hydrologists for complete predictions based on observed parameters. Models proved to be practical for daily life to sustain ecosystems. Blöschl and Zehe (2005) commented the importance of non linearitis of hydrologic systems causing large discrepancies due to small uncertainties in initial and boundary conditions. With the system concepts integrating these complex processes have come into an understanding of river basins or watersheds as dynamic natural systems consisting of interconnected components and processes that form the water cycle (Loucks and van Beek, 2017;Su et al., 2020). The endless recirculation of water in the atmospherehydrosphere-lithosphere is known as the hydrologic cycle. The cycle can be studied according to a particular scale, such as global, regional or basin scale. The basin scale cycle that is known as run off cycle can be considered as an open system of continuous flow of water in the land phase from precipitation to interception, infiltration, runoff, and flows into streams, lakes, reservoirs, soil moisture, groundwater and out of basin transfer, to the return of water vapor through evaporation and transpiration (Bhardwaj, 2019;Stephens et al., 2020). Within a basin, the dynamics of the hydrological processes are governed partially by the temporal and spatial characteristics of inputs and outputs and the land use/land cover conditions. Anthropogenic influences were obvious in these land conversions from natural vegetations to developed areas that the increase of surface runoff and decrease of base flows are associated with the land use change (Makhtoumi et al., 2020). New tools and concepts have been developed that permit a different analysis of processes (mathematical models (van Kempen et al., 2020), remote sensing and Geographic Information System (GIS) (Quinn et al., 2019), communication and telemetering ) and allow watershed to be considered as a management unit (Setiobudi and Sembiring, 2009). A new era of spatial science is also obvious with development of these new tools and concepts. Hydrological characteristics can be generated through spatial analysis, and ready for input into a distributed hydrologic models to define adequately the hydrological response of a watershed (Singh, 2018). This hydrological response then needs to be related back to the specific environmental, climatic and geomorphic conditions.
Predicting the hydrological response of a river basin can be accomplished using regionalization techniques, assuming that new watersheds behave similarly to gauged watersheds (Muharsyah et al., 2020;Pagliero et al., 2019). Because there are many different approaches to regionalization, and because these approaches are often region specific, collaboration between hydrology, GIS and remote sensing, and mathematical techniques for environmental modeling is especially important. In the present paper, some recent development in hydrologic modeling will be reviewed with recognition of large scale hydrologic modeling. Among others these developments indicated the needs of a large scale horizontal routing model at grid scale to be coupled to land surface parameterization schemes that presently (Kauffeldt et al., 2016;Lohmann et al., 1996;Thober et al., 2019) still employed the linear routing model. Extension to non-linear routing scheme will be presented and discussed in this paper.

RECENT DEVELOPMENTS IN HYDROLOGICAL MODELING
Probably the reason behind recent development in large scale hydrologic modeling is the relative political stability in the past two decades almost all over the world, so that large scale international cooperations and assistance flourished. One of such international cooperation is the Mekong Project with Mekong Basin studies (Lu et al., 2020;Pokhrel et al., 2018;Try et al., 2018;Wang et al., 2016). Others are the Large Basin Experiments of the Amazon (LBA) Project (Rodriguez and Tomasella, 2016;Wongchuig Correa et al., 2017) and the Baltic Sea Experiment (Baltex) Project (Gröger et al., 2021;She et al., 2020). It is quite obvious to consider the significance of these large scale experiments and investments that made recent developments in large-scale hydrologic modeling possible. Some questions raised (Herath and Dutta, 2000) were: How do we model large basins, how can we incorporate land use changes as a dynamic processes linked to water resources, what techniques can we use to derive physical watershed characteristics, what different types of models are required for modeling diverse hydrological processes, what are the difficulties in studying large basins? The significance of large basins in any region is also obvious of their contribution to provide water resources the the region indicating them as river basin production systems (Abou Rafee et al., 2019;Gawne et al., 2018;Liersch et al., 2019). The case of Java island, the three largest river basins out of hundreds of other major basins occupied areas approximately 40 percents of the island. A single Memberamo river in West Papua province occupied about 30 percent of the province, and more so in Kalimantan island which is dominated by large river basins. It is good to note also the comment given in Flood Study Report (Beven, 1986) that: `larger catchments would appear to behave more regularly (compared to smaller basins)', indicating a better understanding on larger river hydrology should be more attainable.
Regular grid or raster digital elevation models (DEMs) have become the basis for recent approaches to process modeling of the earth's surface systems, especially hydrological modeling incorporating spatially variable parameters (Lutz, 2018). Raster DEMs can be generated directly from stereophoto maps when these are available, or more recently from satellite imagery, but there remains a significant role for the interpolation of DEMs from scattered point elevation data, perhaps accompanied by streamline data, particularly when the point data include surface specific points such as peaks, pits, saddles and selected points on stream lines and ridge lines (Habib et al., 2020;Xiao-Ping et al., 2016). Hutchinson, (1989) described a morphological approach to the interpolation of digital terrain data which attempts to take into account the special nature of terrain surfaces, and the surface specific points that can be used to sample terrain, as well as potential hydrological applications of the interpolated elevation grid. It has given rise to a procedure which can efficiently calculate raster DEMs with sensible drainage and streamline data. The principal innovation of the procedure is a drainage enforcement algorithm which automatically removes spurious sinks or pits from the fitted grid, in recognition of the observation that sinks are rare in nature. Contribution of this kind of work is to be recognized in recent developments of hydrologic modeling.
With the increased recognition of the importance of feedback mechanisms between land surface processes and climate, there have been sustained efforts to develop more realistic land surface representations to be couple with the general circulation models (GCMs) (Boer et al., 2007;Lestari and Dasanto, 2019;Müller et al., 2021). The current generation of land surface models used in GCMs view the soil column as the fundamental hydrologic unit that effectively ignores the role of topography plays in the development of soil moisture heterogeneity and the subsequent impacts of this soil moisture heterogeneity on watershed evapotranspiration and the partitioning of surface fluxes (Rakovec et al., 2016). Stieglitz et al., (1997) presented an approach to land surface modeling that allows us to view the fundamental hydrologic unit as the watershed rather than the soil column, employing the role of topography in the timing of discharge and the partitioning of discharge into surface runoff and baseflow. The analytic form of TOPography-based hydrological MODEL (TOPMODEL) equations are incorporated into the soil column in a consistent fashion. Soil moisture heterogeneity represented by the saturated low-lands subsequently impacts the partitioning of surface fluxes, including evapotranspiration and runoff (Stoy et al., 2019). The approach was claimed computationally efficient that allow for an improved simulation of hydrologic cycle, and is easily coupled into the existing framework of the current generation of single column land surface models. And because this approach uses the statistics of the topography rather than the details of the topography, it is compatible with the large spatial scale of today's regional and global climate models. This work was further evaluated (Warrach et al., 2002) to incorporate and compare two methods of sub-grid variability in soil moisture and runoff production into Surface Vegetation-Atmosphere Transfer (SVAT) models: (1) the variable infiltration capacity (VIC) model; and (2) a modified TOPMODEL approach. Because neither approach needs to explicitly track surface or subsurface flow within a catchment, they represent computationally efficient ways to represent hydrologic processes within the context of regional and global modeling. The study shows that during low flow periods the baseflow simulation is superior when using the TOPMODEL based runoff formulation, especially during the accession of the hydrograph in autumn. This is due to the fact that it accounts (in a quasi-statistical way) for the water table dynamics. A main drawback of the modified VIC approach, especially for regional and global application, is that with five free parameters, significantly more model calibration is required (Dang et al., 2020). TOPMODEL, on the other hand, only requires the determination of one free parameter and extensive pre-processing of topographic data (Xue et al., 2018), therefore gaining more popularity.

TOPMODEL Concept
The central assumption of the TOPMODEL formulation is that the saturated sub-surface store has an exponential flow law, which provides a satisfactory fit to catchment storm response in a range of circumstances and also provides spatially uniform flow runoff from spatially uniform inputs to the saturated layer (Kirkby, 1986), as may be seen below. The assumed exponential flow law for saturated flow at any site may be written in the form: where qo is the saturated soil discharge on unit hydraulic gradient, and the scaling soil parameter, m is assumed to be spatially uniform along the strip. Expanding the first term in the equation and substituting for S, it may be seen that j is the only timedependent term in the partial differentiation with respect to time, so that, with the exponential store assumption, the continuity equation may be written as: Flow-line strips are defined as following lines of greatest slope, orthogonal to elevation contours. Distance along the strip, x, is measured from the divide in a horizontal direction following the local flow-line direction. The width of the strip is defined by its width, w, at each point along the flow-line. Drainage area per unit strip width is described by the geometry relationship: The geometry of the strip may be described by two other geometrical identities which are derived below. From the equation above, Along a flow strip, the continuity or storage equation for conservation of water mass may be written as: Where: q is the local saturated discharge per unit width, S is the local saturated water storage (defined to be zero at the ground surface and negative for deficits below saturation), t is elapsed time, and, i is the local rate of percolation to the saturated zone.
Substituting the above geometrical relationship and rearranging terms and dividing through by w, then the continuity equation can be written as: This expression provides a basis for routing saturated flow down the length of the hillslope strip, and water balance accounting at each time step allows the evolution of the mean overall storage to be estimated. This can then be redistributed over the catchment at each time step on the basis of spatial distribution of an appropriate index that is derived from topographic and soil characteristics of each point in the catchment. The local values of storage can then be used to identify the surface contributing areas predicted at each time step: the higher the index value, the wetter the point and the more frequently a point will be saturated to a given level, relative to other points in the catchment (Ajami and Sharma, 2018). The use of such an index, which can be considered as an index of hydrological similarity, is important to the simplicity of the TOPMODEL concept because it is not necessary to carry out calculation for every point of space, since every point with the same index value will have the same predicted response given the same local inputs. The TOPMODEL concept can thus be considered as one form of disaggregation approach to modeling the variability of hydrological responses at the subcatchment level. The objective is to find an appropriate catchment scale parameterization for a particular set of circumstances, given the available understanding and measurements. More detailed descriptions of the concept can be obtained from readily available literatures in hydrologic modeling (Kirkby, 1986). Artan et al., 2000 andDobarco et al., (2017) investigated the appropriate spatial scale for a distributed energy balance model by (a) determining the scale of variability associated with the remotely sensed and GIS-generated model input data; and (b) examining the effects of input data spatial aggregation on model response. In order to determine the optimum grid size when partitioning the watershed to model the hydrologic processes in a distributed manner, the guiding criteria should be: (a) minimize the computation time by reducing the number of grid cells; while (b) at the same time maximizing the variation between grids in order to capture the significant patterns in the watershed; and (c) keeping the nonlinear effects of subgrid heterogeneity on the model output to a minimum.

Modeling Grid Size
The question of what a valid grid size would be in any distributed hydrologic models had been an important issue for some time as to relate to the appropriateness of the use of the models to certain field conditions. Earlier distributed watershed models usually would limit the size of the grids to several hectares only so that the above criteria be satisfied, in the same way as the early rainfall-runoff relationships of lumped model types were assumed . The famous rational method or the unit hydrograph analysis approach has an implicit assumption to be valid only for small catchment size, may be up to the order of 100 sq.km and not that small as to only few hectares (Jainet, 2018). And actually this is close to the spatial resolutions of present general cli-mate models that need compatible land surface parameterization schemes. Therefore it is to be recognized now the idea of scale independent in hydrologic modeling just as the case with in geographical information systems (GIS) analysis, though it should be understood that any model parameters cannot be interpreted simply in terms of physical meaning such as in hydraulics, but need to introduce the concept of equivalence of parameterization for the different model scales.

Integrating GIS and Hydrologic modeling
Topography plays an important role in the hydrologic response of a catchment to rainfall and has a major impact on the hydrological, geomorphological and biological processes active in that landscape (Fang et al., 2017). If meaningful hydrologic predictions are to be achieved at the landscape scale, the ability to characterize the spatial variability of hydrologic processes in a simple, yet physically realistic way is of major importance. The automation of terrain analysis and the use of DEMs has made it possible to quantify the topographic attributes of a landscape. One topographic attribute has proven to be particularly important in characterizing hydrologic processes: specific catchment area, which is an approximate measure of runoff per unit width and the convergence and divergence of flow (Yan et al., 2018). Specific catchment area together with other terrain attributes such as slope and profile and plan curvature, have been in different functional forms to describe the spatial distribution of zones of surface saturation, soil water content, runoff, evapotranspiration, erosion and deposition, and catenary soil development (Sommerlot et al., 2016;Tiwari et al., 2017). López-Vicente et al., (2017) explored the sensitivity of spatially distributed predictions of specific catchment area as a function of method of computation and DEM structure and the results indicated that the spatial characterization of hydrologic phenomena using GIS is very much methodologically dependent and that these methodological differences should not be ignored in environmental modeling and data base development.
GIS provides representations of the spatial features of the land surface while hydrologic modeling is concerned with the flow of water and its constituents over the land surface and the subsurface environment. There is obviously a close connection between the two subjects. Hydrologic modeling has been successful in the past in dealing with time progression, and models with many time steps are common, but the spatial disaggregation of the study area has been relatively crude (Shen, 2018). In many cases, hydrologic models assume uniform spatial properties or allow for small numbers of spatial sub-units within which properties are uniform (Widyastuti and Taufik, 2019;Yanto et al., 2017). GIS offers the potential to increase the degree of definition of spatial subunits, in number, in topology, and in descriptive detail, and GIS-hydrologic model linkage also offers the potential to address regional or continental scale processes whose hydrology has not been modelled previously to any significant extent (Tsanakas et al., 2016). Verma et al., (2017) recognized this potential and provided a comprehensive discussion on current state of hydrologic modeling independent of GIS and modeling coupled with GIS. And as new frontiers of linking GIS and hydrologic models that will make modeling more efficient and effective among others are in the followings:  Spatially distributed watershed properties  Partial area flow  Surface water -groundwater interaction  Regional and global hydrology  Spatial patterns of droughts Hydrologic phenomena are driven by rainfall and are thus always time dependent, even though by taking snap-shots at particular points in time or by time averaging over a long periods, a steady state model can be created (Conant et al., 2019). To accomplish a complete linkage between GIS and hydrologic model would require GIS to have time dependent data structure so that through time of the spatial distribution of hydrologic phenomena could be readily observed.

Runoff and Streamflow Simulation Models
Proliferation of recent runoff and streamflow simulation models has been based on account of physical characteristics of watersheds represented by topographic, geomorphologic, soils, vegetation, land use and land cover factors. One such model that is gaining much popularity and wide spread applications is TOPMODEL (Kirkby, 1986), which is a conceptual model based on variable contributing area with the predominant factors determining the runoff production process are represented by topography of the basin and a negative exponential law linking the transmissivity of the soil with the distance to the saturated zone below the ground level. Although conceptual, this model is described as `physically based model` in the sense that its parameters can be measured directly in situ or indirectly obtained from topographic and soil maps. Franchini et al., (1996) performed a detailed analysis to arrive at a better understanding of the correspondence between the model assumptions and the physical reality, in particular, the role of topographic information (topographic index) and the nature of the soil (saturated hydraulic conductivity and its decay with soil depth). Mengelkamp et al. (1997) described the development of a land surface scheme to model the surface energy and water balance (SEWAB) that included individual hydrologic processes and can simulate runoff generation on a wide range of spatial and temporal scales. Calibration and evaluation of the runoff generation processes in SEWAB had been done for small experimental catchment near Cork in Ireland and using data from Cabauw in the Netherlands (Mengelkamp et al., 2001). Local scale studies show that calculating runoff as saturation excess runoff can be appropriate on an annual time scale if net changes in soil moisture storage can be neglected. The hydrograph of the small Irish catchment was analyzed on a 20-minute time scale and characterized by an immediate response to individual rainfall events. The behavior is simulated by explicitly including the ponding and infiltration process for surface runoff generation. For macro-scale hydrologic model, SEWAB is used as a vertical component linked to a horizontal routing scheme (SEROS) and was implemented for large basin of the Odra drainage basin on grid size 18 km. Through calibration of the runoff generation process in SEWAB and of horizontal routing scheme hydrographs at various gauging stations were pro-duced. It is recognized here that land surface scheme of SEWAB needs to be coupled to a horizontal routing scheme to be implemented on large scale hydrologic model indicating the role and possible selection of alternatives of the routing schemes.

LARGE SCALE HORIZONTAL ROUTING MODELS
In previous paper (Lohmann et al., 1996;Nguyen-Quang et al., 2018;Piccolroaz et al., 2016;and Zhao et al., 2017), a large scale horizontal routing model was developed based on the unit hydrograph concept that is to be derived just from measured precipitation and streamflow data on a daily time step and this is to be coupled with land surface parameterization (LSP) scheme like SEWAB. In Lohmann´s paper, the linear transfer function theory was tested to compare the estimated effective precipitation with the runoff predicted by an LSP scheme which should be equal. Given a data series of input X(t) into a linear system and output Y(t) from that system it is in principle straightforward to find a linear tranfer function model connecting the two time series. This transfer function model is characterized by its impulse response function (IRF), called unit hydrograph (UH) by hydrologists. Even though, Lohmann recognized that this transfer process is strongly non-linear.
However, it was assumed that the routing model is linear, causal, stable and time invariant and more precisely expressed by the following equations: and 0 ≤ ≤ for every i With linear model formulation, as proposed by Mateo-Lázaro et al., (2018), where Q S (t) is the slow flow discharge or baseflow, Q F (t) is the fast flow or direct runoff and Q(t) is the total measured stream flow. Assuming a linear time invariant (LTI) relationship between fast component of streamflow and part of precipitation called effective precipitation Peff, a solution can be calculated for the impulse response function UH F for the fast flow and Peff. Both are determined by following integral equation which can be solved with an iterative procedure of equation (1), with t = T as the length of impulse response func-tion. UH(t) is the impulse response function (IRF) of the whole system with the condition ∫ ( ) = 1 0 .
And based on the theory of a cascade of linear reserviors, the Gamma function can represent the IRF of UH(t) as The storage constant k is the same for all n reservoirs. The parameter k and n are subject to calibration for each basin and are functions of basin characteristics.

NON-LINEAR ROUTING SCHEME REVISITED
The new concept of unit hydrograph was well developed few decades ago with the introduction of linear system theory by James Dooge and the nonlinear system techniques by Diskin and Boneh, (1972). The non-linear technique is a logical extention to the above linear approximation that can be done by introducing higher order terms to the impulse response functions (the kernels) as follows.
First, simplify the notations: IRF : UH(t) to become H(t) Input : Peff(t) to become x(t) Output : Q(t) to become y(t) Time : T approaches  Then, equation (8) as a first order approximation can be expressed as In a more general form, the impulse response system can adopt the Volterra expression as an analogy to the Taylor series expansion in the case of scalar function, as follows: The first term on the right represents the slow response system or baseflow and the second term is the linear response system, as commonly known as the unit hydrograph response function representing direct runoff, and the rest are the higher order-m non-linear response systems representing inter-flows of the hydrograph.
The objective of the analysis is, as the case of linear response system analysis, to determine the impulse response functions Hm(t1,t2,...,tm) for significant lower order m as parameter identification problem. One way to accomplish this is by what is called the orthogonal function approximation that is by introducing any known orthogonal function, such as Lagu-erre (continuous form) or Meixner (discrete form) functions. First, necessary properties of the kernels will be stated, then the orthogonal function approximations will be described next.

Properties of the kernels: the assumptions
Diskin and Boneh (1972) elaborated the kernel properties of second order functional series that may be adopted to represent the input-output relationships of a watershed system, provided certain conditions are placed on the kernels involved in the convolution integrals of the series. The conditions imposed on the kernels are due to the nature and definitions of the input and output functions and the general properties of the watershed system. The properties of the watershed system also restrict the range of integration of the convolution integrals in the series. The properties of kernel of the linear subsystem: and for kernel of the second order -non linear subsystem: For higher order non-linear kernels, the following assumptions was suggested to facilitate more efficient computations in estimating the kernels numerically using orthogonal functions approximation: 1. That all non-linear kernels are equal to the sum of products of linear kernel: Where Hq,j(tj)'s are linear kernels which are bounded and Lebesque integrable for all orders (functions of all real variables).
2. All non-linear kernels are symmetries. This assumption was also recognized by Diskin and Boneh (1972) for second order kernel as already indicated above.
The system is considered anticipating or physically realizable system satisfying the Volterra condition: Implying the finite functional relationship: and for non-negative time: t  0, and finite memory system, u, the last expression for first order and m-order responses, respectively take the forms: and for discrete equivalent simply replace d by unity and replace lowercase letters with corresponding capital letters, i.e: (32)

Orthogonal functions approximation
As stated previously, the objective of system identification is to determine the kernels of the response functions HM(S1, S2, …, SM). One way to accomplished this is using the orthogonal functions approximation to be described next.
Example of orthogonal polynomials: the Meixner functions, where: ) . ( )} ; = 0, 1, 2, … The case of first order kernel: assume that we can represent the kernel H1(S) by some linear expansion of the chosen orthogonal polynomials, Pi(S): The case of second or der kernel: H2(S1,S2).
By properties of non-linear kernels, 2 ( 1 , 2 ) = ∑ 1 ( 1 ). 2 ( 2 ) here each of the functions Hq1(S1) and Hq2(S2) are linear kernels that can be expanded using selected orthogonal polynomials as the case of first order kernel: And substituting these into the second kernel formula and re-arranging the summations results: and substituting this into the second order response functional relationship will result: and using the symmetry property of second order kernel, can be solved efficiently for the coefficients {i j}, therefore solving the second order kernel identification. Likewise for higher order kernels. Efficient computational procedures can then be developed to include calibration and verification steps, and ready for prediction.

CONCLUSIONS
1. The introduction of system concept into the complex hydrological processes of watershed system in the past decades has made it possible to develop hydrologic system models that incorporated all the driving factors in an integrating manner. With the system concepts integrating these complex processes have come into an understanding of river basins or watersheds as dynamic natural systems consisting of interconnected components and processes that form the water cycle. Within a basin, the dynamics of the hydrological processes are governed partially by the temporal and spatial characteristics of inputs and outputs and the land use/land cover conditions. New tools and concepts have been developed that permit a different analysis of processes (mathematical models, remote sensing and GIS, communication and telemetering) and allow watershed to be considered as a management unit. Recent hydrological modeling recognized the development of large scale horizontal routing schemes that are compatible to regional and global climate models. 2. The increased recognition of the importance of feedback mechanisms between land surface processes and climate system justified to develop more realistic land surface representations in the forms of large-scale hydrologic models that have reached a mature stage with the incorporation of TOPMODEL concept. This approach provides a basis for routing saturated flow down the length of the hillslope strip, and water balance accounting at each time step allows the evolution of the mean overall storage to be estimated. This can then be redistributed over the catchment at each time step on the basis of spatial distribution of an appropriate index that is derived from topographic and soil characteristics of each point in the catchment. 3. The development of a land surface scheme to model the surface energy and water balance (SEWAB) that included individual hydrologic processes and can simulate runoff generation on a wide range of spatial and temporal scales. For macro-scale hydrologic model, SEWAB was used as a vertical component linked to a horizontal routing scheme (SEROS) where linear transfer function theory was employed to estimate streamflow and compared with measured effective precipitation data. Through calibration of the runoff generation process in SEWAB and of horizontal routing scheme, hydrographs at various gauging stations were produced. It is recognized here that land surface scheme of SEWAB needs to be coupled to a horizontal routing scheme to be implemented on large scale hydrologic model indicating the role and possible selection of alternative of the routing schemes. It was recognized that the nature of the transfer process is highly non-linear. 4. The non-linear routing scheme that had been developed during the early days of the implementation of system theory into hydrologic systems analysis is revisited and considered appropriate as an alternative for modeling grid cell level of current large scale hydrologic modeling.

ACKNOWLEDGEMENT
This work was partly as a recollection of the first author that initially prepared during the visit of the first author to GKSS Research Center -Geesthacth Germany in the month of August 2001, where the SEWAB methodology was first introduced and being applied to the Brantas basin using HYDRO1K data grids. The hospitality of Dr. Heinz-Theo Mengelkamp and Colleagues were highly appreciated. The visit was made possible under the Travel Grant from Unesco Jakarta Office and the approval and support from Dr. Giusseppe Arduino were highly appreciated.