Relating Hydro-Meteorological Variables to Water Table in an Unconfined Aquifer via Fuzzy Linear Regression

Abstract

This study aims to assess the short-term response of groundwater to the main hydro-meteorological variables of drought in a coastal unconfined aquifer. For this purpose, a multiple fuzzy linear regression-based methodology is implemented in order to relate rainfall, streamflow and the potential evapotranspiration to groundwater. Fuzzy regression analysis is recommended when there is a lack of data. The uncertainty of the system is incorporated into the regression coefficients which, in this study, are considered to be fuzzy symmetrical triangular numbers. Two objective functions are used producing a fuzzy band in which all the observed data must be included. The first objective function, based on Tanaka’s model, minimizes the total width of the produced fuzzy band. The second one includes the first while additionally minimizing the distance between the central value of the fuzzy output of the model and the observed value. Validity of the model is checked through suitability measures. The present methodology is applied at the east part of the Nestos River Delta in the Prefecture of Xanthi (Greece), where the observed values of the depth of groundwater level of four wells are examined.

1. Introduction

Groundwater is a vital resource of ecosystems and it is affected by natural and anthropogenic factors. Primarily, groundwater is the main water source in case of drought caused by the variability of precipitation and temperature. The groundwater level of unconfined aquifers is strongly influenced by the variability of the amount of rainfall of an area and the temperature [1,2,3]. Particularly, οn a local scale, these climate changes have greater impact on shallow aquifer systems than on deeper ones [3]. Additionally, shallow aquifers in lowland areas, which constitute an important role in the development of societies, are under more pressure since they are commonly associated with meeting the irrigation needs of crops. Different approaches have been developed to groundwater modeling in which mathematical relationships represent either the physical laws (physical models) or the natural processes of the system (conceptual models). Furthermore, natural processes may be represented through relations coming from the general theory of systems analysis without any consideration of physical laws and empirical relations (statistical-stochastic models).

Several studies, which intend to relate hydro-meteorological variables to groundwater, have been proposed. Viswanathan [4] suggests a multiple linear regression model in order to determine the recharge parameters of a coastal unconfined aquifer. Based on daily rainfall and water table records of the year 1979, he uses a recursive least squares method in order to minimize the difference between the groundwater level of an exploratory well and that estimated of the model. The water table level is estimated as a function of the given series of rainfall events and the prior water table levels. In their research, Ferdowsian et al. [5] investigate statistically the trends of groundwater levels by taking into account the time lag between rainfall and its impact on groundwater. Based on daily rainfall records, two forms of accumulative residual rainfall are estimated and compared in order to be used as one of the two independent variables in a multiple linear regression model. The other independent variable is the number of months since the commencement of the observation. The output (dependent variable) of the model is the water table depth, for which sufficient data from forty-nine (49) wells, both of shallow and deeper aquifers, are collected every three (3) or nine (9) months during a sampling period of 7 to 10 years. The regression coefficients estimated by the model represent the impact and the trend rate of the groundwater level rise or fall over time. In the study of Chen et al. [6], precipitation and temperature are linked to the groundwater level of a carbonate aquifer. On the basis of a groundwater flow and a water budget model, they proposed a statistical-empirical model applying a multiple linear regression in which groundwater level constitutes one function of meteorological variables with a time delay Δt. In another research [7], a multiple linear regression model between rainfall and groundwater is applied in order to investigate the groundwater response to rainfall. The model is powered by monthly records of water table depth and rainfall by considering their values from the previous monthly step. Data from several piezometric stations during the period of 2007 to 2008 are used and the analysis is conducted at representative areas in a regional scale basin. Zhang et al. [8] evaluated the effects of several factors with regard to the fluctuations in water table elevations in the case of shallow aquifers at a local scale. For this purpose, multiple/stepwise multiple regression techniques are used to investigate the linear relationship between identified groundwater level response height and independent factors. Six wells at an experimental site covering a region of 17 km² are used to collect water table elevation data every fifteen minutes on a weekly basis while hydrogeological data and site specific data are also used. In another experimental research [9], a statistical model is suggested for the investigation of the groundwater level response to precipitation, evaporation, river stage and tide level. Daily water table data are selected from twelve (12) wells in a shallow unconfined aquifer of a farmland covering an area of 50 × 150 m² for one year.

All the aforementioned researches are important and useful for the understanding of the interrelationships between the groundwater level and other hydrological and climate factors. The use of such models relies on the availability of sufficient long-term time series of the hydrogeological and hydrometeorological variables examined at each research. Furthermore, analysis in most of these researches is conducted on the basis of a river basin or a regional scale basin.

In the last two decades, fuzzy logic and sets have an increased use in surface water and groundwater hydrological applications examining various problems [10]. Several applications are based on fuzzy IF –THEN rules based methodologies [11,12,13,14], while others are based on hybrid fuzzy multicriteria techniques [15,16,17,18]. These applications are very popular because of the ANFIS toolbox of MATLAB. As far as fuzzy regression methodologies are concerned, a variety of applications werecarried out in modeling both qualitative and quantitative parameters [19,20,21,22,23,24,25,26,27,28,29,30,31].

Concerning the regression methodologies, the research of Bardossy et al. [19] was one of the first applications in groundwater hydrology and surface hydrology dealing with finding a relationship between saturated hydraulic permeability and resistivity. In addition several methodological points and criteria of vagueness were proposed by [19]. Furthermore, Bardossy et al. [32] developed a fuzzy unit hydrograph based model in order to calculate the fuzzy ordinates of the proposed fuzzy unit hydrograph. In the current study, a fuzzy linear relationship between rainfall (R), streamflow (Q), the potential evapotranspiration (PET) and the groundwater depth (D^GW) is investigated in the case of a shallow unconfined aquifer covering a local site of a deltaic environment. Two multiple linear regression models based on the principles of fuzzy logic and sets are implemented. Fuzzy linear regression models may satisfactorily function when there is a lack of data by incorporating uncertainty into the regression coefficients which consider them as fuzzy numbers [19,33]. All of the observed data must be included into the produced fuzzy band, the spread of which is minimized. Appropriateness of the models is checked through the value of the objective function which indicates the total fuzziness. Other suitability measures are also proposed in this work.

2. Materials and Methods

2.1. Basic Consepts of Fuzzy Logic and Sets

Fuzzy logic and sets consist of an extension of Boolean logic, which means that an element x of a given fuzzy set B may take not only two values {1,0}, but also all the values between 1 and 0 including limit values. Thus, an element x may belong to some degree to a fuzzy set B, whereas in the case of Boolean logic an element x either belongs or it doesn’t to a given set (crisp set). Fuzzy methodologies can be autonomous (for example, a set of fuzzy rules based on Mamdani approach, for example, [34]) or hybrid, where uncertainty of complex issues can be incorporated through analysis. Fundamental concepts of fuzzy logic and sets are referenced below:

A fuzzy set $\tilde{B}$ is a mapping $B:X\to \left[0,1\right]$. The membership function $\mu (x)$ of an element $x\in X$ indicates the degree of membership in the $\tilde{B}$.

A fuzzy number $\tilde{Y}$ defined on $ℝ$ is a special kind of fuzzy set satisfying the following properties [35]:

• $\exists x\in X$ such that $\mu (x)=1$ (normal fuzzy set)
• the α-cut, $Y_{[\alpha ]}$, must be a closed interval $\forall \alpha \in (0,1]$
• the support set (strong zero-cut), $Y_{{[0]}^{+}}$, of the fuzzy number $\tilde{Y}$ must be bounded

The α-cut, $Y_{[\alpha ]}$, of the fuzzy number $\tilde{Y}$ (and for any fuzzy set) is a crisp set containing all the elements in the X that have membership value in $\tilde{Y}$ greater than or equal to α:

$$Y_{[\alpha ]}=\{x\in X:Y(x)\ge \alpha \}$$

(1)

$$Y_{{[\alpha ]}^{+}}=\{x\in X:Y(x)>\alpha \}\left(\mathrm{strong}\alpha \text{-cut}\right)$$

(2)

It should be noted that the total fuzziness is taken into account when the strong zero cut, $Y_{{[0]}^{+}}$, is used. More analytically, according to Equation (2) above the 0-cut is an open interval and does not contain the boundaries. For this reason and in order to have a closed interval containing the boundaries, Hanss [36] proposed the phrase worst-case interval W, which is the union of the strong 0- cut and the boundaries.

A symmetric triangular fuzzy number (STFN) (Figure 1) is a special kind of fuzzy number of which the membership function is expressed by the following equation:

$${\mu }_Y(x)=\{\begin{array}{l}1 for x=r\\ \frac{x-(r-w)}{w} if r-w\le x\le r\\ \frac{(r+w)-x}{w} if r\le x\le r+w \\ 0 otherwise\end{array}$$

(3)

where w = the semi-width of the $\tilde{Y}$, r = the central value of the $\tilde{Y}$.

The extension principle is a fundamental principle in fuzzy set theory. In brief, with the use of the extension principle all the operations of the crisp functions can be extended as the fuzzy arithmetic and fuzzy algebraic operations [37,38]. Thus, a crisp function can be performed on a fuzzy number and the result of this operation should also be a fuzzy number.

Let $\tilde{Y}$ be a fuzzy number and let $f$ be a continuous crisp function. Then the result of the performed operation will be a fuzzy number $\tilde{Z}$ with a membership function ${\mu }_Z(\mathbf{x})$, the α-cuts of which can be described as follows [39]:

$$\begin{array}{l}{}^L{Z}_{[\alpha ]}=\min \{f(\mathbf{x}) | \mathbf{x}\in Y_{[\alpha ]}\},\\ {}^R{Z}_{[\alpha ]}=\max \{f(\mathbf{x}) | \mathbf{x}\in Y_{[\alpha ]}\}\end{array}$$

(4)

for 0 ≤ α ≤ 1, where ${}^L{Z}_{[\alpha ]}$ and ${}^R{Z}_{[\alpha ]}$ are the left boundary and the right boundary of the fuzzy number $\tilde{Z}$.

2.2. Fuzzy Linear Regression Analysis

2.2.1. Fuzzy Multiple Linear Regression based on Tanaka’s Model

In this research, a fuzzy multiple linear regression (Equation (5)) is implemented in order to relate the independent variables (R, PET and Q), to the dependent variable (D^GW). Uncertainty caused by the complexity of natural processes can be incorporated through the use of the fuzzy regression based on Tanaka’s approach [40]. This study examines the case of input and output data taken as experimental crisp values while the output of the model and the regression coefficients constitute STFN. Uncertainty is incorporated into the fuzzy regression coefficients, the determination of which leads to a conventional constrained optimization problem [41,42]. The objective function J of this optimization problem is minimized since it indicates the total fuzziness of the solution.

$${\tilde{D}}_j^{GW}= {\tilde{A}}_0 +{\tilde{A}}_R{R}_j+{\tilde{A}}_{PET}PET+ {\tilde{A}}_Q{Q}_j$$

(5)

where R, PET, Q are the crisp data of rainfall, potential evapotranspiration and streamflow, respectively, at the examined point in time j, and ${\tilde{A}}_j =(r_j,w_j)$ are the regression coefficients selected as symmetric triangular fuzzy numbers with central value r and semi-width w.

According to the extension principle the estimated fuzzy output, ${\tilde{D}}_j^{GW}$ will also be a symmetric triangular fuzzy number with a central value r and a semi-width w (Equation (6)).

$$\begin{array}{l}central value: D{ }_{r,j}^{GW}=r_0+r_R{R}_j+r_{PET}PE{T}_j+r_Q{Q}_j\\ semi-width: D{ }_{w,j}^{GW}=w_0+w_R\left|R_j\right|+w_{PET}\left|PE{T}_j\right|+w_Q\left|Q_j\right|\end{array}$$

(6)

Based on the concept of inclusion (Equation (7)), all the observed data must be included into the produced fuzzy band aiming at its minimum width.

$$D_j^{GW,obs}\subseteq {\tilde{D}}_{j,[h]}^{GW}$$

(7)

The concept of inclusion interprets the inclusion constraints (Equation (8)) of the optimization problem according to which an observation of groundwater depth ($D_j^{GW,obs}$) at the examined point in time j is included into the estimated fuzzy groundwater depth (${\tilde{D}}_j^{GW}$) with an associated degree $h\in \left[0,1\right]$. The level h denotes that the observation $D_j^{GW,obs}$ is contained in the support set of the corresponding estimated ${\tilde{D}}_j^{GW}$ with a membership degree greater than h [43].

$$\begin{array}{l}\left(r_0+r_R{R}_j+r_{PET}PE{T}_j+r_Q{Q}_j\right)-\left(1-h\right)\left(w_0+w_R\left|R_j\right|+w_{PET}\left|PE{T}_j\right|+w_Q\left|Q_j\right|\right)={}^L{D}_{j,[h]}^{GW}\le D_j^{GW,obs}\\ \left(r_0+r_R{R}_j+r_{PET}PE{T}_j+r_Q{Q}_j\right)+\left(1-h\right)\left(w_0+w_R\left|R_j\right|+w_{PET}\left|PE{T}_j\right|+w_Q\left|Q_j\right|\right)={}^R{D}_{j,[h]}^{GW}\ge D_j^{GW,obs}\\ w_0,w_R,w_{PET},w_Q\ge 0\end{array}$$

(8)

where ${}^L{D}_{j,[h]}^{GW}$ and ${}^R{D}_{j,[h]}^{GW}$ are the left boundary and the right boundary of the estimated ${\tilde{D}}_j^{GW}$, respectively, while $D_j^{GW,obs}$ is the observation of the groundwater depth at the examined point in time j.

Fuzzy regression analysis based on Tanaka [40] has no error term since the subject of the inclusion constraints is the minimization of the spread of the produced fuzzy band. The objective function J summarizes all the produced semi-widths for all points in time (j = 1, 2, …, k). Thus, suitability of the model is checked by the value of the objective function J (Equation (9), Fuzzy linear regression model I (FLR-1)). A small value of J indicates small fuzzy band and therefore the model has high suitability.

$$\min J=\min \left\{{\displaystyle \sum _{j=1}^k{w}_j^{GW}}\right\}=\min \left\{k{w}_0+{\displaystyle \sum _{j=1}^k\left(w_R\left|R_j\right|+w_{PET}\left|PE{T}_j\right|+w_Q\left|Q_j\right|\right)}\right\}$$

(9)

where k is the number of observations.

2.2.2. Modification of Tanaka’s Model with the Use of a Non-linear Objective Function

When using STFN, the objective function suggested above by Tanaka [40] minimizes the total semi-width for all the fuzzy regression coefficients and for all observations. Thus, under inclusion constraint conditions, it takes into account the distance between the left and the right boundary of each fuzzy output.

In their research, Tzimopoulos et al. [33] use a non-linear objective function based on the least squares model suggested by Diamond [44], which minimizes the distance between the observation (crisp output) and the left and the right boundary of the estimated fuzzy output. Tzimopoulos et al. [33] apply this non-linear objective function in the case of crisp experimental data without taking into account the inclusion constraints and thus, the regression coefficients result crisp numbers. It is worth mentioning that in the case of fuzzy experimental data the non-linear objective function used in the research of Tzimopoulos et al. [33] works well regardless of whether the inclusion constraints are used.

This research applies a modified fuzzy multiple linear regression model based on Tanaka [18] through the use of the non-linear objective function based on Tzimopoulos et al. [33] (Equation (10), Fuzzy linear regression model 2 (FLR-2)) by taking into account the inclusion constraints in the case of crisp experimental data. In contrast with the current study, in their research, Papadopoulos et al. [45] apply this modified fuzzy regression model between only one independent and one dependent variable and they work with probabilities.

$$\min S={\displaystyle \sum _{j=1}^k\{\begin{array}{l}{\left[D_j^{GW,obs}-\left(r_0-w_0+r_R{R}_j-w_R\left|R_j\right|+r_{PET}PET-w_{PET}\left|PE{T}_j\right|+r_Q{Q}_j-w_Q\left|Q_j\right|\right)\right]}^2+\\ {\left[D_j^{GW,obs}-\left(r_0+w_0+r_R{R}_j+w_R\left|R_j\right|+r_{PET}PET+w_{PET}\left|PE{T}_j\right|+r_Q{Q}_j+w_Q\left|Q_j\right|\right)\right]}^2\end{array}}$$

(10)

where the first bracket denotes the Euclidian distance between the observed groundwater depth, $D_j^{GW,obs}$, and the left boundary of the corresponding produced fuzzy output of the fuzzy regression model, ${}^L{D}_j^{GW}$, for level h = 0. The second bracket denotes the Euclidian distance between the $D_j^{GW,obs}$ and the right boundary of the corresponding produced, ${}^R{D}_j^{GW}$, for level h = 0.

Based on the consideration that the total produced fuzziness can be analyzed either around the central values or the observations (Equation (11)) (similar considerations can be found in [46] which aim at different purposes), it is concluded (Equations (11)–(14)) that the objective function S (Equation (10)) takes into account both the distance between the central value and the left-right boundaries, and the distance between the central value and the observation.

$${{\displaystyle \sum _{j=1}^k\left[\left({}^R{D}_j^{GW}-D_j^{GW,obs}\right)+\left(D_j^{GW,obs}-{}^L{D}_j^{GW}\right)\right]}}^2={{\displaystyle \sum _{j=1}^k\left[\left({}^R{D}_j^{GW}-D{ }_{r,j}^{GW}\right)+\left(D{ }_{r,j}^{GW}-{}^L{D}_j^{GW}\right)\right]}}^2$$

(11)

Analyzing Equation (11) and solving for ${\displaystyle \sum _{j=1}^k{\left({}^R{D}_j^{GW}-D_j^{GW,obs}\right)}^2+{\displaystyle \sum _{j=1}^k{\left(D_j^{GW,obs}-{}^L{D}_j^{GW}\right)}^2}}$ and developing the analytical expression from the right, it is concluded, as shown in the following equation [45,46].

$$\begin{array}{l}{\displaystyle \sum _{j=1}^k{\left({}^R{D}_j^{GW}-D_j^{GW,obs}\right)}^2+{\displaystyle \sum _{j=1}^k{\left(D_j^{GW,obs}-{}^L{D}_j^{GW}\right)}^2}}=\\ {\displaystyle \sum _{j=1}^k{\left({}^R{D}_j^{GW}-D{ }_{r,j}^{GW}\right)}^2}+{\displaystyle \sum _{j=1}^k{\left(D{ }_{r,j}^{GW}-{}^L{D}_j^{GW}\right)}^2}+2{{\displaystyle \sum _{j=1}^k\left(D{ }_{r,j}^{GW}-D_j^{GW,obs}\right)}}^2\end{array}$$

(12)

The analytical expression from the left denotes the total sum of the difference between the $D_j^{GW,obs}$ and ${}^R{D}_j^{GW}$ raised to the second power plus the total sum of the difference between the $D_j^{GW,obs}$ and ${}^L{D}_j^{GW}$ raised to the second power. The first and the second terms of the right analytical expression are the total sum of the difference between the boundaries (right and left, respectively) and the central value of the fuzzy groundwater depth ($D{ }_{r,j}^{GW}$) raised to the second power. The last term denotes the (double) total sum of the difference between the $D{ }_{r,j}^{GW}$ and $D_j^{GW,obs}$.

In the case of STFN (Figure 1) the distance between the central value and the right boundary is equal to the distance between the central value and the left boundary (Equation (13)), thus the first and the second terms of the right part of the Equation (12) can be written as follows:

$${\displaystyle \sum _{j=1}^k{\left({}^R{D}_j^{GW}-D{ }_{r,j}^{GW}\right)}^2}={\displaystyle \sum _{j=1}^k{\left(D{ }_{r,j}^{GW}-{}^L{D}_j^{GW}\right)}^2=}{{\displaystyle \sum _{j=1}^{k}D{ }_{w,j}^{GW}}}^2={\displaystyle \sum _{j=1}^k{\left(w_0+w_R\left|R_j\right|+w_{PET}\left|PE{T}_j\right|+w_Q\left|Q_j\right|\right)}^2}$$

(13)

where ${\displaystyle \sum _{j=1}^{k}D{ }_{w,j}^{GW}}$ denotes the total semi-spread of the produced fuzzy band.

Hence, Equation (12) can be re-written as follows:

$$\begin{array}{l}{\displaystyle \sum _{j=1}^k{\left({}^R{D}_j^{GW}-D_j^{GW,obs}\right)}^2+{\displaystyle \sum _{j=1}^k{\left(D_j^{GW,obs}-{}^L{D}_j^{GW}\right)}^2}}=\\ 2{\displaystyle \sum _{j=1}^k{\left(w_0+w_R\left|R_j\right|+w_{PET}\left|PE{T}_j\right|+w_Q\left|Q_j\right|\right)}^2}+2{{\displaystyle \sum _{j=1}^k\left(r_0+r_R{R}_j+r_{PET}PET+r_Q{Q}_j-D_j^{GW,obs}\right)}}^2\end{array}$$

(14)

Consequently, the objective function S includes both the total fuzziness and the (double) distance between the central value of the fuzzy groundwater depth (${\tilde{D}}_j^{GW}$) and the observed groundwater depth ($D_j^{GW,obs}$).

2.3. Suitability Measures

As aforementioned above, validity of the applied fuzzy linear regression models is checked through the values of the objective functions J and S. Lower values of J and S indicate higher suitability of the models. In addition, two more suitability measures are used.

The first one is based on Theil’s inequality coefficient U [47], as presented according to Bliemel [48] (the first of the two formulae) and Botzoris and Papadopoulos [49]. The difference with the conventional U is that the fuzzy output estimated by the two fuzzy linear regression models (which in this study is the fuzzy groundwater depth) is used instead of the (crisp) output estimated by the conventional regression. Based on the extension of principle [36,39], algebraic operations between fuzzy numbers and crisp numbers can be performed producing fuzzy numbers and hence, a fuzzification of U is achieved. Thus, this study proposes, for the first time, the fuzzified Theil’s inequality coefficient, $\tilde{U}$, which is expressed as follows:

$$\tilde{U}=\frac{\sqrt{\frac{1}{k}{{\displaystyle \sum _{j=1}^k\left({\tilde{D}}_j^{GW}-D_j^{GW,obs}\right)}}^2}}{\sqrt{\frac{1}{k}{{\displaystyle \sum _{j=1}^k{\tilde{D}}_j^{GW}}}^2}+\sqrt{\frac{1}{k}{{\displaystyle \sum _{j=1}^k{D}_j^{GW,obs}}}^2}}$$

(15)

where ${\tilde{D}}_j^{GW}$ is the fuzzy output (fuzzy groundwater depth which is produced based on the fuzzy linear regression model) and, $D_j^{GW,obs}$ is the observed groundwater depth at the examined point in time j (which is crisp number).

The problem of Equation (15) is that the fuzzified Theil’s inequality coefficient, $\tilde{U}$ contains some fuzzy inputs, that is, the fuzzy prediction of the water depth ${\tilde{D}}_j^{GW}$. This point can be addressed based on the extension principle. Finally, based on Equation (4), the left (minimum) and the right (maximum) hand sides of the fuzzified Theil’s inequality coefficient, $\tilde{U}$ are calculated as follows for one α-cut:

$$\begin{array}{l}{}^L{U}_{[\alpha ]}=\min \{f(x_1,\dots ,x_k) | x_1\in D_{1 [\alpha ]}^{GW},\dots ,x_k\in D_{k [\alpha ]}^{GW}\},\\ {}^R{U}_{[\alpha ]}=\max \{f(x_1,\dots ,x_k) | x_1\in D_{1 [\alpha ]}^{GW},\dots ,x_k\in D_{k [\alpha ]}^{GW}\}\end{array}$$

(16)

where $D_{j [\alpha ]}^{GW}$, for 0 ≤ α ≤ 1, denotes the α-cuts of the fuzzy groundwater depth and x₁, x₂, …, x₉ are the solutions of the double optimization problem. By using a significant number of α-cuts the fuzzified Theil’s inequality coefficient can be built.

Theil’s U can take values between zero and unit [0,1]. It is always reasonable that their values be close to zero [48,49]. The two fuzzy linear regression models produce another fuzzy groundwater depth value (fuzzy outputs) and therefore another $\tilde{U}$ are estimated for each model. However, since the Theil’s U are produced as fuzzy numbers, it is demanded to choose which of them is greater or smaller in order to decide the more suitable fuzzy regression model. Therefore, a computationally efficient method to compare the fuzzy numbers is presented in the Appendix A.

In the second suitability measure E_dis, which is used in the research of [45], the numerator denotes the Euclidian distance between the observed groundwater depth and the left boundary, the right boundary and the central value of the corresponding fuzzy estimate. The denominator denotes the distance between the observation and the unbiased mean of the historical sample. It is calculated through the following algebraic expression:

$$E_{dis}=1-\frac{{\displaystyle \sum _{j=1}^k\left({\left(D_j^{GW,obs}-{}^L{D}_j^{GW}\right)}^2+{\left(D_j^{GW,obs}-{}^R{D}_j^{GW}\right)}^2+{\left(D_j^{GW,obs}-D{ }_{r,j}^{GW}\right)}^2\right)/3}}{{\displaystyle \sum _{j=1}^k{\left(D_j^{GW,obs}-\overline{D^{GW}}\right)}^2}}$$

(17)

where ${}^L{D}_j^{GW}$, ${}^R{D}_j^{GW}$ and $D{ }_{r,j}^{GW}$ are the left boundary, the right boundary and the central value of the fuzzy groundwater depth, respectively, while $\overline{D^{GW}}$ is the mean value based on the historical sample. The closer E_dis is to unit, the better the model.

3. Results

3.1. Case Study

The study area is located at the south part of the City of Xanthi in the Prefecture of Xanthi, N.E. Greece. It is bounded to the west by the Nestos River and to the south by the Aegean Sea. Its geomorphology is generally considered to be flat with an elevation of a few meters above sea level throughout the entire study area (Figure 2). It is located in a recent sedimentary delta environment of a few tens of meters thick alternate sand, clay and silt layering deposits. It is worth mentioning the occasional presence of organic clay due to the delta marshes. The evolution of the east part of the delta under flooding conditions has been instrumental in forming low potential aquifers in the study area [50,51]. After carefully studying drilling, piezometric and geophysical exploration data of the area, it was concluded that at the north side of the study area, alternate clay and mostly sand layering extend down to a depth of 30 m. A marly layer, 50 m thick, comes in between 30 m to 80 m and below the depth of 80 m, the same clay and sand layers extend again [52,53,54].

The east delta plain extends for 176.4 km², forwhich 106.63 km² are cultivated, while the coastal saline uncultivated lands extend for 45 km². The irrigated lands extend for 89.9 km², while 35 km² of them meet irrigation needs from the Nestos River. The remaining areas meet irrigation needs by pumping groundwater from the unconfined aquifer of the area. The mean annual water consumption is estimated at 27 × 10⁶ m³ [51].

Hydrogeologically speaking, there is a shallow and a deeper aquifer system both of which are formed within the alluvial deposits of the wider study area [53,54]. This study focuses on the shallow hydrological system, which consists of phreatic and of semi-confined aquifers extending to a depth of about 30 m. The transmissivity (T) value of the unconfined aquifer is approximately 1.1 × 10⁻² m²/sec [53,54].

As aforementioned, the inputs of the fuzzy regression models are rainfall (R), the potential evapotranspiration (PET) and streamflow (Q). With respect to R and the PET, their monthly measurements, which refer to the corresponding observations of the groundwater depth (D^GW,obs) for the month j and for the period of October 2006 to October 2008, that is, k = 9, are utilized. It is mentioned that PET is calculated based on the Thornthwaite method. Streamflow data are derived from the three-months mean value, starting from October 2006. As far as the groundwater depth of the shallow unconfined aquifer is concerned, which is the output of the two fuzzy multiple linear regression models, the used measurements were collected from four (4) wells every three months, starting from October 2006 [52].

3.2. Results of the Two Fuzzy Regression Models

Αs aforementioned, for simplicity, the fuzzy linear regression model based on Tanaka’s approach [40] will be symbolized as FLR-1 and that which uses the non-linear objective function S will be symbolized as FLR-2. Likewise, the estimated $\tilde{U}$ based on the results of FLR-1 will be symbolized as ${\tilde{U}}_1$ and the estimated $\tilde{U}$ based on the results of FLR-2 will be symbolized as ${\tilde{U}}_2$.

The results of the FLR-1 model are separately presented in Figure 3 in the case of well 194. For illustration purposes, only R-D^GW and PET-D^GW are presented. The coefficients of the fuzzy regression are presented in

Table 1. The fuzzy coefficients regarding the examined wells based on the FLR-1 model.

Wells	Constant Term		R		PET		Qmean
	$r_0$ (centre)	$w_0$ (semi-width)	$r_R$ (centre)	$w_R$ (width)	$r_{PET}$ (centre)	$w_{PET}$ (width)	$r_Q$ (centre)	$c_Q$ (width)
177	0.7371	0.4752	−0.0048	0	0.0456	0.0101	−0.0281	0
183	1.4440	0.6340	−0.0044	0	0.0416	0	−0.0384	0
186	1.2953	0.6977	−0.0021	0	0.0403	0	−0.0334	0
194	2.2853	0.1977	−0.0011	0	0.0260	0	−0.0442	0.0075

.

As it can be easily observed in the above Figure 3, all the observed data are included into the produced fuzzy band as required by the inclusion constraints. Additionally, as indicated in Table 1, most of the coefficients result incrisp numbers except from the constant term. Particularly, as expected, the coefficients of R and Q are negative while the coefficient of PET has a positive symbol.

Similar findings regarding the regression coefficients are obtained when using the FLR-2 model (

Table 2. The fuzzy coefficients regarding the examined wells based on the FLR-2 model.

Wells	Constant Term		R		PET		Qmean
	$r_0$ (centre)	$w_0$ (semi-width)	$r_R$ (centre)	$w_R$ (width)	$r_{PET}$ (centre)	$w_{PET}$ (width)	$r_Q$ (centre)	$w_Q$ (width)
177	0.7790	0.4743	−0.0059	0	0.0452	0.0103	−0.0283	0
183	1.4824	0.6416	−0.0058	0	0.0413	0	−0.0384	0
186	1.3093	0.7005	−0.0026	0	0.0402	0	−0.0334	0
194	2.2438	0.2914	−0.0007	0	0.0259	0	−0.0432	0.0020

), while the observations of D^GW and its corresponding fuzzy estimates are illustrated (for the case of well 194) in Figure 4.

Furthermore, as it is observed in Table 1 and the Table 2, in both of models the fuzzy coefficient of streamflow seems to be increasing (in terms of absolute value) as the coast is approached, whereas the values of the rainfall’s modeldecreases (in terms of absolute value). Meanwhile, the total fuzziness (J) gradually decreases as the coast is approached, since it takes the highest value in the case of well 177 (the furthest well from the coast) and the lowest value in the case of well 194 (the shortest well from the coast) (

Table 3. Suitability measures and the values of the objective functions J and S for both of the fuzzy linear regression models.

Results based on the FLR-1 model
Wells	J	S	LU	CU	RU	Edis
177	10.294	46.779	0.000	0.1212	0.2803	0.800
183	5.706	12.348	0.000	0.0653	0.1401	0.905
186	6.279	14.293	0.000	0.0700	0.1498	0.874
194	2.847	3.229	0.000	0.0379	0.0823	0.937
Results based on the FLR-2 model
Wells	J	S	UL	UC	UR	Edis
177	10.371	46.650	0.000	0.1196	0.2801	0.802
183	5.774	12.175	0.000	0.0633	0.1397	0.908
186	6.304	14.272	0.000	0.0695	0.1509	0.875
194	2.899	3.154	0.000	0.0370	0.0816	0.939

). In both of the fuzzy regression models, as it easily observed in Table 3, the lowest value of $\tilde{U}$ and the higher value of E_dis are also obtained in the case of well 194. The previous ascertainments can be justified in hydrological terms considering the hydrogeological characteristics of the aquifer under study focusing on the piezometric conditions and the depth values of the groundwater level in the area of the wells 177, 183, 186, and 194 (Figure 2) [52].

In the Table 3 below, the bold values constitute the lowest values of the objective functions J and S and the lowest and highest values that the suitability measures (U and E_dis) correspondingly get in both fuzzy regression models. The fuzzified Theil’s coefficient of inequality, $\tilde{U}$, is described by the left boundary (^LU), the central value (^CU) and the right boundary (^RU).

It is worth mentioning that the Pearson’s r between D^GW-R, D^GW-PET and D^GW-Q, ranges (regarding all the examined wells) from −0.407 up to −0.497, from 0.931 up to 0.955 and from −0.466 up to −0.538 correspondingly.

As aforementioned, in both of the fuzzy regression models, $\tilde{U}$ gets the lowest value in the case of well 194. However, it is difficult to choose the lowest $\tilde{U}$ between the applied two fuzzy linear regression models (for each well) since the $\tilde{U}$ has a fuzzified shape. Hence, the fuzzy numbers comparison method of the Appendix A is used, which is based on the measure R. The values of the ranking measures R (R-values) are presented in

Table 4. R-values of the two fuzzified Theil’s inequality coefficients ${\tilde{U}}_1$ and ${\tilde{U}}_2$.

Ranking Measure/Wells	177	183	186	194
R (for ${\tilde{U}}_1$)	0.1277	0.0669	0.0717	0.0391
R (for ${\tilde{U}}_2$)	0.1266	0.0652	0.0717	0.0384

.

In the above Table 4, the lowest (bold) values indicate that in the case of well 194, ${\tilde{U}}_1$ and ${\tilde{U}}_2$ get the lowest value. In addition, the table information shows that ${\tilde{U}}_2$ is lower than ${\tilde{U}}_1$. The lower the value of R, the closer the membership function of each $\tilde{U}$ becomes to 0 (Figure 5).

In the above Figure 5, the Theil’s inequality coefficients U are illustrated together in the case of well 194. ${\tilde{U}}_2$ (the red shape) is lower than ${\tilde{U}}_1$ (the blue shape) and hence, the FLR-2 model which uses the non-linear objective function S is a little more suitable. Furthermore, ${\tilde{U}}_2$ is lower than ${\tilde{U}}_1$ regarding all the examined wells, while FLR-2 gets a higher E_dis-value than the FLR-1 in all cases of the examined wells.

4. Discussion Points

Based on the results above (Table 1 and Table 2), two fuzzy linear relationships are produced, one for each multiple fuzzy linear regression model (FLR-1 and FLR-2). As expected, rainfall (R) and streamflow (Q) negatively affect the groundwater depth (D^GW), while the potential evapotranspiration (PET) is positively related to D^GW. In both of the two fuzzy regression models, the fuzzy regression coefficient of the streamflow seems to have an inversely proportional relation with the distance of the coast, while the fuzzy regression coefficients of the rainfall grows lower. This could be explained by the geology of the case study since grain size is increased approaching to the coast and therefore the hydraulic interfacebetween the aquifer and the Nestos River may be stronger. This assumption is consistent with the piezometric conditions of (Figure 2) and the observed depths of groundwater of the case study.

In addition, total fuzziness continues to decrease as the distance tothe coast decreases, that is, the FLR-1 and FLR-2 models work better approaching the coast. As can be easily observed in Table 3, in the case of well 194, the total fuzziness J gets its lowest value with respect to the other wells in both fuzzy regression models. Simultaneously, the Theil’s U and the suitability measure E_dis get their lowest and highest values, respectively. Therefore, there is a robust linear (fuzzy) relationship between the aforementioned hydrometeorological variables in the case of well 194. This fact could also be explained by the geology of the case study where significant variability of hydraulic characteristics of the aquifer appears [52]. For instance, a lower permeability in places can negatively affect the short-term influence of streamflow and rainfall on groundwater.

Regarding the suitability measures,the value of Theil’s inequality coefficient U increases when the total fuzziness increases as well, while the higher the values of the suitability measures E_dis are, the lower the fuzziness of the models. Consequently, the suitability measures are consistent. It is observed that the uncertainties of the influence coefficients are negligible, especially for streamflow and rainfall, despite their correlation coefficients being much lower than the corresponding ones of potential evapotranspiration. However, both the higher correlation coefficient of streamflow than the corresponding one of rainfall in all examined wells and the coefficients of streamflow in both of the two fuzzy regression models, show that the contribution of streamflow should not be ignored. Eventually, fuzzy linear relationships instead of crisp relationships are produced.

It is highlighted that in general, both the fuzzy regression models (FLR-1and FLR-2) obtain similar results. However, it is worth commenting on the fact that although the total fuzziness (objective function J) has alower value when using the FLR-1 model (which is obvious, since for the FLR-1 model, J is the objective function), all of the measures, U and E_dis, are a little more suitable for the FLR-2 model regarding all the examined wells. Since the solution of the FLR-2 model has a better performance according to the majority of the criteria and because of the fact that the objective function of the FLR-2 model includes both the distances between central values-observed values and the objective function of FLR-1 (Equation (14)), the authors propose that the use of the FLR-2 model shall be preferred.

In that point, it is desirable to point out the difference between the crisp-fuzzy linear relationship. In fuzzy linear regression based on the Tanaka model the observed depth of groundwater (D^GW) is included in the support of the corresponding fuzzy estimate and hence, a membership degree corresponds to it. Whereas, statistical regression creates a single line in which the observations are either coincident or they are not. This is shown in the following example:

Let us consider the point with the black cycle (4th point) of the multiple fuzzy linear regression based on the Tanaka model (FLR-1) performed in the case of well 177 and illustrated in Figure 6 below.

The membership function of the fuzzy depth of groundwater is described by Equation (3). Based on Equations (6) and (8), Equation (3) can be rewritten as follows:

$$\begin{array}{l}{\mu }_Y(D_j^{GW})=\{\begin{array}{l}1-\frac{-D_j^{GW,obs}+\left(r_0+r_R{R}_j+r_{PET}PE{T}_j+r_Q{Q}_j\right)}{\left(w_0+w_R\left|R_j\right|+w_{PET}\left|PE{T}_j\right|+w_Q\left|Q_j\right|\right)} \\ if \left(r_0+r_R{R}_j+r_{PET}PE{T}_j+r_Q{Q}_j\right)-\left(w_0+w_R\left|R_j\right|+w_{PET}\left|PE{T}_j\right|+w_Q\left|Q_j\right|\right)\le D_j^{GW,obs}\le \left(r_0+r_R{R}_j+r_{PET}PE{T}_j+r_Q{Q}_j\right)\\ \\ 1-\frac{D_j^{GW,obs}-\left(r_0+r_R{R}_j+r_{PET}PE{T}_j+r_Q{Q}_j\right)}{\left(w_0+w_R\left|R_j\right|+w_{PET}\left|PE{T}_j\right|+w_Q\left|Q_j\right|\right)} \\ if \left(r_0+r_R{R}_j+r_{PET}PE{T}_j+r_Q{Q}_j\right)\le D_j^{GW,obs}\le \left(r_0+r_R{R}_j+r_{PET}PE{T}_j+r_Q{Q}_j\right)+\left(w_0+w_R\left|R_j\right|+w_{PET}\left|PE{T}_j\right|+w_Q\left|Q_j\right|\right) \\ 0 otherwise\end{array}\\ w_0,w_R,w_{PET},w_Q\ge 0\end{array}$$

(18)

where the brackets of the nominators denote the central value of the fuzzy depth of groundwater and the brackets of the denominators denote its semi-width. The up and down inequalities denote the left and right boundaries, respectively, while the term $D_j^{GW,obs}$ is the observation j of groundwater depth.

In the case of well 177 it holds: Y = 1.01, X1 = 150.60, X2 = 49.75, x3 = 33.66, where X1, X2, and X3 are the measurements of rainfall, potential evapotranspiration and streamflow, respectively, while Y is $D_j^{GW,obs}$. Based on the results of fuzzy linear regression of Tanaka (FLR-1), the left hand boundary is equal to 0.3563 m, the right hand boundary is equal to 2.3147 m and the central value is equal to 1.3370 m. Replacing these values into Equation (18), the membership degree which corresponds to the value of the observation can be estimated:

$$1-\frac{\left(-1.01+1.3370\right)}{0.9777}=0.6655$$

Let us also consider another point in Figure 6 above (the last of the right point with the yellow cycle). For that point it holds: Y = 8.90, X1 = 1.50 (R), X2 = 158.51 (PET), x3 = 7.59 (Q). According to the results, the left-, and right-hand boundaries and central value are 5.6686 m, 9.8208 m and 7.7447 m, respectively, while $D_j^{GW,obs}$ is equal to 8.90 m. Its membership degree is equal to

$$1-\frac{(8.90-7.7447)}{2.0761}=0.4435$$

It should be noted that in the case that the modified version of Tanaka’s model (FLR-2) is used, the membership degree, which corresponds to the 4th observation, is equal to 0.8210, while the membership degree which correspond to the last observation (9th point) is equal to 0.4400.

In addition, the outputs of the fuzzy regression based on the Tanaka approach (dependent variable) may have different uncertainties expressed by its own membership function, whereas in statistical regression all the outputs have the same error [19]. Last, in fuzzy linear regression there is no theoretical obstacle to take errors in independent variables into account, while fuzzy measure (or measures of vagueness) are commonly used. Therefore, the two approaches should not be interpreted by each other in the same way [19].

5. Concluding Remarks

This study investigates the short-term response of groundwater (in terms of the depth of groundwater) to rainfall, streamflow and the potential evapotranspiration, seeking a fuzzy linear relationship among these hydrometeorological variables. The area under investigation is a shallow unconfined aquifer at the east part of the Nestos River Delta in the Prefecture of Xanthi, Greece where a small size of observed groundwater depth of four wells are related to the hydrometeorological records at the same point in time.

Fuzzy regression analysis is preferred in this application because of the lack of data. Two fuzzy multiple regressions based on Tanaka’s model and a modified version of it are applied. The modified version uses a non-linear objective function which additionally considers the distance between the observed groundwater depth and the central value of the corresponding fuzzy estimate. In both cases, the problem concludes to a constrained optimization problem where all the observed data must be included within the produced fuzzy band.

In fuzzy linear regression based on Tanaka’s approach there is no error term. Suitability of the model is checked through the value of the objective function, while other fuzzy measures may be used. This study uses a fuzzified version of Theil’s inequality coefficient U, which is estimated for each fuzzy regression model. Since Theil’s inequality coefficient U has a fuzzy form, a ranking measure, R, is applied in order to compare Theil’s inequality coefficient with different fuzzy regression models. Another suitability measure, E_dis, takes into account the Euclidian distance between the observed groundwater depth and the corresponding fuzzy estimate as well as the unbiased estimation of the mean value of the historical sample.

Based on suitability measures, the two fuzzy multiple linear regression models performed well in the case of well 194 and thus a fuzzy linear relationship is achieved. The fuzzy regression coefficients show that the groundwater depth is negatively related to rainfall and streamflow and positively related to the potential evapotranspiration. Groundwater seems to be differentially affected by streamflow and rainfall as the coast is approached, while the total fuzziness decreases when the distance of the coast decreases. Correlation which is separately carried out between the dependent variable and the independent ones shows that the groundwater seems to have a stronger linear dependence on streamflow than on rainfall, while the highest correlation coefficient is that of potential evapotranspiration. Meantime, the uncertainties of the influence coefficients of rainfall and streamflow are negligible. However, the contribution of streamflow should not be ignored. Although similar results are obtained from both fuzzy linear regression models, based on the suitability measures, the model which uses the non-linear objective function is preferred.