Geostatistical Characterization of the Transmissivity: An Example of Kuwait Aquifers

: Management of groundwater resources is critical to arid countries like Kuwait. One crucial information that is often not available to groundwater modelers is transmissivity of aquifers. Kriging is used to characterize the transmissivity of the two most commonly used aquifers in Kuwait (the Dammam and Kuwait Group aquifers). The transmissivity of the two aquifers is represented as a random spatial function where heterogeneity is described by the probability distribution and variogram of sample values. A structural analysis was performed which consisted of the construction and interpretation of sample variograms and the selection of model variograms to best ﬁt the structure of the log-transformed transmissivity. The analysis indicated that the transmissivity of the Dammam aquifer is highly anisotropic and has a longer range of inﬂuence and larger variance (39 km; 333,487) than those for the Kuwait Group aquifer (10 km; 44,613). The model variograms were used in the kriging analysis to estimate the spatial average value of the transmissivity of the two aquifers in Kuwait. The estimated transmissivity of the Kuwait Group aquifer ranges from <50 to about 800 m 2 day − 1 , and is controlled principally by the thickness of the saturated zone in the aquifer. The estimated transmissivity in Dammam aquifer ranges from ~100 to 2200 m 2 day − 1 , with high transmissivity ﬁelds coinciding with both the areas of major anticline and fault structures and the areas of aquifer recharge by fresh water, indicating the enhancement of the transmissivity on these structures by dissolution. Low transmissivity ﬁelds are located in the areas where the aquifer waters are stagnant. The performed analysis can be used to aid in the development of numerical models for sustainable management of the aquifer system in Kuwait, deﬁning artiﬁcial recharge sites, and in optimizing the locations of future development wellﬁelds. on broader spatial the effectiveness of kriging in the large-scale extrapolation of transmissivities. This case study highlights, as an example, the two principal aquifers in Kuwait, i.e., Kuwait Group aquifer and Dammam aquifer. Mapping T will aid in the optimum development of groundwater and wellfield management; i.e., understanding the flow system and flow patterns, optimizing placement of development wells, selecting areas for artificial recharge, and conducting numerical modeling studies.


Introduction
Transmissivity (T) is often used in specifying the groundwater extraction potential of an aquifer [1]. Hence, estimating the T is considered an essential parameter for developing local and regional water management plans, modeling groundwater flow and solute transport [2], and selecting areas for additional hydrologic tests [3]. Transmissivity depends both on the properties of the porous matrix and the fluid passing through the strata. A pumping test is the most reliable and standard method to determine the T value in aquifers. However, these pumping tests are expensive, and often, their results are localized spatially. This makes it imperative to extrapolate the information on a larger scale for developing sustainable water management plans. The geostatistical techniques are the most widely used for extrapolation of the measured values to a broader spatial domain [4]. Kriging is a

Geohydrology
Kuwait forms part of the hyper-arid region, with brackish groundwater as the only source of natural potable water. The country obtains its required supply of groundwater from two principal aquifers, i.e., Dammam (Dm) and Kuwait Group (Kg) through established wellfields ( Figure 1). The Dm aquifer is developed in the chalky and dolomitic limestone unit, of the middle part of the Dammam Formation (Middle Eocene), and occurs under confining conditions across most of Kuwait ( Figure 2). The thickness of Dm Formation is from 60 m to 200 m in the northern part of Kuwait. The potentiometric aquifer surface (in the initial stage) varies from about 100 m above sea level in the extreme southwest of Kuwait to a few meters near the eastern coast, with a relatively gentle pre-

Geohydrology
Kuwait forms part of the hyper-arid region, with brackish groundwater as the only source of natural potable water. The country obtains its required supply of groundwater from two principal aquifers, i.e., Dammam (Dm) and Kuwait Group (Kg) through established wellfields ( Figure 1). The Dm aquifer is developed in the chalky and dolomitic limestone unit, of the middle part of the Dammam Formation (Middle Eocene), and occurs under confining conditions across most of Kuwait ( Figure 2). The thickness of Dm Formation is from 60 m to 200 m in the northern part of Kuwait. The potentiometric aquifer surface (in the initial stage) varies from about 100 m above sea level in the extreme southwest of Kuwait to a few meters near the eastern coast, with a relatively gentle pre-development gradient (before 1970) of about 0.8 m km −1 [11]. The direction of groundwater flow is generally towards the east and northeast. Water quality distribution in the Dm aquifer follows its potentiometric trends typically; the total dissolved solids (TDS) increase from~2500 mg L −1 in the southwestern areas of Kuwait to~40,000 mg L −1 near the coast and Kuwait Bay, to >100,000 mg L −1 in northeastern Kuwait. The zone to the northeast, where the TDS values exceed those found in seawater, represents a non-flushed portion of the aquifer filled with connate water [12,13].
The Kg aquifer is developed in the unconsolidated gravel and sand units of the Kuwait Group of Miocene-Pleistocene age and occurs under unconfined conditions throughout Kuwait. Although the upper and lower permeable units of the Kg are separated at places by a semi-pervious sequence of silt, sand, and clay layers (Figure 2), the entire sequence is considered as one hydrographic unit for the aquifer development. The saturated thickness of the aquifer ranges from zero at the southwestern corner of Kuwait, where the aquifer is dry, to over 300 m to the northeast. This general trend in the saturated thickness is interrupted by localized reductions in the coastal areas of the anticlines, and in some areas, complete desaturation (structures at southeast and east of Kuwait). The potentiometric surface of the Kg aquifer follows closely that of the Dm aquifer, and hence, its flow direction. However, the potentiometric values in the former are lower by about 5 m, with head differences increasing towards the coastal area. The Kg aquifer contains water similar in quality to that of the Dm aquifer [14], and ranges from~3000 mg L −1 in the southwest to~130,000 mg L −1 in the northeast over a distance of about 150 km [15].  [11]. The direction of groundwater flow is generally towards the east and northeast. Water quality distribution in the Dm aquifer follows its potentiometric trends typically; the total dissolved solids (TDS) increase from ~2500 mg L −1 in the southwestern areas of Kuwait to ~40,000 mg L −1 near the coast and Kuwait Bay, to >100,000 mg L −1 in northeastern Kuwait. The zone to the northeast, where the TDS values exceed those found in seawater, represents a non-flushed portion of the aquifer filled with connate water [12,13]. The Kg aquifer is developed in the unconsolidated gravel and sand units of the Kuwait Group of Miocene-Pleistocene age and occurs under unconfined conditions throughout Kuwait. Although the upper and lower permeable units of the Kg are separated at places by a semi-pervious sequence of silt, sand, and clay layers (Figure 2), the entire sequence is considered as one hydrographic unit for the aquifer development. The saturated thickness of the aquifer ranges from zero at the southwestern corner of Kuwait, where the aquifer is dry, to over 300 m to the northeast. This general trend in the saturated thickness is interrupted by localized reductions in the coastal areas of the anticlines, and in some areas, complete desaturation (structures at southeast and east of Kuwait). The potentiometric surface of the Kg aquifer follows closely that of the Dm aquifer, and hence, its flow direction. However, the potentiometric values in the former are lower by about 5 m, with head differences increasing towards the coastal area. The Kg aquifer contains water similar in quality to that of the Dm aquifer [14], and ranges from ~3000 mg L −1 in the southwest to ~130,000 mg L −1 in the northeast over a distance of about 150 km [15].  A layer of cherty limestone unit ranging in thickness from 1.5 to 30 m overlying the top of Dm aquifer and clay, averaging 10 m thick, separates the Kg and Dm aquifers ( Figure 2). The entire hydrologic system in Kuwait is bounded on the top by the water table of the Kg aquifer, and on the  base by the relatively impermeable nummulitic limestone with lignites and shale intercalations of the  Dammam Formation, along with the underlying massive anhydrides (20-200 m) of the Rus Formation (Lower Eocene).

Transmissivity Data Sets
Transmissivity data for the Dm and Kg aquifers were obtained from the Ministry of Electricity and Water (MEW), Kuwait. MEW has a record of aquifer tests conducted in wellfields in Kuwait in addition to a number of exploration wells ( Figure 3). The data set consists of 206 measurements for the Dm aquifer and 216 measurements for the Kg aquifer. For both aquifers, transmissivity data used in this study are for those conducted in fully penetrating wells. The transmissivity data are usually reported in imperial gallons per day per foot (Ig day −1 ft −1 ), and were converted to square meter per day (m 2 day −1 ) units.  Kuwait is bounded on the top by the water table of the Kg aquifer, and on the  base by the relatively impermeable nummulitic limestone with lignites and shale intercalations of the  Dammam Formation, along with the underlying massive anhydrides (20-200 m) of the Rus Formation (Lower Eocene).

Transmissivity Data Sets
Transmissivity data for the Dm and Kg aquifers were obtained from the Ministry of Electricity and Water (MEW), Kuwait. MEW has a record of aquifer tests conducted in wellfields in Kuwait in addition to a number of exploration wells ( Figure 3). The data set consists of 206 measurements for the Dm aquifer and 216 measurements for the Kg aquifer. For both aquifers, transmissivity data used in this study are for those conducted in fully penetrating wells. The transmissivity data are usually reported in imperial gallons per day per foot (Ig day −1 ft −1 ), and were converted to square meter per day (m 2 day −1 ) units. The two units of Kg aquifer (upper and lower) are considered as one hydrostratigraphic unit for the aquifer development, and most of the pumping tests conducted in the aquifer have been carried out with screened intervals covering both the units, and T is an average of the two units. The two units of Kg aquifer (upper and lower) are considered as one hydrostratigraphic unit for the aquifer development, and most of the pumping tests conducted in the aquifer have been carried out with screened intervals covering both the units, and T is an average of the two units.

Statistical Analysis
The T data were statistically analyzed, after cleaning the data of the outliers, following the procedure suggested by Woodbury and Sudicky [17]. The three outliers from Kg aquifer and four from Dm aquifer where the T values were <20 m 2 day −1 were not considered in statistical analyses. A single value of >4500 m 2 day −1 for Dm aquifer was also omitted from statistical analyses. This approach of Woodbury and Sudicky seems to provide more realistic results in both determining population statistics and variogram estimation. Table 1 presents a comparison of the corresponding statistics for the T in the two aquifers. T of Dm aquifer has a higher mean and higher variance compared to the Kg aquifer. The higher variation in T of the Dm aquifer is indicative of the karst within this aquifer. The histogram of Kg and Dm aquifers appears to be moderately skewed to the right (Figure 4a,b), as indicated by the values of the coefficient of skewness. The two histograms indicate marked asymmetry, especially the Dm aquifer, with tails towards the higher values of transmissivity, suggesting log-normal distribution. This behavior for transmissivity histograms appears to be consistent with those observed in geostatistical transmissivity studies elsewhere [5,[17][18][19][20].

Statistical Analysis
The T data were statistically analyzed, after cleaning the data of the outliers, following the procedure suggested by Woodbury and Sudicky [17]. The three outliers from Kg aquifer and four from Dm aquifer where the T values were <20 m 2 day −1 were not considered in statistical analyses. A single value of >4500 m 2 day −1 for Dm aquifer was also omitted from statistical analyses. This approach of Woodbury and Sudicky seems to provide more realistic results in both determining population statistics and variogram estimation. Table 1 presents a comparison of the corresponding statistics for the T in the two aquifers. T of Dm aquifer has a higher mean and higher variance compared to the Kg aquifer. The higher variation in T of the Dm aquifer is indicative of the karst within this aquifer. The histogram of Kg and Dm aquifers appears to be moderately skewed to the right ( Figure  4a,b), as indicated by the values of the coefficient of skewness. The two histograms indicate marked asymmetry, especially the Dm aquifer, with tails towards the higher values of transmissivity, suggesting log-normal distribution. This behavior for transmissivity histograms appears to be consistent with those observed in geostatistical transmissivity studies elsewhere [5,[17][18][19][20]. In this study, the logarithm of the data was considered for geostastical analyses, the T for any unknown location was then extrapolated by merely taking the exponential of the estimated value; this approach was successfully used in earlier studies [19,21]. The variogram shows a better correlation with the measured values [22]. The approach was tested for both Kg and Dm aquifers (Figure 5a,b). The percentage cumulative frequency distribution shows a reasonable correlation with a straight line (normality condition), in spite of the fact that this data was skewed, indicating validity of the use of log transformed T.  In this study, the logarithm of the data was considered for geostastical analyses, the T for any unknown location was then extrapolated by merely taking the exponential of the estimated value; this approach was successfully used in earlier studies [19,21]. The variogram shows a better correlation with the measured values [22]. The approach was tested for both Kg and Dm aquifers (Figure 5a,b).
The percentage cumulative frequency distribution shows a reasonable correlation with a straight line (normality condition), in spite of the fact that this data was skewed, indicating validity of the use of log transformed T.

Structural Analysis
The estimation of the sample semivariogram and selection of an appropriate theoretical semivariogram is the fundamental basis of any geostatistical analysis [23]. Semivariograms describe the spatial correlation structure of the regionalized variables, which is defined as where γ(h) is the variogram function, and Z(x) and Z(x + h) are two of the random variables comprising the random field. The random field is assumed to satisfy the intrinsic hypothesis; i.e., the is second-order stationary, meaning that its mean function is constant, and its covariance function depends only on the separation distance h (Olea, 1984; Philip and Kitanidis, 1989). The experimental variogram and the theoretical model of the log(T) of both aquifers were computed by using the GS+ Geostatistics for the Environmental Sciences version 7.0, 2004.

Experimental Variogram
There are numerous methods that are used for computation of the experimental variogram including: the "classical" estimator [24], the Cressie-Hawkins estimator [25], the squared median of the absolute deviation (SMAD) estimator [26], and Omre estimator [27]. In this study, we have used the classical estimator for calculating the experimental variogram after screening outliers from the database, as suggested by Woodbury and Sudicky [17].
The classical estimator equation is [22] (ℎ) where n(h) is the number of data pairs separated by the vector, or lag distance, h, and Z(xi) are measured values of log transmissivity at coordinates xi. The estimator γ(h) is the experimental variogram. The classical estimator is an optimal estimator of the variogram, provided that the pairs Z(xi) and Z(xi + h) are bivariate normal and independent. The experimental semivariogram for specified distances and direction were obtained for both the Kg and Dm aquifers (Figure 6a,b) utilizing an approach proposed by Istok et al. [28]. The

Structural Analysis
The estimation of the sample semivariogram and selection of an appropriate theoretical semivariogram is the fundamental basis of any geostatistical analysis [23]. Semivariograms describe the spatial correlation structure of the regionalized variables, which is defined as where γ(h) is the variogram function, and Z(x) and Z(x + h) are two of the random variables comprising the random field. The random field is assumed to satisfy the intrinsic hypothesis; i.e., the difference is second-order stationary, meaning that its mean function is constant, and its covariance function depends only on the separation distance h (Olea, 1984; Philip and Kitanidis, 1989). The experimental variogram and the theoretical model of the log(T) of both aquifers were computed by using the GS+ Geostatistics for the Environmental Sciences version 7.0, 2004.

Experimental Variogram
There are numerous methods that are used for computation of the experimental variogram including: the "classical" estimator [24], the Cressie-Hawkins estimator [25], the squared median of the absolute deviation (SMAD) estimator [26], and Omre estimator [27]. In this study, we have used the classical estimator for calculating the experimental variogram after screening outliers from the database, as suggested by Woodbury and Sudicky [17].
The classical estimator equation is [22] where n(h) is the number of data pairs separated by the vector, or lag distance, h, and Z(x i ) are measured values of log transmissivity at coordinates x i . The estimator γ(h) is the experimental variogram. The classical estimator is an optimal estimator of the variogram, provided that the pairs Z(x i ) and Z(x i + h) are bivariate normal and independent. The experimental semivariogram for specified distances and direction were obtained for both the Kg and Dm aquifers (Figure 6a,b) utilizing an approach proposed by Istok et al. [28]. The experimental variogram for the two aquifers shows a clear sill, and the presence of spatial correlation at the shortest lag.
experimental variogram for the two aquifers shows a clear sill, and the presence of spatial correlation at the shortest lag.

Model Semivariogram
The exponential model was used to fit the experimental variogram for both Kg and Dm aquifers. An iterative process was applied to determine the optimum length of the variogram, by reducing the lag distance, which is the range of the maximum and minimum distance between points over which the experimental variogram is calculated. The lag distance was reduced from 50.0 km to 25.0 for Kg, and from 90.0 km to 45.0 km for Dm aquifer. A regular lag class distance interval was also systematically reduced, from 6.0 km to 1.5 km, with each active lag distance for both aquifers. The exponential models were found to be more applicable for both aquifers. The simple exponential model equation used was as follows [29]: where C0 is the nugget variance, and C0 + C the sill. C0 and C represent the random and spatial components of variation, respectively, a0 is the model's range parameter, which is distinguished form a, the range of influence, or the separation distance over which spatial dependence is apparent, for the exponential model a, equals 3a0. The residual sums of squares (RSS), r 2 , and proportion C/(C0 + C) were used to select the fitted model. The proportion C/(C0 + C), is the proportion of the spatially structured variance (C) to the sample variance (C0 + C). This value will be 1.0 for a variogram with no nugget variance (where the curve passes through the origin); conversely, it will be 0.0, where there is no spatially dependent variation at the range specified, i.e., where there is a pure nugget effect. Table 2 lists the fitted parameters to the experimental variogram, while Figure 6a,b displays the fitted model variogram of log(T) for the Kg and Dm aquifers, respectively. It can be seen from the fitted variogram models that the spatial characteristics of two aquifers are significantly different; the Dm aquifer has a larger range (i.e., the distance beyond which log(T) is uncorrelated) and sill values compared to Kg aquifer. The range of influence is about 39 km for Dm aquifer and 10 km for Kg aquifer, respectively. These results are in general agreement with observations made for variograms on the larger spatial scale; ranges can lie between less than 1.0-10.0 km for alluvial aquifers, and can reach up to 20.0 km for limestone and chalk aquifers [19]. This behavior can be attributed to the different nature of transmissivity in the two aquifers. The Kg aquifer is characterized by lithology variation within a short distance that affects its permeability, and transmissivity. While Dm aquifer, has a uniform carbonate lithology, with significant structural elements and karst development, the spatial correlation of the developed transmissivity is expected to be larger.

Model Semivariogram
The exponential model was used to fit the experimental variogram for both Kg and Dm aquifers. An iterative process was applied to determine the optimum length of the variogram, by reducing the lag distance, which is the range of the maximum and minimum distance between points over which the experimental variogram is calculated. The lag distance was reduced from 50.0 km to 25.0 for Kg, and from 90.0 km to 45.0 km for Dm aquifer. A regular lag class distance interval was also systematically reduced, from 6.0 km to 1.5 km, with each active lag distance for both aquifers. The exponential models were found to be more applicable for both aquifers. The simple exponential model equation used was as follows [29]: where C 0 is the nugget variance, and C 0 + C the sill. C 0 and C represent the random and spatial components of variation, respectively, a 0 is the model's range parameter, which is distinguished form a, the range of influence, or the separation distance over which spatial dependence is apparent, for the exponential model a, equals 3a 0 . The residual sums of squares (RSS), r 2 , and proportion C/(C 0 + C) were used to select the fitted model. The proportion C/(C 0 + C), is the proportion of the spatially structured variance (C) to the sample variance (C 0 + C). This value will be 1.0 for a variogram with no nugget variance (where the curve passes through the origin); conversely, it will be 0.0, where there is no spatially dependent variation at the range specified, i.e., where there is a pure nugget effect. Table 2 lists the fitted parameters to the experimental variogram, while Figure 6a,b displays the fitted model variogram of log(T) for the Kg and Dm aquifers, respectively. It can be seen from the fitted variogram models that the spatial characteristics of two aquifers are significantly different; the Dm aquifer has a larger range (i.e., the distance beyond which log(T) is uncorrelated) and sill values compared to Kg aquifer. The range of influence is about 39 km for Dm aquifer and 10 km for Kg aquifer, respectively. These results are in general agreement with observations made for variograms on the larger spatial scale; ranges can lie between less than 1.0-10.0 km for alluvial aquifers, and can reach up to 20.0 km for limestone and chalk aquifers [19]. This behavior can be attributed to the different nature of transmissivity in the two aquifers. The Kg aquifer is characterized by lithology variation within a short distance that affects its permeability, and transmissivity. While Dm aquifer, has a uniform carbonate lithology, with significant structural elements and karst development, the spatial correlation of the developed transmissivity is expected to be larger.

Directional Variograms (Anisotropy Analysis)
The directional semivariogram for four directional classes (0 • , 45 • , 90 • , and 135 • ) was computed for both the Kg and Dm aquifers. There Kg aquifer shows anisotropy in T values in 90 • (Figure 7). However, in Dm aquifer, a very clear anisotropy is evident (Figure 8), and there are distinct differences in the structure in all the four directions; the two directional variograms (representing 0 • and 45 • ) rise more sharply than the other two directions (90 • and 135 • ). This confirms the presumption of earlier workers that the Dm aquifer possesses highly anisotropic transmissivity due to karst networks [12]. Moreover, the dominantly observed orientation of the anisotropy becomes more significant when related to the structural geology of the Dammam Formation. The general trend of the major elongated anticline and synclinal structures in southern and central Kuwait (Figure 9) is in the north-northwest direction [30]. These structures are associated with major faulting, fracturing, shear zones as well as very small closure anticlines, running in the north-northeast direction [31,32]. This corresponds to the T from the directional variogram.

Directional Variograms (Anisotropy Analysis)
The directional semivariogram for four directional classes (0°, 45°, 90°, and 135°) was computed for both the Kg and Dm aquifers. There Kg aquifer shows anisotropy in T values in 90° (Figure 7). However, in Dm aquifer, a very clear anisotropy is evident (Figure 8), and there are distinct differences in the structure in all the four directions; the two directional variograms (representing 0° and 45°) rise more sharply than the other two directions (90° and 135°). This confirms the presumption of earlier workers that the Dm aquifer possesses highly anisotropic transmissivity due to karst networks [12]. Moreover, the dominantly observed orientation of the anisotropy becomes more significant when related to the structural geology of the Dammam Formation. The general trend of the major elongated anticline and synclinal structures in southern and central Kuwait (Figure 9) is in the north-northwest direction [30]. These structures are associated with major faulting, fracturing, shear zones as well as very small closure anticlines, running in the north-northeast direction [31,32]. This corresponds to the T from the directional variogram.

Transmissivity Kriging
Kriging is an effective interpolation technique that considers both distance and degree of variation between data points for estimating value at any unknown location [7,34]. The kriging weights were obtained by fitting semi-variogram models, in this study block ordinary kriging was used. The linear, Gaussian, spherical, and exponential models were tested with the lag distance to select the best fitting.
A block ordinary kriging using the 16 nearest neighbors within a radius of 167 and 170 km for Kg and Dm aquifer, respectively, was performed to estimate the transmissivity values at an unknown location in both Kg and Dm aquifers. The log(T) structural models, given by Equation (4), were used in the estimation of the transmissivity function distribution. Furthermore, the variance of the error of the estimate was obtained. The computation was carried out using the three-dimensional kriging module of the GS+ software, which uses the following kriging equations [18]: where λ j 0 is the weighting function of each point, j, to be used in the estimation of point x 0 ; γ(x i − x j ) is the variogram value for the distance between x i and x j , and µ is the Lagrange multiplier. The variance of the error of estimation for point x 0 is given as

Cross Validation
Cross validation was carried out to reexamine the fitted models by systematically changing the lag distance, uniform distance, and the model type. Four types of variogram models, including linear, Gaussian, spherical and exponential, were reexamined with different lag distance. Scatter plots of the true values versus the estimated values are illustrated in Figure 10a,b. These show regression coefficients of 0.935 and 0.901, for Kg and Dm aquifers, respectively.
The coefficient of the least squares model describes the linear regression equation between the actual and estimated values. The regression coefficient will be 1.0 for the best fit, shown by the blue line in Figure 10a,b. These values were the highest values obtained among the examined models and the lag distances. Accordingly, the exponential model was selected based on the criteria and the difference between the actual and estimated values (Table 3).  The r 2 is low for both Kg (0.365) and Dm (0.418). It is worth mentioning that high r 2 is certainly a yardstick for the fitness of the model, but not the sole criterion, especially for predictive purposes. The low root mean square error (RMSE) is of more practical usefulness in a predictive model [35]. A useful rule of thumb is that the test set RMSE should be less than 10% of the range of target property value. The RMSE was calculated for both the Kg and Dm aquifers ( Figure 11 and Table 3). The mean of the error of Kg and Dm were 2.309 and −14.930, respectively. The means are not zero, although they are well within the 10% range of the target value.
The lower RMSE for Kg aquifer can be explained with a lesser variation in T values i.e., 61 to 1252 m 2 d −1 compared to Dm aquifer, where the RMSE was −14.93 and the T value ranged from 22 to 2237 m 2 d −1 . This higher variation in T within Dm aquifer is due to karst development, and the fractures within the carbonate aquifer; field observations have also confirmed these inferences and frequent loss of circulation was observed in drilling reports. The T values for Dm aquifers are underestimated, but within very close range of the measured values. The r 2 obtained after cross validation for the simulated values for Kg and Dm aquifers were 0.418 and 0.423, and they were within the range of 0.366 to 0.428 for the different active lag distances, class distance, and models, and the standard error prediction defined by relation SDx (1 − r 2 ), where the SD = standard deviation of the actual data. The r 2 is low for both Kg (0.365) and Dm (0.418). It is worth mentioning that high r 2 is certainly a yardstick for the fitness of the model, but not the sole criterion, especially for predictive purposes. The low root mean square error (RMSE) is of more practical usefulness in a predictive model [35]. A useful rule of thumb is that the test set RMSE should be less than 10% of the range of target property value. The RMSE was calculated for both the Kg and Dm aquifers ( Figure 11 and Table 3). The mean of the error of Kg and Dm were 2.309 and −14.930, respectively. The means are not zero, although they are well within the 10% range of the target value.
The lower RMSE for Kg aquifer can be explained with a lesser variation in T values i.e., 61 to 1252 m 2 d −1 compared to Dm aquifer, where the RMSE was −14.93 and the T value ranged from 22 to 2237 m 2 d −1 . This higher variation in T within Dm aquifer is due to karst development, and the fractures within the carbonate aquifer; field observations have also confirmed these inferences and frequent loss of circulation was observed in drilling reports. The T values for Dm aquifers are underestimated, but within very close range of the measured values. The r 2 obtained after cross validation for the simulated values for Kg and Dm aquifers were 0.418 and 0.423, and they were within the range of 0.366 to 0.428 for the different active lag distances, class distance, and models, and the standard error prediction defined by relation SDx(1 − r 2 ), where the SD = standard deviation of the actual data.

Results
A spatial map has been generated from the T values obtained by geostatistical techniques for Kg and Dm aquifers (Figures 12 and 13). The transmissivity distribution map of Kg aquifer ( Figure 12) indicates that the transmissivity of the aquifer increases from <50 m 2 day −1 in the southwestern corner of Kuwait to about 400 m 2 day −1 to the northeast, which appears to be following the general trend in the saturated thickness of the aquifer. This general trend is interrupted by several local anomalies where the aquifer has low transmissivity due to presence of a gatch layer. The low transmissivity anomalies (<50 m 2 day −1 ) are attributed to the reduction of the aquifer-saturated thickness on top of the anticline structures ( Figure 3). The high transmissivity anomalies (>700 m 2 day −1 ) are attributed to the thickening of the saturated aquifer at the synclinal structure, and may be also due to the presence of extended faults from deeper strata.

Results
A spatial map has been generated from the T values obtained by geostatistical techniques for Kg and Dm aquifers (Figures 12 and 13). The transmissivity distribution map of Kg aquifer ( Figure 12) indicates that the transmissivity of the aquifer increases from <50 m 2 day −1 in the southwestern corner of Kuwait to about 400 m 2 day −1 to the northeast, which appears to be following the general trend in the saturated thickness of the aquifer. This general trend is interrupted by several local anomalies where the aquifer has low transmissivity due to presence of a gatch layer. The low transmissivity anomalies (<50 m 2 day −1 ) are attributed to the reduction of the aquifer-saturated thickness on top of the anticline structures ( Figure 3). The high transmissivity anomalies (>700 m 2 day −1 ) are attributed to the thickening of the saturated aquifer at the synclinal structure, and may be also due to the presence of extended faults from deeper strata. The spatial map for T of Dm aquifer depicts a ranges <100 m 2 day −1 to ~2200 m 2 day −1 (Figure 13). The high T values are located in the southern parts of Kuwait, and coincide with the areas of major anticline, lineations, and fault. Furthermore, these high transmissivity fields are located in the areas where the aquifer is recharged by fresher water from the southwest. On the other hand, low T fields are found to the north and northeast of these structures, and where the water in Dm aquifer is stagnant. This might be used to confirm that the Dm aquifer transmissivity is enhanced by preferential dissolution of the original rock matrix along weak planes (such as faults, joints, etc.), which are concentrated at the areas of anticlines. The spatial data coverage of Dm aquifer was much larger than Kg aquifer (Figure 3). The T obtained from kriging for the two aquifers has a certain degree of correlation. Co-kriging was the next logical step in this geostatistical analysis, with the aim that T of Kg aquifer could be better The spatial map for T of Dm aquifer depicts a ranges <100 m 2 day −1 to~2200 m 2 day −1 (Figure 13). The high T values are located in the southern parts of Kuwait, and coincide with the areas of major anticline, lineations, and fault. Furthermore, these high transmissivity fields are located in the areas where the aquifer is recharged by fresher water from the southwest. On the other hand, low T fields are found to the north and northeast of these structures, and where the water in Dm aquifer is stagnant. This might be used to confirm that the Dm aquifer transmissivity is enhanced by preferential dissolution of the original rock matrix along weak planes (such as faults, joints, etc.), which are concentrated at the areas of anticlines. The spatial map for T of Dm aquifer depicts a ranges <100 m 2 day −1 to ~2200 m 2 day −1 (Figure 13). The high T values are located in the southern parts of Kuwait, and coincide with the areas of major anticline, lineations, and fault. Furthermore, these high transmissivity fields are located in the areas where the aquifer is recharged by fresher water from the southwest. On the other hand, low T fields are found to the north and northeast of these structures, and where the water in Dm aquifer is stagnant. This might be used to confirm that the Dm aquifer transmissivity is enhanced by preferential dissolution of the original rock matrix along weak planes (such as faults, joints, etc.), which are concentrated at the areas of anticlines. The spatial data coverage of Dm aquifer was much larger than Kg aquifer (Figure 3). The T obtained from kriging for the two aquifers has a certain degree of correlation. Co-kriging was the next logical step in this geostatistical analysis, with the aim that T of Kg aquifer could be better The spatial data coverage of Dm aquifer was much larger than Kg aquifer (Figure 3). The T obtained from kriging for the two aquifers has a certain degree of correlation. Co-kriging was the next logical step in this geostatistical analysis, with the aim that T of Kg aquifer could be better estimated using the widespread Dm aquifer T. There were 78 identical locations that were tested. The T values for these were not correlated, and the correlation coefficient was 0.047 ( Figure 14). This is essentially due to the difference in intrinsic characteristics of the two aquifers compared.
Water 2018, 10, x FOR PEER REVIEW 14 of 16 estimated using the widespread Dm aquifer T. There were 78 identical locations that were tested. The T values for these were not correlated, and the correlation coefficient was 0.047 ( Figure 14). This is essentially due to the difference in intrinsic characteristics of the two aquifers compared.

Discussion
The applied geostatistical techniques are useful for extrapolating the T over large spatial domains. The use of kriging for T estimation proved to be valid both for Kg aquifer that contains sand, silt, and clay (essentially a variable lithology), and the Dm aquifer having nearly uniform limestone lithology. The use of geostatistical tools to derive T values correlate well with the measured values, and can reduce the uncertainty in wellfield development, and can complement the expensive well testing.
The T of the two aquifers is represented as a random spatial function where heterogeneity is described by the probability distribution and variogram of sample values. A structural analysis was performed which consisted of the construction and interpretation of sample variograms and the selection of model variograms to best fit the structure of the log-transformed transmissivity. The analysis indicated that the Dm aquifer is highly anisotropic, and has a longer range of influence and more substantial variance (39 km, 333,487) than the Kg aquifer (10 km, 44,613).
The model variograms were used in the kriging analysis to estimate the spatial average value of the T for the two aquifers in Kuwait. The estimated transmissivity distribution in the two aquifers was validated by geological and hydrological data. The estimated Dm aquifer T ranges from about 100 to 2200 m 2 /day, with high transmissivity fields coinciding with both the areas of major anticline and fault structures, and the areas of aquifer recharge by fresh water. Low T fields are located in the areas where the aquifer waters are stagnant. The estimated T for Kg aquifer ranges from <50 to about 800 m 2 day −1 , and is controlled principally by the saturated thickness of the aquifer.
The performed structural analysis and the estimated transmissivity and their associated variance fields for the two aquifers can be a critical input for developing a groundwater model that can aid in sustainable development and management of the wellfields in Kuwait.
The exponential model and lag distance are useful to estimate T values in a similar geological setting, both in isotropic and anisotropic formations. The extrapolated T values using kriging can be used for development of numerical modeling studies (deterministic and stochastic) for the aquifer system in the region. The T values can be used effectively in defining artificial recharge or underground water storage areas in the aquifers and in optimizing placement of future wellfield developments.

Discussion
The applied geostatistical techniques are useful for extrapolating the T over large spatial domains. The use of kriging for T estimation proved to be valid both for Kg aquifer that contains sand, silt, and clay (essentially a variable lithology), and the Dm aquifer having nearly uniform limestone lithology. The use of geostatistical tools to derive T values correlate well with the measured values, and can reduce the uncertainty in wellfield development, and can complement the expensive well testing.
The T of the two aquifers is represented as a random spatial function where heterogeneity is described by the probability distribution and variogram of sample values. A structural analysis was performed which consisted of the construction and interpretation of sample variograms and the selection of model variograms to best fit the structure of the log-transformed transmissivity. The analysis indicated that the Dm aquifer is highly anisotropic, and has a longer range of influence and more substantial variance (39 km, 333,487) than the Kg aquifer (10 km, 44,613).
The model variograms were used in the kriging analysis to estimate the spatial average value of the T for the two aquifers in Kuwait. The estimated transmissivity distribution in the two aquifers was validated by geological and hydrological data. The estimated Dm aquifer T ranges from about 100 to 2200 m 2 /day, with high transmissivity fields coinciding with both the areas of major anticline and fault structures, and the areas of aquifer recharge by fresh water. Low T fields are located in the areas where the aquifer waters are stagnant. The estimated T for Kg aquifer ranges from <50 to about 800 m 2 day −1 , and is controlled principally by the saturated thickness of the aquifer.
The performed structural analysis and the estimated transmissivity and their associated variance fields for the two aquifers can be a critical input for developing a groundwater model that can aid in sustainable development and management of the wellfields in Kuwait.
The exponential model and lag distance are useful to estimate T values in a similar geological setting, both in isotropic and anisotropic formations. The extrapolated T values using kriging can be used for development of numerical modeling studies (deterministic and stochastic) for the aquifer system in the region. The T values can be used effectively in defining artificial recharge or underground water storage areas in the aquifers and in optimizing placement of future wellfield developments.