## 1. Introduction

Modelling time distribution of soil moisture is a key issue for evapotranspiration and biomass evaluation and is often adopted for deriving drought awareness indices. Water budget models help computing time evolution of soil moisture provided hydroclimatological and soil information. While runoff time series are often used to drive water budget models calibration, it may conduct to false conclusions about the other model outputs such as percolation and evapotranspiration fluxes, in the absence of vegetation response observations. Thus a lot of uncertainty is attached to the calibrated model parameters and may constitute a handicap against model application. The aim of this study is to propose a methodology to cope with vegetation information inside the calibration process of a water balance model using a qualitative approach. A review of evapotranspiration estimation through water balance modelling is reported in Kebaili Bargaoui (2011). In section 2, we present the data used to apply this methodology. In section 3, we present the methodology of uncertainty quantification using kernel distribution of model parameters. In section 4, resuling kernels are provided as well as a sensitivity study of results to the choice of soil parameters evaluation method.

## 2. Data

Two watersheds are studied: the Wadi Sejnane watershed (North Tunisia) and Wadi Chaffar watershed (South Tunisia). They are of comparable moderate sizes (respectively 376 km^{2} and 250 km^{2}). They have distinguishable occupation and climate. Sejnane basin is a forest basin under subhumid climate. Comparatively, for the Chaffar basin, vegetation cover comprises mainly olives under an arid climate. The soil type is principally sandy for Chaffar basin while a dominance of clay soils is outlined for Sejnane basin.

Potential evapotranspiration series are computed using the Turc formula based on monthly solar radiation and mean air temperature observed series at surrounding meteorological stations. A mean daily value is obtained for each month. Runoff (mm) series are estimated using observed daily stream discharges at the basin outlets with standard gauging methods. A ten year calibration period from September 1989 to August 1999 is considered for Chaffar basin including daily basin average rainfall evaluating using Thiessen method based on a network of 10 raingauges. A three year calibration period from September 1988 to August 1990 is available for Sejnane basin based on a rainfall network of 14 stations. Both basins have water tables. However piezometric data are not included in the study. Fig 1a and Fig. 1b report the time series of ETP, rainfall and runoff during the calibration periods.

## 3. Methodology

The water budget lumped BBH model presented by Kobayashi et al. (2001) is performed at daily time scale. Table 1 reports model equations. Mean daily rainfall and mean daily potential evapotranspiration are model inputs. Soil moisture content W (mm) and actual evapotranspiration ETR (mm) are model results of interest as well as runoff Rs (mm), percolation (Gd >0 mm) and capillary rise (Gd <0 mm). Seven parameters control model input-output transformations: thickness of active soil layer (*D* mm), effective soil porosity *p*; parameter related to the field capacity (*a* mm); parameter representing the decay of soil moisture (*b* mm); parameter representing the daily maximal capillary rise (*c* mm); parameter representing the moisture retaining capacity (0< η <1); parameter representing the stomatal resistance of vegetation to evapotranspiration (0< σ <1). The parameter W_{max} (mm) which represents the total water-holding capacity is a key parameter of the model.

According to Kobayachi and al. (2001) a/W_{max} is “nearly equal to or somewhat smaller than the field capacity”. After Teshima et al., (2006), *b* is a measure of soil moisture recession that depends on hydraulic conductivity and active soil layer depth *D*. In Iwanaga et al. (2005) a sensitivity analysis of BBH model applied to an irrigated area in semi-arid region suggest that soil moisture RMSE is most sensitive to σ, η and c. All parameters are subject to calibration using soil, vegetation as well as climatic and hydrologic information.

Moreover, we have introduced pedo transfer functions in the model in order to reduce the number of parameters to be calibrated on the basis of hydrometeorological series (Bargaoui and Houcine, 2010). It is worth noting that Kobayachi et al. (2001) adjusted soil humidity profiles measurements for BBH model calibration. As such observations are not often available; it seems an important task to adapt the original model using pedotransfer submodels especially when dealing with ungauged or partially gauged basins. To that purpose, three key soil characteristics are considered: saturated hydraulic conductivity *K*_{s}, soil water retention curve shape parameter *B* and field capacity S_{FC}. Assuming the percolation function as an exponential decay function, the leakage *L(s)* is identified according to Guswa et al. (2002) model as reported in (Eq. 1) where *s* is the ratio *W*/*W*_{max}. Consequently, parameters a, b, c of the original BBH model are obtained by identification (Eq. 2, 3, 4) using the three soil parameters *K*_{s}, S_{FC}, and *B*.

Rawls et al. (1982) model is adopted for estimating *K*_{s} while S_{FC} is derived according to two different models: Cosby and Saxton model which was recently adopted by Zhan et al., (2008) and Cosby et al. (1984) model. Effectively, this is suggested as a way to take into account uncertainty related to soil parameters. For the basin of Chaffar, because of lack of detailed information, we assume that *p* as well as *K*_{s} and S_{FC} parameters are those corresponding to the dominant soil class. For the Sejnane basin, a spatial mean of soil class properties is adopted using the spatial repartition of soil types as well as the area they cover within the basin. On the other hand, for the two cases, *B* =9 is adopted according to Rodriguez-Iturbe and al. (1999).

### 3.1. Model calibration

Finally, only the set of parameters (*D*, σ, η) remains subject to calibration through fitting observed and predicted runoff time series. The daily time step is adopted to run the model while annual, monthly and decadal time steps are adopted for its fitting. Many trials are firstly performed to adjust *D* choosing simply between three alternatives: *D* =1000 mm; *D* = 500 mm; *D* =300 mm which represent common values adopted in water balance models. Then, once *D* is fixed, the set of parameters (σ, η) is selected according to annual absolute relative runoff error AARE. Based on the idea of equifinality (Beven, 1993), a threshold value AARE_{ s} related to AARE is adopted for eliminating poor solutions using a grid of candidate solutions with Δ σ = Δ η= 0.01. Hence, only those pairs for which AARE > AARE_{ s} are selected and analyzed in the following.

Eq. (5) reports the objective function. It quantifies the absolute relative runoff bias during the calibration period:

where y_{oi} is the annual observed runoff (mm) for year *i*; y_{si} is the annual computed runoff (mm) for year *i*; N is the number of years of the calibration period. Additionally, for the selected solutions Nash coefficient R_{N,M} evaluated on the monthly basis as well as Nash coefficient R_{N,D} evaluated on the decadal basis are reported. The assumption of existence of capillary rise response is tested through the calibration process. It is further believed that if performance criteria (AARE, R_{N,M}, R_{N,D}) are better in presence (or in absence) of capillary rise assumption, then the assumption is retained.

### 3.2. Uncertainty quantification

Thus, the adoption of a fixed value for the threshold AARE_{s} will give rise to a number of acceptable solutions (σ, η). Here, the marginal kernel which represents a non parametric estimation of the statistical distribution of a given random variable (here the parameters σ and η) is adopted to represent parameter uncertainty. Similarly, the kernels of resulting outputs are computed in order to analysis the effect on model outputs especially evapotranspiration which is the variable of interest. A Gaussian kernel is adopted to perform the analysis.

### 3.3. Including vegetation information

It is now proposed to accurate (σ, η) kernel distribution by introducing the ratio Kv of mean annual actual evapotranspiration to mean annual potential evapotranspiration. In effect, as noticed by Eagleson (1994) after works of Ehleringer (1985), ecologists recognize three types of vegetation selection and adaptation in response to environmental stress due to water shortage (Type 1: desert annual grasses and humid climate trees; Type 2: semi-arid and sub-humid trees and shrubs; Type 3: perennial desert plants). Considering actual evapotranspiration as surrogate of vegetation productivity, three typical curves of Kv versus the inverse of the soil moisture are drawn by Eagleson (1994). Here, we assume the interval 0.45< Kv < 0.55 (mean Kv = 0.5) for type 2 (Sejnane basin) and 0.15< Kv < 0.25 (mean Kv = 0.2) for type 3 (Chaffar basin) which correspond to the values reported into the graph of Eagleson (1994) in case of weak environmental stress. Effectively, such an hypothesis is justified by the fact that the calibration periods represent mean water conditions for the two basins.

So, kernels of parameters and evapotranspiration conditional to the above conditions will also be drawn in order to evaluate the effect of including vegetation information supplementary to runoff observations on model results.

## 4. Results

Table 2 reports the 5 parameters which are not subject to fitting on basis of hydroclimatological series.

The thresholds AARE_{ s} = 5% and AARE_{ s} = 20 % have been applied respectively to Sejnane basin and Chaffar basin. Effectively, it was assumed that owing the more important time variability of rainfall and runoff for Chaffar basin series, it was more indicated to enlarge the threshold AARE_{ s} for this basin. The analysis of simulation results and runoff performance criteria relatively to the hypothesis of taking account or not for capillary rise (CR) results in not taking it into account for Chaffar basin (CR=0) while taking it into account for Sejnane basin (CR≠0).

### 4.1. Sejnane basin results

Fig. 2a reports the kernels corresponding to η in case of Sejnane basin. The result is sensitive to the choice of the parameter S_{FC} while kernels do not change with the change of the class of Kv in both assumptions on S_{FC}. Fig. 2b reports the resulting kernels of σ. It is worth noting that in both assumptions on S_{FC}, kernels are of uniform type reflecting the importance of uncertainty about σ. Conversely, the conditioning of results to the appropriate class of Kv (0.45 < Kv < 0.55) in relation with the vegetation and climate conditions of Sejnane watershed, reduces the uncertainty on σ and leads to two different intervals of variability for σ (smaller value of σ under the assumption of smaller value of S_{FC}).

More generally, the comparison of model error variances on monthly and decadal time scales suggests that the assumption of S_{FC} = 0.45 is more suitable for this basin. In effect, Fig. 3 which reports variances corresponding to the selected (σ,η) sets with AARE_{s}= 5% under the two assumptions on S_{FC}, shows that smaller variance values are achieved in the case where S_{FC} = 0.45.

Fig. 4 reports AARE values as well as the values achieved by the other performance criteria (R_{N,M} and R_{N,D}). Results are sorted according to Kv and dispatched by class of Kv. Four classes are considered: 0.35 < Kv < 0.45 ; 0.45 < Kv < 0.55 ; 0.55 < Kv < 0.65 ; Kv > 0.65. It is worth noting that the parameter sets (σ,η) which result in 0.45 < Kv < 0.55 exhibit the best AARE performance criteria while the other criteria are less sensitive to Kv conditions. Such a result might constitute a justification of adopting Eagleson (1994) Kv versus environmental stress condition variable within model fitting. Fig. 5 which reports the kernels of predicted actual evapotranspiration shows a clear reduction in uncertainty due to the inclusion of the constraint about vegetation and climate type. As well, it is noticeable that the kernel is less sensitive to the choice of S_{FC} when including such a constraint.

### 4.2. Chaffar basin results

Fig. 6a and
Fig. 6b report respectively R_{N,M} and R_{N,D} values obtained in case where S_{FC} =0.166 and CR≠0. They are reported according to the corresponding Kv. It is noticeable that Kv values with 0< Kv <1 result from such simulations. Negative values of R_{N,M} are often encountered suggesting very poor performances. Also, values of R_{N,D} are sometimes very low. Better results are obtained when assuming CR=0 (Fig. 7a and
Fig. 7b) with Kv values lying only in the interval (0.1< Kv < 0.2) which is more coherent with vegetation and climate information (type 3 curve).

Finally Fig. 8 reports η kernels in the two cases (CR=0 and CR≠0). It is noticed that the distribution of resulting evapotranspiration is sensitive to the model assumption about CR. The introduction of the constraint about Kv reduces a little the spread of the kernel distribution. The kernels of σ are reported in Fig. 9. It is noticeable that they are of uniform type in the interval (0,1) : U(0,1) in the case where CR=0 and U(0.5, 1) in the case where CR≠0. For the case CR=0, the constraint about Kv reduces the uncertainty and results in a uniform distribution U(0, 0.5). Fig. 10 reports the kernel distribution of evapotranspiration in the case CR=0. The constraint about Kv highly reduces the uncertainty about this output.

## 5. Conclusions

The methodology developed herein aimed to integrate the type of vegetation response within the calibration process of a water budget model at basin scale and daily time step. From developments using two different watershed of moderate size under two different climatic and vegetation conditions, it results in reducing the uncertainty about the parameters σ representing the resistance of vegetation to evapotranspiration and the parameter η representing the moisture retaining capacity. Hence, the uncertainty about actual evapotranspiration predictions has been also reduced due to such an analysis. This methodology is easily transferable to other water balance models as well as vegetation and climate situations.