Transfer Function-Noise Modeling of the Dynamic Relationship between Botswana's Monthly Maximum Temperature and Evaporation

Corresponding Author: Oliver Moses Okavango Research Institute, University of Botswana, Maun, Botswana Tell: +267 6817203 Fax: +267 6861835 Email: omoses@ori.ub.bw Abstract: In Botswana, the rate of evaporation over open water surfaces is high because the country is characterised by hot and semi arid climatic conditions. The purpose of the study is to investigate the correlation between Botswana's monthly maximum temperature and evaporation and between monthly maximum wind speed and evaporation, for the stations Gaborone, Francistown and Maun and to fit transfer function-noise models between the highly correlated elements. The data were obtained from the Botswana Department of Meteorological Services. Cross correlation analysis revealed that for the selected weather stations, monthly maximum evaporation is much more strongly correlated to monthly maximum temperature than to monthly maximum wind speed. Monthly maximum temperature and evaporation were modelled with transfer function-noise models. The results of the study are crucial mainly to those in water management, agriculture and other similar industries.


Introduction
An unending circulation of water within the atmosphere starts with large quantities of water evaporating from water bodies (Ahrens, 1994). Evaporation depends on several factors, which include meteorological variables (Xu and Singh, 1998). Evaporation and the meteorological variables that affect it are of interest to various applications such as in water resources management, agriculture and climatology. Evaporation from waterbodies accounts for most of the moisture in the atmosphere (≈ 90%), while the remaining 10% is accounted for by transpiration and sublimation (Moran and Morgan, 1986;Ahrens, 1994;AMPJ, 2016). The contribution from transpiration is far greater than that from sublimation. Injection of moisture into the atmosphere is necessary for precipitation formation. Precipitation is the main way by which water is transferred from the atmosphere back to the surface of the earth. In regions with little rainfall, evaporation losses contribute considerably to the decrease of water levels in water bodies such as dams (Kumar et al., 2013).
Since meteorological variables influence evaporation, it is crucial to study dynamic relationships between these variables and evaporation. Dynamic relationships can be studied using Transfer Function-Noise Models (TFNMs). These models have been successfully employed in studying dynamic relationships between climatic and hydrological variables (Box and Jenkins, 1970;1976;Lungu, 1991;Lungu and Sefe, 1991;Box et al., 1994;Jain and Lungu, 2003;Manzione et al., 2010). When long series are used as input for TFNMs, output series of extensive length can be simulated (Manzione et al., 2010). The success in the application of these models in the past has motivated their application in this study. In this paper, these models are employed to study the dynamic relationships between evaporation and meteorological variables in Botswana, using data from the Botswana Department of Meteorological Services.
The data are for the stations Maun, Francistown and Gaborone. For this technique, it is crucial to select only input variables that have the strongest influence on evaporation for the constructed models to be simple enough to be understood (Condorena and da Costa, 2011). Literature reveals that there are four main meteorological parameters that affect evaporation, namely, solar radiation, wind speed, relative humidity and air temperature (Singh, 1992;Kumar et al., 2013). Based on the selection criteria of the input variables for the adopted technique and on what literature reveals about meteorological parameters in relation to evaporation, this study does not consider precipitation as one of the main meteorological parameters that affect evaporation. However, other authors have studied relationships between rainfall and extreme temperatures in Botswana  and between temperature and rainfall variability in South Africa (Nkuna and Odiyo, 2016). The present study found that out of the four main meteorological elements that affect evaporation, only air temperature and wind speed had reliable data. This study therefore investigates the correlation between Botswana's monthly maximum evaporation and temperature and between monthly maximum evaporation and wind speed and fits TFNMs between the strongly correlated climatic elements (Box and Jenkins, 1970;1976;Box et al., 1994;Huang and Wu, 2014;Yu et al., 2014;Deng et al., 2015;Park and Koo, 2015). Fitting of TFNMs between the stated climatic elements has never been done in Botswana; hence this study fills this gap.
It has been pointed out above that evaporation losses contribute significantly to the decrease of water levels in water resources. The water resources that are highly affected by evaporation and are of interest in this study are the country's Okavango Delta and some major dams which include Shashe, Letsibogo, Gaborone and Bokaa. These water resources provide fresh water to the nation. In addition, the Okavango Delta is one of the country's main tourist attractions. It is also a Ramsar and a World Heritage Site. The models constructed in this study can be used to predict the output series (evaporation) with a lead time of a few months, which will be valuable information, mainly to those responsible for the management of the limited water resources in Botswana, agriculture and other related works (Manzione et al., 2010;Kenabatho et al., 2015).

The Study Area
The study area is Botswana, which is a landlocked country (Fig. 1). It lies between latitudes 17 and 27 0 S and between longitudes 20 and 29 0 E. Drought is common and rainfall is highly variable in terms of both space and time. Rainfall occurs mainly during the period October to March. Annual mean rainfall ranges between a maximum of 650 mm in the north (Kasane) and a minimum of about 250 mm in the southwest (Bhalotra, 1984;1987). Mean monthly maximum temperatures range between 32 and 35°C in the northern half of the country in October and January, while mean monthly minimum temperatures range between 2 and 7°C in July in the southwest (Bhalotra, 1984;1987). Wind speeds have a low annual mean, which is about 11 km h −1 (Larsson, 1986), but daily maximum values can exceed 75 km h −1 . Botswana's climate can be summarised as hot and semi arid to arid (Pike, 1971;Bhalotra, 1984;1987;Moses, 2007). The country's weather conditions are influenced by synoptic scale weather systems such as the Atlantic Ocean high pressure cell, Indian Ocean high pressure cell, cut-off lows, frontal systems, surface lows, high pressure cells, Inter Tropical Convergence Zone (ITCZ) and westerly and easterly troughs (Moses and Parida, 2016).
The study area is represented by three meteorological stations located in three areas of national significance, which are Francistown, Maun and Gaborone. Maun is situated in the Okavango Delta (northwestern part of the country). The Delta is a source of fresh water, supports livelihoods and ecological functions in the northwestern region. It receives most of its water from Angolan high lands, but it loses approximately 98% of its water through evapotranspiration (USAID, 2013). Moreover, the Delta is a tourist site of high economic value. Tourism is the country's second largest economic sector after mining. Gaborone and Francistown are also areas of national significance since they are the country's first and second largest cities respectively. Dams shown in Fig. 1, which include Shashe and Letsibogo (near Francistown), Gaborone and Bokaa (Gaborone), have been constructed to supply these cities with fresh water. Evaporation from these water resources has persistently been a major problem, which limits their optimal performance (SMEC, 1987;Moses, 2007).

Data
Monthly maximum Temperature (MT), monthly Maximum Evaporation (ME) and monthly Maximum Wind speed (MW) observational data for the selected stations were provided by the Botswana Department of Meteorological Services. The geographical positions of these stations, their data types and periods are given in Table 1. The stations started their operations in the 1960's. Only the data periods indicated in the table had quality data. Periods later than those indicated in the table had huge data gaps and were therefore not considered in the analysis.

Transfer Function-Noise Models
TFNMs were employed in the study of the dynamic relationships between evaporation and meteorological elements that affect it. They were fitted to the highly correlated climatic elements, i.e., they were fitted either between MT and ME, or between MW and ME. Cross correlations, discussed in the sub section below, were used to determine the strengths of the correlation coefficients and to identify the TFNMs. To represent these models with an equation, suppose that pairs of observations (X t ,Y t ) are available at equispaced intervals of time and that X measures the level of an input to a system and that the level of X influences the level of a system output Y. When X is changed from one level to another, it has no immediate effect on the output, but produces a delayed response with Y eventually coming to equilibrium at a new level. Such a change is called a dynamic response, which is described by transfer function models. These models can be expressed as (Box and Jenkins, 1970;1976;Jenkins, 1979;Box et al., 1994;Huang and Wu, 2014;Park and Koo, 2015): where, the operators δ r (B) and ω s (B) are polynomials of orders (r,s), which are expressed respectively as: where, B is the backward shift operator, X t is a measure of the level of an input to a system which influences the level of a system output Y t , b is the delay parameter representing the number of complete time intervals before a change in X t begins to have an effect on Y t .
For stability, the roots of the characteristic equation δ r (B) = 0, with B regarded as a variable, must lie outside the unit circle. Equation 1 is of order (r,s,b). The order of this equation is identified from the data during the data based identification process (Romanowicz et al., 2010). The orders of r and s range from 0 to 2 (Box and Jenkins, 1970). In the real world, there are other influences that affect Y other than X, whose net effect is to corrupt the output predicted by the transfer function model. These other influences on Y, which also need to be accounted for by TFNMs, are called disturbances or noise. For a series that has been differenced (which removes trends and seasonality) to attain stationarity, TFNMs have the form (Box and Jenkins, 1970;1976;Jenkins, 1979;Box et al., 1994): With: where, n t is the net effect of the noise which corrupts the system at the output, yt = ∇ d Y t and x t = ∇ d X t are stationary processes with zero means, a t is the uncorrelated random variable, d is the order of nonseasonal differencing necessary to make the series stationary, ∇ = (1-B) is the backward differencing operator, ϕ(B) is the autoregressive operator, θ(B) is the moving average operator. Identification of TFNMs was simplified by prewhitening the input of the system. As the number one step of the identification process, the pre-whitening procedure of Box and Jenkins (1970) involves fitting an Autoregressive Moving Average (ARMA) model to the differenced input series. Autocorrelation Functions (ACFs) and Partial Autocorrelation Functions (PACFs) were used to fit the ARMA model to the differenced input series (Box and Jenkins, 1970;1976;Lungu, 1991;Lungu and Sefe, 1991;Box et al., 1994;Bierkens and Knotters, 1999;Jain and Lungu, 2003;Lungu et al., 2003;Huang and Wu, 2014;Kenabatho et al., 2015;Park and Koo, 2015). The 95% confidence limits of the autocorrelation and partial autocorrelations functions were approximated by ±2/n 1/2 (Harvey, 1993;Diggle and Chetwynd, 2011). An ARMA model has the form (Box and Jenkins, 1970;1976;Box et al., 1994): or when rearranged: where, α t is the uncorrelated input series. The rearranged ARMA model (Equation 4b) shows that the correlated input series x t is transformed to the uncorrelated white noise series α t by the transformation ϕ(B)θ −1 (B). The advantage of using the differencing method to attain stationarity is that if the same order of differencing is applied to both the input and the output series, the resulting transfer functions between the differenced series are identical to those between the original series. Assuming that the transformation applied to the input series x t can be applied to the differenced output series y t , yields: where, β t is the uncorrelated output series. The ARMA model represented by Equation 4a is a special case of the general multiplicative Autoregressive Integrated Moving Average (ARIMA) models, of order (p,d,q)x(P,D,Q)s, that account for both the non-seasonal and seasonal dependencies in the time series. These models have the general form (Box and Jenkins, 1970;1976;Box et al., 1994;Huang and Wu, 2014;Kenabatho et al., 2015;Park and Koo, 2015): where, ϕ p (B) is the autoregressive operator of order p, θ q (B) is the moving average operator of order q, D is the degree of seasonal differencing, s is the periodicity, ∇ s = (1-B s ) is the simplifying operator, Φ P (B s ) is the seasonal autoregressive operator of order P, Θ Q (B s ) is the seasonal moving average operator of order Q, z t is the series of observations.

Cross Correlations
Cross correlations at lag zero, which are a measure of the average correlations (Rummel, 1976), were used to determine the strengths of the relationships between the dependent and the independent time series, i.e., between the differenced series of MT and ME and between MW and ME. As mentioned above, TFNMs were fitted between the strongly correlated climatic elements. Cross correlations at various lags were used in the identification of the models, i.e., the orders (r,s,b) of the operators δ r (B) and ω s (B) in Equation 1, 2a and 2b and were also used to estimate the parameters of these operators. The advantage of cross correlations is that they can characterise lagged relationships between two series (LC, 2016) and can also discern relationships between variables that may possibly be difficult to discern because of strong trends in the data (Nkuna and Ndiyo, 2016). For the pre-whitened input series α t (from Equation 4b) and the corresponding output series β t (from Equation 5), the cross correlation between the two series at lag k, r αβ (k), is given by the equation (Chatfield, 2004;Montgomery et al., 2008): where, α and β are the averages of the α t and β t series respectively, while σ α and σ β are their standard deviations, n is the length of the two series, i is the i th value. Considering one series as the input and the other series as the output from a dynamic linear system, the delayed response is computed through the impulse response function (LC, 2016). The cross correlation at lag k is directly proportional to the impulse response function ν'(k) and is expressed as (Box and Jenkins, 1970;1976;Box et al., 1994): where, 2 s α and 2 s β are estimates of the variances 2 α σ and 2 β σ , respectively.
Approximate standard error of r αβ (k) were estimated by: Approximate standard errors provide a statistical check whether r αβ (k) is considerably different from zero (Pong-Wai, 1979).

Analysis of Mean Monthly Data
Mean monthly Maximum Temperature (mean MT), mean monthly maximum wind speed (mean MW) and mean monthly Maximum Evaporation (mean ME) were analysed to give an overview of their monthly variation. Variations of mean MT, mean ME and mean MW for Maun are shown in Fig. 2a-c respectively. Variations of mean monthly values of the other stations were similar to these figures. From Fig. 2a, it can be seen that mean MT vary from month to month, they are at their lowest in June and July and they are at their highest in January and December. Fig. 2b shows that mean ME also vary from month to month in a similar manner as mean MT in Fig. 2a. Figure 2c reveals that the lowest values of mean MW are attained in February, while the highest values are attained in October.

Cross Correlation Coefficient Analysis
The estimated cross correlation coefficients at lag zero, between MT and ME series and between MW and ME series, are shown in Table 2. The coefficients between MT and ME for Gaborone, Francistown and Maun are 0.66, 0.60 and 0.58, respectively. In contrast, the coefficients between MW and ME for Gaborone, Francistown and Maun are much smaller, being 0.28, 0.12 and 0.14 respectively. The estimated confidence limits of the coefficients are ±0.05 for Francistown and Gaborone and ±0.06 for Maun. To make it easier to interpret the coefficients, the percentages of variance in common between MT and ME and between MW and ME were computed by squaring the coefficients and then multiplying them by 100 (Rummel, 1976). Consequently, the percentages of variance in common between MT and ME for Gaborone, Francistown and Maun are 44, 36 and 34%, respectively. On the other hand, the percentages of variance in common between MW and ME for Gaborone, Francistown and Maun are 8, 1 and 2% respectively. The high coefficients between MT and ME imply stronger correlations between the variables, while the low coefficients imply weaker correlations. In other words, on the basis of the estimated magnitudes of the coefficients, MT has a much greater influence on ME than MW. However, all the correlation coefficients are positive, indicating positive associations between the variables. TFNMs were fitted between the highly correlated MT and ME series. The coefficients between MT and ME and between MW and ME are higher for Gaborone than for the other two stations (Table 2). This could be attributed to drier conditions in Gaborone compared to the other two stations; hence Gaborone has higher evaporation rates. Maun and Francistown are closer to higher rainfall areas of Angola and Zambia, which have a positive influence on their surrounding air moisture level.
A study by Kumar et al. (2013), which estimated evaporation from climatic factors in India, found correlation coefficients of 0.7625 and 0.6612 between MT and ME and between MW and ME, respectively. Clearly these coefficients are higher than those found in this study and they indicate that MW has almost the same level of influence on ME as MT in India. This implies that the level of influence of meteorological factors on evaporation differs from region to region.

Prewhitening of the Series
In the prewhitening process, the MT series was taken as the input series, while the ME series was taken as the output series. The first step was to fit an ARMA model (Equation 4b) to the MT series using Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF). Figure 3 is a plot of the ACF and PACF for Francistown's MT series. The estimated 95% confidence limits of the ACF and PACF are ±0.1. ACFs and PACFs for the other stations' input and output series resembled Fig. 3. This figure shows that the series were made stationary by the order of differencing d = 1 and D = 0. Since the autocorrelation at lag 0 is always equal to 1 (Shaw, 1983;Box and Jenkins, 1970;1976;Box et al., 1994), the figure shows that only the autocorrelation at lag 1 is statistically significant since it shoots beyond the 95% confidence limits. This suggests that for the MT series, the suitable ARMA model (Equation 4b) can be well represented by a first order autoregressive model: Comparing Equation 4b and 10 yields: Equation 11 is the transformation that was found to be suitable for transforming the input MT series to the uncorrelated pure white noise series α t , i.e., prewhitening the input series. In this equation, the estimated values of the autoregressive parameter φ that minimise 2 t α Σ are given in Table 3. By the assumption made in Equation 5, the transformation represented by Equation 11 was also applied to the output ME series to transform it to the uncorrelated white noise series β t .

Cross Correlation Function Analysis
Cross Correlation Functions (CCFs) were used to identify the delay parameter b, the orders r and s of the operators δ r (B) and ω s (B) in Equation 1, 2a and 2b and were also used to estimate the parameters of these operators. CCFs were between the uncorrelated α t series (prewhitened MT series) and β t series (prewhitened ME series). The variances 2 s α and 2 s β for the uncorrelated series α t and β t respectively, as well as approximate standard errors (Equation 9) of CCFs were computed and recorded in Table 3. CCFs for Francistown, Gaborone and Maun are plotted in Fig.  4a-c respectively. Figure 4a generally shows that the cross correlations smoothly decay towards zero with increasing lag. The cross correlation at lag 0 is the only one that is significantly greater than the approximate standard errors contained in Table 3. The cross correlations at lag 8 and 12 were considered as outliers, which were inevitable due to the fact that the same transformation that was used to transform the MT series was also used to transform the ME series. Francistown's CCF suggests that the delay parameter b in Equation 3a is equal to zero. The CCF for Gaborone ( Fig. 4b) also, generally indicates that the cross correlations smoothly decay towards zero with increasing lag. The cross correlations at lag 0, 1 and 9 overshoot the approximated standard errors, but the cross correlations at lag 1 and 9 were treated as outliers as already discussed in the case of Fig. 4a. Hence the delay parameter b for Gaborone was also taken to be zero. For the same reasoning as in the cases of Fig. 4a and 4b, the cross correlation at lag 0 in Fig. 4c (Maun's CCF) was considered as the only statistically significant cross correlation. As a result, the delay parameter b for Maun was also taken to be zero. The delay parameter b being zero implies that for the data considered in the analysis, there is no period of delay for a change in temperature to have an effect on evaporation.
The determined orders (r,s,b) were (1,1,0). This means that the parameters of the operators δ r (B) and ω s (B) were δ 1 , ω 0 and ω 1 . These parameters were estimated and recorded in Table 3. Having determined these parameters, the transfer function model in Equation 3a simplifies to:

The Noise Component
To determine the equation for the noise component n t , in Equation 12, the estimated values of the parameters δ 1 , ω 0 and ω 1 were used to generate the series: which was then used to estimate the n t series as (y t -y' t ).
Autocorrelation and partial autocorrelation functions of the n t series were computed and recorded in Table 4. They were used to identify the noise models for the selected stations. The identified noise models were autoregressive models of order 1, of the form: ( ) The estimated values of the autoregressive parameter ξ in Equation 14 are recorded in Table 5.

The Final Transfer Function-Noise Models
By Equation 12, 14 and the estimated parameters in Table 3 to 5, the transfer function-noise models fitted to the strongly correlated monthly maximum temperature and evaporation series for Francistown, Gaborone and Maun are Equation 15a-c respectively: The transfer function-noise models (Equation 15a-c) fitted to the data, were checked for adequacy using the criteria of Box and Jenkins (1970) and were found to be adequate. The models involve a small number of parameters, their cross correlation functions decay after a finite lag k (Fig. 4a to 4c) and absolute values of the autoregressive parameter ξ (Table 5) of the noise models are less than 1.

Correlations between
Botswana's monthly maximum temperature and evaporation and between monthly maximum wind speed and evaporation for Maun, Francistown and Gaborone were assessed. The correlations between monthly maximum temperature and evaporation were found to be much stronger than those between monthly maximum wind speed and evaporation. The more strongly correlated monthly maximum temperature and evaporation were modelled using transfer function-noise models given as Equation 15a-15c, for Francistown, Gaborone and Maun respectively. When checked for adequacy using the Box and Jenkins (1970) criteria, the models were found to be adequate. Such models can be used to predict the evaporation output series with a lead time of a few months. The results of the study are crucial primarily to those responsible for the management of the limited water resources in Botswana, agriculture and other related applications.