The Irrigation Cooling Effect as a Climate Regulation Service of Agroecosystems

Agroecosystems provide a range of benefits to society and the economy, which we call ecosystem services (ES). These services can be evaluated on the basis of environmental and socioeconomic indicators. The irrigation cooling effect (ICE), given its influence on the land surface temperature (LST), is an indicator of climate regulation services from agroecosystems. In this context, the objective of this study is to quantify the ICE in agroecosystems at the local scale. The agroecosystem of citrus cultivation in Campo de Cartagena (Murcia, Spain) is used as a case study. Once the LST was retrieved by remote sensing images for 216 plots, multivariate regression methods were used to identify the factors that explain ICE. The use of a geographically weighted regression (GWR) model is proposed, instead of ordinary least squares, as it offsets the spatial dependence and gives a better fit. The GWR explains 78% of the variability in the LST, by means of three variables: the vegetation index, the water index of the crop, and the altitude. Thus, the effects of the change in land use on the LST due to restrictions on the availability of water (up to 1.22 °C higher for rain-fed crops) are estimated. The trade-offs between ICE and the other ES are investigated by using the irrigation water required to reduce the temperature. This work shows the magnitude of the climate regulation service generated by irrigated citrus and enables its quantification in agroecosystems with similar characteristics.


Introduction
The main function of agroecosystems is traditional food and fiber production. However, other goods and services are offered by these systems, which are called ecosystem services (ES) [1]. These services can be defined as the contributions of an ecosystem to human wellbeing, which are generally classified into provisioning, support, regulation, and cultural services [2]. ES provide a framework to identify and analyze all agricultural outputs [1]. Nevertheless, not all ES contributions are positive, as there are also the so-called ecosystem disservices (EDS). EDS can be defined as the ecosystem-generated functions, processes, and attributes which result in perceived or actual negative impacts on human wellbeing [3]. In agroecosystems, there are often trade-offs between ES and EDS. For example, an increase in food production requires a greater demand for water [4]. These trade-offs become even more important in semi-arid environments, where there is water scarcity affecting the ES supply of agroecosystems [5]. Therefore, it is necessary to know the ES and EDS derived from the irrigation of semi-arid agroecosystems [6,7]. Thus, water consumption will be an EDS in these areas; however, at the same time, it will produce an ES by allowing, among others, the regulation of temperature extremes [8].
between the LST and its influencing factors in the context of the urban heat island [40][41][42][43][44] but not in the generation of the ICE.
In this context, this paper proposes the quantification of the ICE in agroecosystems at the local scale, through the study of spatial patterns of LST reduction produced in irrigated agriculture. OLS and GWR models are used to test which method allows for its correct quantification [40][41][42][43][44].
In addition to the existence of a cooling effect in irrigated agroecosystems in the study area, the factors determining the ICE are studied [36]. To this end, multivariate regression techniques are used, paying special attention to the character of the sample. Modelling makes it possible to quantify the effect of changes in land-use and land-cover on this climate phenomenon and, in short, on the ES of climate regulation. This quantification is important to obtain temperature data that allows the subsequent joint valuation of ES in these agroecosystems [16,17].
In order to achieve this objective, irrigated citrus fruit crops in an agricultural region of the Western Mediterranean with limited water resources are studied. In the last 40 years, the area has undergone a profound change in land use, moving from dry to irrigated systems, especially in the case of citrus. This has also happened in other areas of the Mediterranean basin [45], Florida [46], and California [47], which makes this case study of interest beyond the local scale.

Study Area
The agricultural region of Campo de Cartagena, located in the Region of Murcia (Southeastern Spain; see Figure 1), is affected by a semi-arid climate, with an average annual temperature and rainfall of 18 • C and 350 mm, respectively. The use of water resources from the Tagus-Segura transfer (122 hm 3 /year) and groundwater (89.44 hm 3 /year) [48] has led to this region undergoing, in recent decades, a very important change in both land-use and current production orientations [49]. Thus, this region increased its irrigated agricultural area by 57.13% over the period of 1984-2016. Currently, there is a total of 47,430 ha of irrigated land, while rain-fed land has declined to 13,896 ha [50]. Of the irrigated land in the study area, there is a notable extension (8000 ha) dedicated to citrus cultivation (lemon, orange, mandarin, and grapefruit).

Methodology
In order to quantify the ICE through the study of spatial LST patterns, information was gathered from 216 citrus plots registered in Campo de Cartagena in 2016 from the Information System on Land Occupancy in Spain (SIOSE) [51]. These plots have a drip irrigation system and an average area of 3.10 ha. All plots had the same average water supply managed by Irrigation Communities. For each

Methodology
In order to quantify the ICE through the study of spatial LST patterns, information was gathered from 216 citrus plots registered in Campo de Cartagena in 2016 from the Information System on Land Occupancy in Spain (SIOSE) [51]. These plots have a drip irrigation system and an average area of 3.10 ha. All plots had the same average water supply managed by Irrigation Communities. For each plot, both the LST and the different factors influencing it were measured. The factors which could explain the variability in LST were chosen through a review of previous related work [33,52]: air temperature (Tair), Normalized Difference Vegetation Index (NDVI), Normalized Difference Water Index (NDWI), distance to the coast (Dist-coast), altitude, slope, and dominant orientation. The techniques used for factor estimations were remote sensing and the use of spatial databases with GIS tools.
This work complements the study of Zabala et al. [8], which attempted to identify all the ES of Mediterranean irrigated agriculture. One of the ES identified as relevant was climate regulation. Considering water consumption as a limiting factor, changes in the supply of this ES can be observed. These changes will occur when moving from an irrigated to a rain-fed agroecosystem, as a result of water scarcity in the area. However, climate regulation ES was not quantified previously, this work being the first to obtain climate regulation ES by retrieving LST from satellite images. In fact, local climate regulation ES is the ES least commonly quantified in agroecosystems. In addition, water consumption was also taken into account, with an analysis of the trade-offs with other ES.

Land Surface Temperature retrieval and its influencing factors
The LST ( • C) was retrieved from the information contained in satellite images of the Landsat 8-OLI-TIRS remote sensor of the United Stated Geological Survey (USGS). The image under consideration has multispectral bands with a spatial resolution of 30 m (in 100 m thermal bands) and the coordinate system ETRS89. The image was captured under clear conditions (0% cloudy) in summer (27 July 2016 at 10:44:01 a.m.) and was corrected radiometrically, passing values from Digital Numbers to top-of-atmosphere values. A single-satellite image was selected, in summer dates, as is usual in other studies [14,26,44,[53][54][55], due to the great impact that irrigation has on extreme temperatures and the scarce influence on other seasons of the year. Thus, the use of a single-satellite image was intended to carry out a spatial comparison between land uses, rather than a temporal analysis.
For the retrieval of the LST, the Chander and Markham technique [56] was used. This technique was proposed by the USGS [57], does not require climatic variables, and is generally used to retrieve this variable from satellite information using the expression: where T B is the brightness temperature, W is the wavelength of emitted radiation, P is a constant (14,380), and LSE is the emissivity correction factor proposed by Sobrino [58]. The first two variables are provided by the image used.
To explain the LST factors, information was obtained from the variables NDVI and NDWI, obtained from the same satellite image: where NIR is near-infrared (band 5), R is visible red (band 4), and SWIR is short wavelength-infrared (band 6). The combination of these bands allowed us to measure the intensity of the green vegetation and the amount of water it contains [59]. The values of these indices range from −1 to 1; the higher the value, the healthier and denser the vegetation. Differences in values are determined by the irrigation of the area, after confirming the absence of rainfall days prior to the satellite image selected. For the previous month, no rain events (0 mm) occurred in the study area [60]. Average Tair in • C was based on the interpolation of Tair values from meteorological stations in the study area [60].
Average altitude (m), average slope ( • ), and orientations (−1: shaded spot, 0: flat, 1: sunny spot) were obtained from the digital elevation model [61]. Given the different shapes of the plots, the averages of these plots were selected to obtain the average altitudes and slopes. In the case of the orientations, the most frequent orientation value of the plot was selected.
Average Dist-coast of the irrigated citrus plots (km) [51]; the study region is adjacent to the Mediterranean Sea, which is largely affected by the sea currents and wind. For this study, Dist-coast was considered as a proxy for the oceanic and atmospheric states [43].
The information used in this study, as described in Table A1, is fully recorded in a data repository [62], where data from the variables used in the multiple linear regression models for the 216 citrus plots are shown.

Spatial Regression Model
The OLS method is the most frequently used method to estimate a linear regression model. This is due to its simplicity and the optimal character of the estimated parameters for cross-sectional data sets. This technique has been used to study samples distributed in space, assuming that the relationships are spatially constant [41]. A global regression model is expressed as: where, for each observation i, y i is the dependent variable, x i1 to x in are the independent variables, β 0 is the intercept, β 1 to β n are the estimated coefficients, and ε i is the sampling error, these parameters being invariant. Nevertheless, OLS estimation with spatial data often fails to comply with one of the basic hypotheses of this technique, such as the independence of residuals. The residuals are spatially autocorrelated due to the spatial dependence of the study variable. This leads to OLS estimators being inefficient, and hence, the model suffers from poor specification.
To determine the degree of spatial concentration or dispersion of the data, an indicator of spatial autocorrelation is applied, such as Moran's I index [63,64]. In the presence of spatial autocorrelation, it is necessary to apply a technique specifically designed to treat and model it, such as GWR [38][39][40][41][42].
The GWR technique makes it possible to estimate multivariate regression models by incorporating the spatial dependence of the data. The GWR model makes the estimates of coefficients a function of the location and quantifies them separately and independently for each sample unit, passing from a global analysis to a local one [63]. The GWR model can be expressed as: where (µ i , v i ) indicates the co-ordinates of each observation (which act as weights), and β(µ i , v i ) are the coefficients of the model (which vary with the location). Thus, for each observation, the estimate provides both a value for each coefficient and a measure of the adjustment, allowing the local variation and the effects of the variables to be studied individually. A weighting scheme known as adaptive Gaussian Kernel [64] was used in the allocation of GWR model weights. It assigns a higher weight to the closest observations, which is zero if the distance between the two observations exceeds the bandwidth distance, calculated for each case by cross-validation [62].
In order to compare the quality of the models, the corrected Akaike information criterion (AICc) and the corrected determination coefficient (R 2 c) [37,41,64] were used. In addition, for prediction purposes, the root mean square error (RMSE) was used to compare the measure of performance of the final models [65].

Spatial Distribution of the LST and its Influencing Factors
Firstly, a satellite image was selected from the Landsat 8 remote sensor on 27 July 2016. The image corresponds to the period in which the highest temperatures of the year occurred in the study area [48]. From this satellite image, the LST was retrieved. Figure 2a shows the spatial distribution of the LST, with maxima in the central-western zone and minima in the east and south. This distribution is consistent with that of the agricultural uses by comparing LST from irrigated and dry lands (Figure 2b). On 27 July 2016, more than one-third of the study area was above 35 • C, while the 216 irrigated citrus plots in the agricultural zone ( Figure 2c) had an average LST of 32.47 • C.

LST Models
Based on the factors selected to explain the ICE, a multiple linear regression model was first fitted using an OLS procedure. Using the backward stepwise selection procedure, the model was found to have three significant factors: NDVI, NDWI, and the distance of the plots to the coast (Distcoast). This specification explained practically three-quarters of the LST variability in irrigated citrus plots (R 2 c = 0.722), without problems of multicollinearity or stationarity [66] (Table 1).  We analyzed seven factors, whose descriptive statistics and spatial distributions are shown in Table A1 and Figure A1, respectively. The citrus plots had mean altitudes and slopes of 58.56 m and 1.50 • , respectively. The plots with the highest altitudes and slopes were located in the western and southern parts of the study area. Of the plots, 75.25% had a flat orientation, the average Dist-coast was 9.65 km (the plots furthest from the sea were 31 km away and the nearest were 0.5 km away), and the average Tair was 25.88 • C (with minimum and maximum temperatures of 24.88 and 26.627 • C, respectively). From the spectral indices, an average NDVI of 0.34 and an average NDWI of 0.07 were obtained (with maximum values up to 0.67 and 0.31, respectively).

LST Models
Based on the factors selected to explain the ICE, a multiple linear regression model was first fitted using an OLS procedure. Using the backward stepwise selection procedure, the model was found to have three significant factors: NDVI, NDWI, and the distance of the plots to the coast (Dist-coast). This specification explained practically three-quarters of the LST variability in irrigated citrus plots (R 2 c = 0.722), without problems of multicollinearity or stationarity [66] (Table 1).  Table 1 presents the results of the first GWR model, with improved adjustment relative to the OLS model (higher R 2 c and lower AICc). Although the use of this model diminished the spatial autocorrelation, it was not completely corrected and there was still a spatial pattern clustered in the residuals, indicating that the model was mis-specified.
Given that Dist-coast, one of the three significant variables, itself had a spatial pattern of clustering (Moran´s I index = 0.695, P-value = 0.000), the estimation was performed again with the two remaining factors, both by OLS (OLS-2) and by GWR (GWR-2), obtaining the results presented in Table 2.
The OLS-2 model had a worse fit (R 2 c = 0.544) and continued to show spatial autocorrelation, while the GWR-2 model improved the fit (R 2 c = 0.77, AICc = 490.11) and neutralized the spatial dependence effect (random pattern). The results show that the Dist-coast variable caused a clustered spatial distribution of the residuals from the local model. Therefore, it was eliminated from the set of seven initial factors considered. Finally, the models were estimated with the remaining six, obtaining the results shown in Table 3. In this model, both OLS and GWR estimation show three significant explanatory factors: the two already present in the previous specifications-NDVI and NDWI-plus the altitude above sea level of the plot (Altitude). While the problems of spatial autocorrelation persisted in the OLS estimation, the residuals had a random pattern in the GWR estimation. In turn, the GWR model showed the best fit (R 2 c = 0.78, AICc = 480.17) of all the models presented.
The analysis of variance (ANOVA) of the GWR model (Table 4) confirmed that the GWR specification showed the best fit, where there was a significant reduction in the magnitude of the residuals. In addition, lower RMSE values in the GWR model indicate smaller differences between predicted and observed values than in the OLS model. It should be noted that non-linear specifications were tested, and the interactions between factors were not significant. Therefore, the GWR model met all the criteria of econometric validation and explained 78% of the variance of the LST, based on three variables. Hence, a higher NDVI or NDWI or lower altitude reduces the surface temperature of the citrus plots. Based on the mean values of the model coefficients, each additional NDVI point reduced the LST by 2.15 • C, each NDWI point reduced the LST by 7.78 • C, and each 100 m of altitude above sea level increased the temperature of the plots by 1.50 • C.
The GWR model also makes it possible to evaluate the changes in local relationships between the LST and its explanatory factors at all points in the sample. By studying the local fit for each of the citrus plots of the study area, local R 2 values of the GWR-3 model were obtained (as shown in Figure 3a). Figure 3a also shows that the plots closer to the coast, especially to the northeast, explained better the variation in the modelled LST. In addition, the variation in the local regression coefficients of the GWR-3 model (Figure 3b-e) reflected the spatial heterogeneity of the adjustments. Meanwhile, the representation of the standard residuals of the GWR-3 model (Figure 3f) confirmed the random spatial pattern.

LST Predictions based on Changes in Land Use and Land Cover
The structural utility of the GWR-3 model has been proved. This model allows us to identify the factors that explain the ICE of the agricultural system. In addition, this model has a predictive utility that allows determination of the effects of land-use changes on target ecosystem services [65].
For example, a decrease in available irrigation flows in the study zone [48] would lead to the substitution of irrigation crops for other rain-fed crops. In this case, the values of the spectral indices (NDVI and NDWI) are modified, while the plot altitudes remain invariant.
The values of both indices for the most representative rain-fed crops in the area (almond, herbaceous crops, and grassland) were obtained by remote sensing for plots close to these crops on the reference date. From these, the LST was predicted for the 216 plots in the sample, assuming that plots were occupied by each of these crops. This gave the results shown in Table 5, where the average value predicted by the model for the current situation of irrigated citrus is also included. The difference in mean temperature was significantly higher for rain-fed crops, with respect to the reference of irrigated citrus. The average increase ranged from 1.20 • C to 1.22 • C, which indicates an increase in LST close to 4%.

Discussion
In this paper, the ICE of citrus cultivation was determined using multiple linear regression models. It was confirmed, using Moran´s I index, that the estimates made using OLS were not appropriate, given the spatial dependence of the observations. The use of the GWR technique, after the selection of appropriate factors, seems appropriate to eliminate the spatial dependence bias. The GWR gave a model with a better fit to the observations. Thereby, GWR models yielded significantly better predictions of the LST than OLS models, in agreement with the results of other studies [29,40,44,67,68] which used local models to explain the LST variable.
Thus, in the study of the ICE phenomenon, it has been observed that OLS models do not accurately identify the spatial variation in the LST in a heterogeneous agricultural environment, such as the one studied here. Compared to these conventional regressions, GWR models have the potential to allow for a better understanding of the ICE and the associated factors. Thereby, the spatial variation for each of the citrus plots considered (Figure 3a) could be visualized, with a better explanation of the reduction in the LST (higher local R 2 ) in plots close to the coast, and with a trend for the LST to decrease in a NE-SW direction. Likewise, the three factors found to influence the ICE differed in their behaviors among the plots. The intensity of temperature reduction was greater in the NE sector of this agricultural zone, when the NDVI was taken into account (Figure 3c), and in the central-eastern part considering the NDWI (Figure 3d). In terms of altitude (Figure 3e), a greater increase in the LST was observed in central areas. Thus, greater vegetation cover and greater water content of the vegetation favor the provision of climate regulation ES.
The inverse relationship of the LST and NDVI shown in this study is consistent with the results of numerous papers [69,70], confirming the role of vegetation in the reduction of surface temperatures through the transfer of latent heat from the surface to the atmosphere by evapotranspiration.
In addition, the LST values showed an inverse relationship with the NDWI [71], demonstrating the importance of the water held in crops in reducing the surface temperature.
Regarding the third factor-the altitude of the plot-the result was, in principle, counter-intuitive, where higher elevation of a plot above sea level was associated with an increase in the LST. Nevertheless, given the special characteristics of the region, with low levels of altitude (most of the plots coincide with the coastal plain) and the proximity to the coast, this effect was expected. Moreover, given the low average altitude in the area, there have been other studies [41] that confirmed a direct relationship between elevation and the LST in low-altitude areas.
The effects of changes in land use and land cover on the ICE were determined. Differences of up to 1.22 • C between irrigated citrus crops and dry land in the study area were found. These differences in the LST were due to the different thermal capacity, roughness, and surface albedo that each type of land use has [69]. These differences (between irrigated and dry land) were in line with differences found in other studies: Yang et al. [34] showed that irrigation cooled daytime LST by 1.15 K in China, and Lemus-Canovas et al. [72] showed that green areas in cities can decrease the LST by up to 2.5 • C, compared to urban areas. Furthermore, Liu et al. [73] showed that irrigated rice fields experienced a temperature reduction of 1.85 K, compared to rain-fed maize fields.
In order to evaluate the relative importance of the ICE, a comparative study of the ICE and the rest of the ES and EDS present in this agroecosystem was conducted. Using the results of a recent study [8] which identified significant ES (climate regulation, crop productivity, and biodiversity) and EDS (groundwater pollution and water supply) in this same agroecosystem and study area, the trade-offs were calculated. When land-use changes from dry land to irrigated citrus crops happen, the climate regulation ES (reducing the LST by 1.22 • C) is not the only one affected. Such a land-use change would lead to an increase in agricultural productivity, with an average value of 23,000 €/ha, and a reduction in biodiversity by 40% of bird species in this area, in spite of an increase in the diversity of soil organisms. On the other hand, EDS would also be produced, such as water consumption, with values around of 6,000 m 3 /ha and a groundwater pollution increase of 225 mg NO 3 /L.
In this way, it is possible to consider the effects that restrictions on the availability of water in the area may have on the LST, given the scarcity of this resource [74]. Irrigation water restrictions would transform the current citrus plots into dry land, as well as altering the cold effect phenomenon and the rest of ES and EDS provided by the agroecosystem.
Although this change in surface temperature may seem small (3.62% cooling), it has great importance in arid and semi-arid areas such as Campo de Cartagena, given its influence on agricultural yields. Higher temperatures shorten the crop cycle and, consequently, alter the local climatic suitability for specific crops [75]. Likewise, as is well known in agroecosystems, the provision of shade and reduction of temperatures allows for higher water infiltration and retention rates, regulation of the air quality, reduction of greenhouse gas emissions, and a greater variety of ecological niches [76]. Therefore, variation in the LST has an impact on the activity and diversity of soil organisms [77], the avifauna [78], and biodiversity in general [75]. In addition, several studies [10,79] have shown that the ICE of crops can generate a higher bioclimatic quality, comfort zones, and a reduction in the use of heating appliances in cities, thus saving energy and contributing to a decrease in greenhouse gases.
Thus, the decrease of LST by more than 1 • C is a significant reduction, within the ES framework [14,31,34,[71][72][73]80,81]. This is an important value to provide climate regulation ES in irrigated lands. Climate regulation is a crucial service in the set of ES provided by agroecosystems [17]; especially so under an uncertain future where climate change [10] makes the development of effective mitigation and adaptation strategies even more important. This is evident in the priorities of the Common Agricultural Policy of the European Union, which are intended to promote local measures to preserve the ES of agricultural ecosystems.

Conclusions
In this paper, we evaluated the provisioning of climate regulation ES by agroecosystemsspecifically, the ICE phenomenon-based on monitoring of the LST. To this end, a study considering the citrus crops in Campo de Cartagena, Southeastern Spain, was developed. Firstly, the existence of the ICE was confirmed by means of LST variations. Then the factors that determine LST were studied using the GWR, improving the results of OLS. The GWR showed the influence of the NDVI, the NDWI, and the altitude in the generation of this climatic phenomenon, while a reduction of up to 1.22 ºC in the LST with respect to dry land was quantified. This reduction represents a 3.62% decrease in LST, with respect to dry land.
The factors that influence the ICE have provided new knowledge to develop local policies to improve the provision of climate regulation ES, which are fundamental in mitigating the effects of climate change-one of the most important socioeconomic challenges facing society. A practical implication of the results is to provide a quantitative way to take into account the LST in the design of agro-environmental schemes, in the same way that biodiversity protection or CO 2 storage are evaluated. Furthermore, it allows us to take into account seasonal crops planning for climate regulation ES. This is made possible by analyzing the conditions that can influence the service flow provided. These implications provide agricultural and policy planners and managers with a new vision in the provision of climate regulation ES through the design of citrus plots with abundant vegetation and water content.
This study is not exempt of limitations. The LST varies over time and only a fixed thermal image was used. Previous studies [30,82] showed differences in LST between days and seasons; however, some factors (such as NDVI) may not have an effect on the LST during night-time or in the winter season. Therefore, the results of this spatial analysis should be complemented by further studies that analyze the LST multitemporal effect. Moreover, a spatial comparison between different agricultural areas under different climatic conditions would be recommended.
On the other hand, a recent study [8] on this same agroecosystem and study area evaluated the synergies/trade-offs between agroecosystem ES, indicating the significance of climate regulation ES compared to other ES. Thus, the quantification of the ICE is a necessary step in the full evaluation of the agroecosystem ES. Land-use change from dry land to citrus crops resulted in a reduction in temperature; however, at the same time, it modifies other related ES, the most representative synergies being increases in the crop productivity (by 18,852 €/ha/ • C) and water consumption (by 4918 m 3 /ha/ • C). Thus, socioeconomic benefits produced by the agricultural use of water are, therefore, broader than the traditional ones of food production or employment generation.
These results can be generalized to other semi-arid agricultural regions dedicated to irrigated citrus cultivation, as in the Spanish Levant or larger areas of the Mediterranean coast. This work establishes baseline values for the LST and its differences with other agricultural uses, in order to compare these results with other zones. Thus, this methodology constitutes an appropriate tool to complete the evaluation of agroecosystem ES, given that it allows the study and quantification of LST-one of the basic pillars (although less studied) of the climate regulation ES. Funding: This research was funded by AGRISERVI project, grant number AGL2015-64411-R financed by "Spanish Ministry of Economics and Competitiveness (MINECO)", and the 20912/PI/18 project, financed by "Fundación Séneca-Agencia de Ciencia y Tecnología de la Región de Murcia". The authors are also grateful for the financing of the first author´s pre-doctoral research by the "Spanish Ministry of Education, Culture and Sport (MECD) (FPU16/03562)".

Conflicts of Interest:
The authors declare no conflict of interest.