The Role of Irrigation Expansion on Historical Climate Change: Insights From CMIP6

To produce food for a growing world population, irrigated areas have increased from approximately 0.63 million km2 of land in 1900 to 3.1 million km2 of land in 2005. Despite this massive expansion, irrigation is still overlooked in most state‐of‐the‐art Earth system models (ESMs) involved in the Coupled Model Intercomparison Project phase 6 (CMIP6). To our knowledge, only three CMIP6 models represent irrigation activities: CESM2, GISS‐E2‐1‐G, and NorESM2‐LM. Here, we investigate the role of irrigation on historical climate at global and regional scales by analyzing trends of key surface climate variables in CMIP6 simulations during 1900–2014. The three models including irrigation show distinct behavior from the 15 models without irrigation over intensively irrigated areas (i.e., >50% of grid cell area is equipped by irrigation): both annual (months that correspond to monthly air temperature higher than 274 K) mean latent heat flux (LHF) and soil moisture increase over time, in contrast to the models without irrigation that show no trend or even a negative trend. The positive LHF trend over intensively irrigated areas in the irrigation‐on models is consistent with three satellite‐based LHF products. While annual (considering the warmest month in a year) warming trends are found in these regions for most of the no‐irrigation models, the increase in LHF induces a cooling trend for the models with irrigation, which was not confirmed by observational air temperature data sets. These findings, supported by satellite data, indicate the importance of improved representation of land management in the next generation of ESM.

, and have increased fivefold from approximatively 0.63 million km 2 of land in 1900 to 3.1 million km 2 of land in 2005 (Siebert et al., 2013). Irrigation induces biogeochemical and biogeophysical effects on climate (Mahmood et al., 2014;Sacks et al., 2009), and its expansion can alter precipitation patterns at the regional (DeAngelis et al., 2010) and global scales (Guimberteau et al., 2012;Puma & Cook, 2010). Irrigation amounts per unit irrigation area have been increasing over time just like irrigation extent has been Increasing (Siebert et al., 2015). Thiery et al. (2017) have shown that irrigation quantities are applied preferentially during hot (Mediterranean) or dry (South Asia) days.
The additional transpiration and evaporation associated with irrigation tend to cool the local climate . Moreover, irrigation perturbs the local surface energy and water budgets (Haddeland et al., 2006;Ozdogan et al., 2010;Pokhrel et al., 2011) and modulates local to regional atmospheric processes and landatmosphere interactions (Qian et al., 2020) via the impact of irrigation on the partitioning between sensible and latent heat fluxes (LHFs; Bonfils & Lobell, 2007;Kanamaru & Kanamitsu, 2008;Saeed et al., 2009;Sorooshian et al., 2014;Yang et al., 2020). Irrigation affects water resource availability over heavily populated irrigation hotspots (Vicente-Serrano et al., 2019). Despite the clear climate and water cycling impacts of historical and future irrigation (de Vrese & Stacke, 2020;Hurtt et al., 2020;Lawrence et al., 2016), it has generally not been considered in evaluations of drivers of observed climate change . The role of irrigation on land-atmospheric coupling and land surface model initialization has been demonstrated in several studies (H. Kim & Lakshmi, 2019;Lawston et al., 2015;Lu et al., 2015Lu et al., , 2017Pryor et al., 2016). Lu et al. (2017) found that irrigation led to a significant decline in coupling strength between the irrigated land and the atmosphere in the Midwestern United States. Lu et al. (2015) showed that the inclusion of realistic irrigation physics led to a better description of the land-atmospheric coupling. Lawston et al. (2015) demonstrated that the inclusion of irrigation in land surface models had the potential to improve forecasts, providing a better tool to adapt to increasing crop demands.
Earth observations and Earth system models (ESMs) have been widely used to study the past and potential future evolution of the Earth system processes that determine the exchanges between land surfaces and the atmosphere (Sellers et al., 1997). Incorporating irrigation among other land management practices (e.g., conservation agriculture, reservoir operation, and groundwater pumping) into ESMs has been highlighted as an important research direction (Hirsch et al., 2018;Mizuochi et al., 2021;Pokhrel et al., 2016). Despite intensive improvements in the representation of soil and vegetation processes in land models in ESMs (Cheruy et al., 2020), and although changes in the water cycle due to land and water management are now of global importance, most state-of-theart ESMs do not yet account for irrigation as they do for land cover change, aerosols, and other forcings (Cook et al., 2015).
Consequently, persistent biases in a number of key ESM variables (e.g., soil moisture and temperature) were observed in several regions (Al-Yaari et al., 2019;Cheruy et al., 2014). The existence of systematic hot and dry biases (precipitation deficit) in summer is, for instance, a known deficiency of many models participating in Phase 5 of the Coupled Model Intercomparison Project (CMIP5; Cheruy et al., 2014), particularly in regions of transition between humid and dry climates, where surface/atmosphere interactions contribute to the dryness of the future climate (Berg et al., 2014). This deficit could be due to the difficult characterization of the variables of the critical zone (Lahoz & de Lannoy, 2014), surface-atmosphere couplings (Cheruy et al., 2014), or, as mentioned above, the poor representation (or even an absence of representation) of irrigation activities in space and time in land modules within the ESMs (Qian et al., 2020).
Compared with that in CMIP5, the number of models participating in CMIP phase 6 (CMIP6) has increased greatly, and most models used an updated version (Cook, Mankin, et al., 2020;Danabasoglu et al., 2020;Hajima et al., 2020;Müller et al., 2018;Swart et al., 2019;Wu et al., 2019). However, as in CMIP5, most ESMs included in CMIP6 still neglect irrigation activities in climate projections. This negligence of irrigation activities in land models is due to a number of challenges ; for instance, the applied amount of irrigation water needs to be realistic, and there is uncertainty in the timing, techniques, and efficiency of irrigation practices. This, in turn, may lead to significant uncertainties in projections of future climate responses to irrigation.
However, progress has been increasingly made in this line of research to characterize changes in climate in response to irrigation. For instance, the Land Use Model Intercomparison Project (Boysen et al., 2020;Lawrence et al., 2019) was launched recently within the framework of CMIP6 to assess the biogeophysical effects of land use and land cover change (LULCC) in a multimodel framework. The majority of previous studies quantified 10.1029/2022EF002859 3 of 19 the impact of irrigation on weather, climate, and hydrology at local and regional scales based on a single model (Cook et al., 2015;Gormley-Gallagher et al., 2020;Thiery et al., 2017), focusing on temperature extremes (Hauser et al., 2019;Hirsch et al., 2017), or focused on incorporating irrigation into land-only hydrological models (Pokhrel et al., 2016).
To the best of our knowledge, no previous study has investigated the role of irrigation on historical climate using an ensemble of ESMs. Here, we take advantage of the newly available CMIP6 simulations from the most recent generation of climate models to re-examine the impact of irrigation activities at the global scale by investigating long-term trends of five essential climate variables (ECVs) that characterize the Earth's near-surface climate (evapotranspiration [ET], or the corresponding surface energy term, the LHF; precipitation; soil moisture; radiation components; and air temperature) over the 1900-2014 period. Within this study, we (a) compare models including and excluding irrigation (hereinafter referred to as CMIP6.irr and CMIP6.noirr, respectively) using CMIP6 historical (fully coupled) simulations; and (b) evaluate simulated trends against existing observational data sets.

Trend Analyses
For all CMIP6 variables except air temperature, we computed for each year the mean of the months with a monthly mean temperature higher than 274 K, assuming it approximately captures the growing season (Jiang et al., 2017;Prentice et al., 2011), when irrigation can be applied. One bias that the models could have is that irrigation quantities, even though realistic at the annual timescale , might be applied in different months/seasons than the ones where irrigation is actually happening (Jha et al., 2022). For the air temperature, the annual maximum of monthly air temperature (T max ) was considered instead, as previous studies have shown more significant irrigation impacts on hot extremes (Chen & Dirmeyer, 2020;Pitman et al., 2012;Thiery et al., 2020).
The impact of irrigation on the average of fluxes has already been studied for the models with irrigation versus their versions without irrigation (Cook et al., 2015;Thiery et al., 2017), so here we focus on the trends. To investigate the impact of irrigation on historical climate, the magnitude, or the rate at which trends are increasing/ decreasing is compared between CMIP6.irr and CMIP6.noirr within the 1900-2014 period. The widely used hydrometeorological time series (Tabari et al., 2011) Sen's slope nonparametric estimator (Sen, 1968) is used to estimate the magnitude of the trend (change per unit time) of the areas equipped for irrigation as well as the CMIP6 ECVs. The concept is to calculate the slope between all the pairs of ordered time points and then take the median of all these slopes as a unique estimator, representative of the entire time series. For a sample of N pairs: where x i and x j denote the corresponding data values at times i and j. The N values of S k are ranked from smallest to largest, and the median of Sen's slope is calculated as follows: A positive S med represents an increasing trend, and a negative sign represents a decreasing trend. Its value gives a quantitative estimate of the trend of the time series.
This test has low sensitivity to abrupt breaks due to inhomogeneous data series, does not require the assumption of normality, is robust against outliers, and can handle missing values (Chandler & Scott, 2011). The nonparametric Mann-Kendall test (Kendall, 1957) was applied to assess the strength and significance level of the trend. The trend is considered statistically significant (i.e., a significantly decreasing or increasing trend) if the p-value <0.05.

Irrigation Data Sets
To study historical irrigation expansion, we use the Historical Irrigation Dataset (HID) version 1.0 (Siebert et al., 2015). This data set provides gridded maps of the area equipped for irrigation expressed in hectares per grid cell based on national and subnational statistics at a 5 arc-minute resolution and a temporal resolution of 10 years from 1900 to 1980, increasing to 5 years from 1985 to 2005.
In addition, we use the latest version (5.0) of the FAO "Global Map of Irrigation Areas" (Siebert et al., 2013). This map (Figure 1a) quantifies the areas equipped for irrigation circa 2005 at a 5 arc-minute resolution (equal to 0.083°, i.e., ∼8 to 9 km approximately). It is a static map of irrigated areas, which are expressed as the fraction of each 5 arc-minute pixel area that is equipped for irrigation. This map was developed through a collaboration between the United Nation's FAO Land and Water Division and the Rheinische Friedrich-Wilhelms-Universität in Bonn, Germany. It relies on the data collected at national and sub-national level about irrigation water withdrawal per year (m 3 ), and is distributed via FAO's global water information system, AQUASTAT. The FAO "Global Map of Irrigation Areas" has been used as input to most continental to global scale models accounting for irrigation (e.g., de Rosnay et al. (2003), with an earlier version), and serves as a basis and/or benchmark for many recent studies aiming at proposing new irrigation maps (Meier et al., 2018).
It should be noted that this map is not exactly the same as the map in the historical irrigation data set for 2005 (HID). To compile the FAO map, the best and most precise estimate available for a year circa 2005 was used, whereby the reference year differs across countries. In the HID data set, in contrast, the authors tried to exactly match the situation in 2005 (see Figure S1 in Supporting Information S1 for the HID and FAO maps circa 2005).
The spatial values were classified into six levels at a 1° by 1° resolution, as illustrated in Figure 1a. The approximate fraction of the number of land pixels for each class (excluding Greenland and Antarctica) is illustrated in Figure 1c (right axis).
The evolution of the areas equipped for irrigation is illustrated in Figure 1b throughout the 1900-2005 period for each irrigation class (Figure 1a). Since 1950, there has been an abrupt increase in irrigation areas, mainly for grid-cell areas equipped for >5% irrigation. Figure 1c (left axis) shows the significant Sen's slope in the HID irrigation data sets: all irrigated areas expanded from 1900 to 2005, with the highest trends occurring in areas that were heavily irrigated. In irrigated grid cells where the irrigated grid cell fraction is greater than 50%, there is a rate of increase of 542 ha (spatially averaged) per decade from 1950 to 2005. In the subsequent analysis, the influence of this increasing trend on climate variables is presented and discussed. We interpolate the irrigation data sets onto 1° by 1° grid using bilinear interpolation (e.g., Stanfield et al., 2016).

Overview of CMIP6 Historical Simulations
CMIP6 is a primary source of data for climate change studies. In this study, historical simulations from 18 fully coupled climate models (listed in Table S1 in Supporting Information S1; constrained by data availability) were downloaded from the Earth System Grid Federation (ESGF) database (available at https://esgf-node. llnl.gov/search/cmip6/). These simulations were run under standardized protocols with uniform output and a consistent structure (Eyring et al., 2016). In this study, we use historical fully coupled (ocean-atmosphere-land) runs that are driven as much as possible by historical forcing data sets based on observations over 1850-2014 (Eyring et al., 2016). They include time-evolving data sets of atmospheric composition, solar forcing, and gridded LULCC, following the Land-Use Harmonization database (LUH2; Hurtt et al., 2020). However, not all ESMs represent the effects of LULCC with the same level of complexity, and most models are not able to make full use of the information provided in the LULCC data sets produced for CMIP6 (Hurtt et al., 2020;Lawrence et al., 2016).
The CMIP6 historical simulations start in 1850 until 2014, from which we examined the 1900-2014 period, which overlaps with the availability of historical irrigation data sets. These simulations with and without irrigation will be hereinafter referred to as CMIP6.irr and CMIP6.noirr, respectively. Historical simulations by CMIP6 models offer several members, and we focus on a single member for each model, the ones archived as "r1i1p1f1" within the ESGF databases, that is, with realization number 1, initialization method indicator 1, perturbed physics number 1, and forcing number 1. Other ensemble members were used for some models to check the internal variability among the different members.
To address the impact of irrigation on the simulated climate, we investigate five ECVs: soil moisture (hereinafter referred to as SM; in mm), annual maximum of monthly air temperature (T max ; in K), LHF (in W/m 2 ), precipitation (hereinafter referred to as Pr; in mm/day), and net radiation (hereinafter referred to as R net ; in W/m 2 ). Net radiation is computed as the balance of four radiation components (all in W/m 2 ), as follows: R net = SW net + LW net = (SW down − SW up ) + (LW down − LW up ). LW net is the net longwave surface radiation, and SW net is the net shortwave surface radiation, with down and up referring to surface downwelling and upwelling radiation, respectively. We focus on monthly means for SM, LHF, Pr, and R net and on the annual maximum of air temperature (T max ) and interpolate all data over the 1900-2014 period onto a common 1° by 1° grid using bilinear interpolation (e.g., Stanfield et al., 2016).

Main Features of the CMIP6.irr Models
Among the 18 selected climate models (Table S1 in Supporting Information S1), only three of them, namely CESM2, GISS-E2-1-G, and NorESM2-LM, include time-varying irrigation. These three CMIP6.irr models all respond to prescribed time-varying irrigated areas, which are prescribed as input data from the CMIP6 LULCC database (Hurtt et al., 2020). They were taken from the HYDE 3.2 database (Goldewijk et al., 2017), where irrigation is based on the HID data set (Siebert et al., 2015), shown in Figure S1 in Supporting Information S1. Therefore, the delineation of the areas where irrigation is occurring (current extent as well as historical expansion) is documented in the three CMIP6 models from the best available data sets. While CESM2 uses the Community Atmospheric Model version 6 (Danabasoglu et al., 2020), NorESM2 uses CAM6-Nor with some specific modifications made to the atmospheric component: an independent module for the life cycle of particulate aerosols, changes in the moist convection scheme and the local moist energy formulation, and a representation of aerosol-radiation-cloud interactions, as stated by Seland et al. (2020). Readers are directed to Seland et al. (2020) for more information about more differences between CESM2 and NorESM2-LM. GISS-E2.1-G is the coupled atmosphere-ocean-sea-ice-land NASA Goddard Institute for Space Studies (GISS) ESM (model E2.1), with G referring to the "GISS Ocean" model, which is an updated and improved version of GISS-E2-R used in CMIP5 Miller et al., 2021). GISS-E2.1-G describes three types of crops (Y. Kim et al., 2015), but it does not have dynamic vegetation, and all the parameters describing vegetation structure, including LAI, are prescribed at a monthly time step LAI (Cook, McDermid, et al., 2020).
As described by Cook, McDermid, et al. (2020) and Miller et al. (2021), irrigation-as the water required to ensure optimal crop growth-in GISS-E2-1-G is prescribed from an updated version of the time-varying irrigation water demand (IWD) data set of Wisser et al. (2010). These IWD data were generated by combining the best empirically constrained global state-of-the-art estimates of irrigation applications: a land-only terrestrial water balance model, the University of Frankfurt/FAO Global Maps of Irrigated Areas, and crop-specific calendars and water demand coefficients, which consider regional cropping practices. The IWD is calculated between 1900 and 2005 using the irrigation areal extent from Siebert et al. (2015) and is estimated between 1871 and 1900 by extrapolating backwards from the trend between 1900 and 1929. Irrigation is held constant after 2005 using the flux from 2004. Irrigation water to satisfy the irrigation demand is first sourced from lakes and rivers within the same irrigated grid cell. However, if the prescribed IWD cannot be fully satisfied by this flux, additional water is added from fossil groundwater, which represents an addition to the total mass of water within GISS-E2-1-G. Shukla et al. (2014) demonstrated that the simulated variability of summer monsoon circulation indices was improved by the inclusion of irrigation activities in GISS-E2-1-G. Nevertheless, irrigation in GISS-E2-1-G is applied on the entire vegetated grid cell fraction instead of only the crop fraction, causing all vegetation to respond to irrigation, which may lead GISS-E2-1-G to be too sensitive to irrigation, thus amplifying ET and the subsequent climate response.
Unlike GISS-E2-1-G, irrigation is internally calculated in CESM2 and NorESM2-LM. The irrigated crop water demand is triggered by a soil moisture deficit when there are growing crops (Lombardozzi et al., 2020;Sacks et al., 2009). The water for irrigation is initially sourced from river water storage (Danabasoglu et al., 2020); however, if this source is not fully sufficient for the demand, water is pulled diffusely from the oceans (Lawrence et al., 2018(Lawrence et al., , 2019. Similarly to GISS-E2-1-G, this methodology effectively assumes that irrigation is never water-limited, which would be the case in some situations in reality. Note also that the timing of irrigation is tied to the crop life cycle, which is prognostic and known to be biased in CLM5 in some regions, especially over India, where crop planting can occur too early in the year during the dry season (Lombardozzi et al., 2020). Furthermore, the irrigation parameterization does not distinguish between irrigation techniques (e.g., drip vs. flood), which may further bias the results. Thiery et al. (2017) evaluated the simulated irrigation amounts within CESM1.2.2. Locally, they may be overestimated or underestimated compared to observations and hence may overestimate or underestimate the climate response to irrigation. At the regional scale, the simulated irrigation amounts show a reasonable match with observed rates, with higher skills in South Asia. Most current irrigation implementations take a certain amount of water and apply it to the surface (either above or below the canopy). This water can either infiltrate, evaporate or run off. But in case of rice paddies, for instance, small basins are created where the water is actually stagnating. Current ESMs do not capture this effect at the moment, and therefore likely underestimate the effect of paddy irrigation on near-surface climate.
A common feature of the three climate models is that they do not consider groundwater as a possible water source for irrigation via groundwater pumping. In CESM2 and NorESM2-LM, irrigation is never limited by water availability, and if water needed for irrigation exceeds what is available in rivers, it is drawn from the oceans to ensure water conservation in the coupled climate system. Such unlimited irrigation is a way to circumvent the lack of groundwater pumping, since it helps achieving high irrigation rates even where groundwater is heavily used. Felfelani et al. (2021) recently developed a prognostic groundwater module for CLM5, accounting for aquifer pumping, lateral groundwater flow, and conjunctive use of groundwater for irrigation, and it will probably be incorporated in a future CESM version. Moreover, while dynamic lake area and dam management have recently been implemented in CLM5, these parameterizations currently do not influence irrigation water availability in the model (Vanderkelen et al., 2021).

Observation-Based ECVs
The observation-based products listed below serve as a benchmark between the two different CMIP6 categories: CMIP6.irr and CMIP6.noirr. These products were interpolated to the 1° by 1° grid used for the CMIP6 models.

GLASS
The Global Land Surface Satellite (GLASS) ET product is estimated using the Bayesian model averaging method (Yao et al., 2014). This approach merges five process-based ET algorithms, including the Moderate Resolution Imaging Spectroradiometer (MODIS) ET product algorithm, the remote-sensing-based Penman-Monteith ET algorithm, the Priestley-Taylor-based ET algorithm, the modified satellite-based Priestley-Taylor ET algorithm, and the semiempirical Penman ET algorithm of the University of Maryland. The GLASS ET products cover the 1982-2015 period with a spatial resolution of 0.05° by 0.05° and an 8-day period. GLASS ET is provided in mm/ month, which is almost equivalent to W/m 2 (as one needs ∼30 W/m 2 to evaporate 1 kg/m 2 of water).

GLEAM
The Global Land Evaporation Amsterdam Model (GLEAM) is one of the remote sensing-based ET and SM products (Martens et al., 2017) comprising a set of algorithms with different inputs of remote sensing observations and meteorological reanalyses. In this study, the monthly GLEAM v3.0 ET product with a spatial resolution of 25 km between 1982 and 2014 was considered. GLEAM ET is provided in mm/month, which is almost equivalent to W/m 2 . In GLEAM, potential ET is computed using the Priestley-Taylor formula based on near-surface temperature and net surface radiation measurements. Second, potential ET is converted to actual ET using a multiple stress factor as a function of remotely sensed soil moisture and vegetation optical depth (Martens et al., 2017).

Fluxnet-MTE
This global LHF product has been produced by upscaling the current Fluxnet-MTE (multitree ensemble) observations of LHF using machine learning techniques and model tree ensembles (Jung et al., 2011). The gridded Fluxnet-MTE data set is derived from continuous in situ measurements of Fluxnet-MTE, remote sensing, and meteorological observations based on model tree ensembles. These data have been extensively validated using observations from eddy covariance systems, and the product has been widely used (Thiery et al., 2015;Zeng & Zhang, 2020). These data sets are provided in a monthly time step in MJ/m 2 /d (0.0864 MJ/m 2 /d = 1 W/m 2 ) at a spatial resolution of 0.5° during the 1982-2011 period. It should be noted that a new product was recently released, namely, FLUXCOM (Jung et al., 2019). However, caution was warranted by Jung et al. (2019) when studying trends due to issues with sensor age-based drift in MODIS reflectance. Therefore, this product was not considered in this study.

Climate Research Unit Gridded Time Series (CRU)
We use the updated version of the CRU TS v4 monthly near-surface air temperature data over the period 1901-2014. The spatial resolution of the CRU data sets is 0.5° by 0.5° over the global land surface except for Antarctica. It is derived by the interpolation of monthly climate anomalies using angular-distance weighting from extensive networks of weather station observations (Harris et al., 2020). These data were used to screen out GISS, GLEAM, and Fluxnet-MTE monthly values that correspond to CRU monthly air temperatures lower than 274 K.

Willmott Gridded Time Series (Willmott)
We use the University of Delaware monthly mean surface air temperature data (Willmott et al., 2001) over the period 1900-2014. The spatial resolution of the Willmott data sets is 0.5° by 0.5° over the global land surface. It is derived extensively from the archive of Legates & Willmott (1990) and from the Global Historical Climate Network (GHCN2).

CMIP6 Simulations
Historical anomaly time series (relative to a 1900-1920 baseline period) of spatially averaged values of all pixels within the 50%-100% irrigation class for the variables LHF, SM, Pr, T max , R net , LW net , and SW net are shown in Figure 2. Most of the pixels of this irrigation class are situated into South Asia (see Figure 1a), and we can see that the abrupt change (e.g., LHF and SM) starts from approximately 1960, when a massive irrigation expansion took place as part of the Green Revolution in the 1970s in India (Mishra et al., 2018). Irrigation expansion was enabled by the important rise in the number of groundwater wells (Ambika & Mishra, 2020). In the CMIP6.irr models, the increase in LHF (with a rate of 3.4 W/m 2 /decade) and SM (with a rate of 0.72 mm/decade) resulted in a cooling of the surface (by a rate of −0.16 K/decade) during the 1960-2014 period. This is in line with a recent observational study (Ambika & Mishra, 2020), which demonstrated that irrigation modulates SM and enhances LHF, causing noteworthy cooling in India. This evaporative cooling mechanism was described by Boucher et al. (2004) along with increased absorption of solar radiation. Moreover, Roy et al. (2007) showed a 3-4°C daytime cooling in irrigated areas between the pre-Green Revolution (unirrigated) and post-Green Revolution (irrigated) periods. There is a slight increase in Pr (with a rate of 0.06 mm/day/decade) simulated by CMIP6.irr compared to CMIP6. noirr. The mostly irrigated class mostly corresponds to India, which is a recognized hotspot of land-atmosphere coupling (Koster et al., 2004). The response of rainfall to land conditions is therefore strongly modulated by atmospheric boundary layer and convective parameterizations (Cheruy et al., 2013), but the CMIP6 experimental design does not readily permit assessment of these feedbacks.
Based on SW net , there was a solar dimming trend for both irrigated and non-irrigated CMIP6 model types. Cook et al. (2015) related this dimming to irrigation-induced increases in cloud cover. Moreover, the dimming trend over India was recently confirmed using the Global Energy Balance Archive data (GEBA), which is based on ground-based stations measuring energy fluxes (Moseid et al., 2020). The latter study validated NorESM2-LM, CanESM5, MIROC6, CESM2, and CNRM-ESM2-1 SW net simulations against the GEBA observations and found similar dimming trends among all models and observations. They found a dimming trend, which is also consistent with several other observational studies (  This dimming is mainly attributed to an increase in aerosol optical depth (Soni et al., 2016). However, the SW net trends in CMIP6.irr and CMIP6.noirr are indistinguishable; therefore, there does not appear to be a link between dimming and irrigation. In contrast, there is an increase in the LW net simulated by CMIP6.irr (with a rate of 2.9 W/m 2 /decade), which is reflected in R net anomalies. The variability (shaded area in Figure 2, which represents ± one standard deviation) of CMIP6.irr is larger than that of CMIP6.noirr. Figure 2 is reproduced in Figure S2 in Supporting Information S1 with showing the different values of the three models of CMIP.irr. This variability comes from GISS.E2.1.G as the other two models, which share the same land model and irrigation scheme, tend to agree on each other. Moreover, over this irrigation class, the mean of all variables during 2005-2014 for the three models in CMIP6.irr and the ensemble CMIP6.noirr is shown in Figure S3 in Supporting Informa tion S1. The mean values of LHF and SM of the three models (CMIP6.irr) are higher than the ensemble CMIP6.noirr. For T max , GISS.E2.1.G and NorESM2-LM are smaller than CMIP6.noirr. While NorESM2-LM has the lowest value for Lwnet, CESM2 has the lowest values for R net and Swnet. This variability among the CMIP6. irr models may be explained by the different parametrizations of the irrigation implementation.
Time series of spatially averaged values of all pixels within the other irrigation classes (i.e., 0%-50%) for LHF, SM, Pr, T max , R net , LW net , and SW net anomalies (i.e., subtracted from 1900 to 1920) of CMIP6.irr on and CMIP6. noirr during the period 1900 to 2014 are depicted in Figures S4-S8 in Supporting Information S1. As an example, the behavior of CMIP6.irr and CMIP6.noirr LHF over all irrigation classes is illustrated in Figure 3. Over regions where irrigated area is essentially zero (i.e., 0%-0.1%), the temporal dynamics of time series anomalies are very similar in both CMIP6.irr and CMIP6.noirr, with an increase since circa 1960. Over this class, the trends of all studied variables are consistent between CMIP6.irr and CMIP6.noirr ( Figure S4 in Supporting Information S1). Figure 3 shows that LHF starts to diverge between CMIP6.irr and CMIP6.noirr in the least irrigated class, with slightly lower values of LHF in CMIP6.irr than in CMIP6.noirr. The difference in LHF between the two ensembles becomes clear in the third irrigation class (i.e., 5%-10%), and this difference between CMIP6.irr and CMIP6. noirr LHF increases with the fraction of equipped areas for irrigation, which is confirmed for most variables by Figures S4-S8 in Supporting Information S1. This implies that the differences shown in Figure 2 between CMIP6.irr and CMIP6.noirr over the intensively irrigated regions can be attributed to the irrigation processes. It should be noted that the monthly mean of another variable named "daily maximum temperature" (T daymax ) showed similar anomalies ( Figure S9 in Supporting Information S1) as T max . In this figure, however, CMIP6.irr includes only one model, GISS.E2.1.G, as the other two models do not provide this variable. Given the similarity of the results using T daymax or T max , and as most of CMIP6 models do not provide "T daymax ," T max will be presented in the subsequent analyses. Figure S10 in Supporting Information S1 reproduces Figure 3 showing the different values of the three models in CMIP6.irr. The three models tend to agree over irrigation classes up to 50% but GISS. E2.1.G, again, shows higher values particularly over the recent years over the irrigation class 50%-100%.
To better understand the behavior of the CMIP6 models individually and since the substantial change of the anomalies shown in Figure 2 started in 1960, Sen's slope was quantified (per pixel) over the 1960-2014 period for the LHF, SM, Pr, T max , R net , LW net , and SW net variables simulated by all CMIP6 models used in this study. The resulting Sen's slope values were spatially averaged based on the FAO irrigation map classes (see Figure 1a). Over this period, while the values of LHF and SM have increased in CMIP6.irr models (i.e., CESM2, GISS-E2-1-G, and NorESM2-LM), they have decreased or have no trend for all other CMIP6.noirr models (as illustrated in Figure 4). This clearly suggests that irrigation has a nonnegligible effect on the water and surface energy balance over intensively irrigated areas. The largest increases are observed over irrigated grid cells where the fraction irrigated is greater than 50%. Over this irrigation class, T max decreased for CMIP6.irr (a slight cooling trend), whereas it increased for the majority of the other CMIP6.noirr models. This divergence between CMIP6.irr and CMIP6.noirr over these regions (mostly grid cells located in South Asia) can most likely be attributed to irrigation, as it is one of the major differences between these two categories of simulations.
For the historical evolution of R net , the effect of irrigation is only detectable over highly irrigated regions. CMIP6. irr simulations begin to diverge from CMIP6.noirr above the irrigation fraction class 20%-50%. There are similar trends of both the CMIP6.irr and CMIP6.noirr models for SW net . This indicates that irrigation has a much larger impact on surface turbulent flux partitioning (latent heat) than on surface shortwave radiation (i.e., clouds), which is in line with previous studies (Cook et al., 2015). No clear influence of irrigation on precipitation was detected. Previous studies (Alter et al., 2015;Huber et al., 2014;Qian et al., 2020) have pointed out that irrigation increases precipitation over or downwind of the irrigated regions, but we do not evaluate precipitation trends associated with increasing irrigation in this study, partly because we anticipate that it will be difficult to identify a clear signal within the context of generally more noisy precipitation trends. Although the representation of irrigation is identical in NorESM2-LM and CESM2, the simulated trends such as the T max trends are not. This suggests that the impact of irrigation on the simulated climate can differ between models despite identical irrigation implementation. These differences may originate from substantial differences in the model setup, structure and atmospheric parametrizations (de Vrese & Hagemann, 2018). Although CESM2, GISS-E2-1-G, and NorESM2-LM implement an irrigation scheme, GISS.E2.1.G shows the strongest increase in LHF, SM, and LW net over all irrigation classes, followed by NorESM2-LM and CESM2. This is not surprising owing to the different quantities of water applied for irrigation, as there are no common data sets of irrigation forcing, defined standards, or implementation protocols for irrigation in the three climate models (see a full description of the irrigations schemes and implementations in the models in Section 2.3.2).
Potentially, the observed trends in CMIP6.irr may be smaller than the internal variability of CMIP6.noirr trends in CMIP6.irr may be masked in CMIP6.noirr by internal variability within the individual ensemble members. Under this hypothesis, the positive trend of LHF, for instance, observed over the 1960-2014 period in models with irrigation, may just be one particular realization ("r1i1p1f1") among all the possible ones. Luckily, several CMIP6 models provide multiple-member ensembles. To examine how the modeled trends can differ due to simulated internal climate variability, we assessed the ensemble members of one CMIP6.irr model (CESM2) and one CMIP6.noirr model (IPSL-CM6A-LR). In addition, LUMIP includes specific runs (Histo-noLu) to evaluate individual land-use impacts, including a simulation without irrigation changes (i.e., irrigated area and all other LULCCs held constant at 1850 levels). This run does not isolate irrigation, although at the high-irrigated area level, it is probably mostly an irrigation signal. This enables direct comparison of the same model with irrigation on and off.
Trends over the 1960-2014 period from 11 members of the CESM2 fully coupled historical ensemble are depicted in Figure 5. All coupled historical ensemble members agree well for the LHF, T max , and SM and disagree with the Histo-noLu simulation. The behavior of CESM2 Histo-noLu is similar to CMIP6.noirr (Figure 4). Similar to CESM2, trends over the 1960-2014 period from nine IPSL-CM6A-LR coupled historical ensemble members are depicted in Figure S11 in Supporting Information S1. Based on the results of IPSL-CM6A-LR and CESM2, trends are consistent across the ensembles for most fields, indicating that the forced trend due to climate change or climate change plus irrigation in the case of CESM2 dominates over the internal variability.

Observational Validation
Our results show that the differences in LHF and T max between the CMIP6.irr and CMIP6.noirr simulations can be highly linked to irrigation influences. Therefore, to validate whether they make the CMIP6.irr climate models more realistic, we compared the simulated time series of LHF and T max to those of observation-based data sets.
We first focused on heavily irrigated regions (i.e., 50%-100% class in Figure 1a) and used three observation-based estimates (Fluxnet-MTE, GLASS, and GLEAM) for LHF, with a common availability period over 1983-2014. The validation period for T max is longer, and we consider the 1960-2014 period (used in Figure 4) against CRU and Willmott. For both variables, simulated and observed time series are offset to share the same mean over the beginning of the analyzed period for easier trend comparison. Figure 6a shows the resulting anomalies (relative to 1983-1989) along with those of the averaged CMIP6.irr and CMIP6.noirr models. While the three observation-based data and the CMIP6.irr models showed an increase in LHF, the CMIP6.noirr models showed no increase or decrease but remained steady during the 1983-2014 period. CMIP6.irr (3.62 W/m 2 ) is closer to the three observations (Fluxnet = 2.38, GLASS = 4.15, and GLEAM = 3.02 W/m 2 ) than both CMIP6.all (0.88 W/m 2 ) and CMIP6.norir (0.33 W/m 2 ) in terms of magnitude. The better match between CMIP6.irr and observations is in line with a recent observational study that showed increasing trends of LHF over western India (Pan et al., 2020). The robust positive LHF trend in intensely irrigated regions obtained across a suite of independent observational products adds confidence to the finding that accounting for irrigation increases the ESM realism in terms of near-surface climate representation in intensely irrigated regions. It should be noted that none of the LHF observations used here to evaluate the CMIP6 models account explicitly for irrigation but that in all cases, the effects of irrigation on LHF are indirectly considered via the observations used to retrieve LHF (e.g., remotely sensed soil moisture and LAI). Furthermore, looking at the time series anomaly of the ensemble mean of all CMIP6 models with and without irrigation (CMIP6.all in Figure 6a), the separation of the CMIP6 models according to irrigation gave better agreement with the observations than taking the ensemble mean of all CMIP6 models. This indicates that, in contrast to previous studies (Flato et al., 2013;Loew et al., 2016), the ensemble mean, over intensively irrigated areas, is not always a better  Figure 1a). The shaded area represents ± one standard deviation of each product.
choice. Therefore, caution should be exercised when considering the ensemble mean of all different CMIP6 models. Figure 6b shows the annual maximum of monthly air temperature time series anomalies (relative to 1940-1959) for CMIP6.irr and CMIP6.noirr along with CRU and Willmott over heavily irrigated regions (i.e., 50%-100% class in Figure 1a). Both Willmott and CRU observations show almost the same variability with relatively small differences. The CMIP6.irr models correctly capture the weak decrease in air temperature that started around 1960 (tend to agree slightly better with Willmott). The entire period is subject to global warming, which becomes clearly discernable in both kinds of simulations after 1980 in the non-irrigated areas ( Figure S4 in Supporting Information S1), and in the other irrigation's classes for the CMIP6.noirr simulations only ( Figures S5-S8 in Supporting Information S1). The CMIP6.irr simulations, in contrast, show more complex evolutions, with a decrease of T max between ca. 1950 and 1970, probably caused by irrigation intensification, while the T max increase afterward suggests that the cooling from irrigation is offset by global warming.
However, despite the variability in the observed T max , we can estimate that the latter starts to diverge from the CMIP6.irr simulations around 1980, when the observed T max is rather stable but over the last decade, showing a strong warming trend. The CMIP6.irr models simulate a lower T max than the observations, with a decreasing trend, which is not realistic. The CMIP6.irr simulations are not realistic either, with an excessive T max over the last two decades, resulting from a too early warming trend. In highly irrigated areas where the cooling effect of irrigation conflicts anthropogenic warming, our results suggest that the observed T max could be better captured by climate models with irrigation, but with reduces rates. This analysis is consistent with the fact that GISS-E2-1-G overestimates irrigation demand, and all three models assume unlimited water supply. But given that the simulated LHF is well captured by the CMIP6.irr models in the highly irrigated areas, the mismatch of simulated and observed T max trends may also be due to biases in the three climate models. It could be also that the irrigated models can be wrong in their magnitude and may overestimate irrigation and the resulting ET increase, thus cool the surface layers more than a with a more realistic irrigation. To expand this analysis, Figure 7 shows the 1983-2014 trends for both LHF and T max over all irrigation classes. Except for the least irrigated class, divergence is detectable starting with irrigation classes 5%-10%: the simulated LHF trends agree better with observational data in CMIP6.irr than in CMIP6.noirr. They remain underestimated in CMIP6.irr, however, for the most highly irrigated class (50%-100%), in which the mean trend of CMIP6.irr simulations (2.57 W/m 2 /decade) very satisfactorily matches the observed ones (GLASS with 2.63, GLEAM with 2.54, and Fluxnet-MTE with 2.24 W/m 2 /decade).
As already mentioned, the analysis of the trends of annual maximum of monthly temperature (T max ) over 1983-2014 is complicated by anthropogenic warming. In agreement with global land estimates (Dunn et al., 2020;, it is approximately +0.2 K/decade according to the observation data (CRU and Willmott), based on the least irrigated class (<5%) where the simulated temperature trends are hardly sensitive to irrigation (larger trends for CMIP6.irr, probably owing to large model dispersion, as shown in Figure 4 for T max ). This posi tive trend in observed T max significantly decreases with irrigation intensity, to be less than 0.05 K/decade in the most intensively irrigated class. This response of the warming trend to irrigation intensity is hardly captured by climate models that ignore irrigation (CMIP6.noirr), especially given the large model spread. In the models implementing irrigation (CMIP6.irr), the warming trend is reduced with extensive irrigation, but was excessive compared to observations, that may lead to an unrealistic cooling trend during 1983-2014 in the mostly highly irrigated class.
The reduced warming trends found in CMIP6.irr and observations over highly irrigated areas are consistent with the alleviation of hot extremes by expanding irrigation already found by Thiery et al. (2020) by combining observed temperature and irrigation data with historical simulations with a previous version of CESM2. A likely explanation is that irrigation increases LHF, which alleviates extreme temperatures, a relationship that has been observed in India (Bonfils & Lobell, 2007;Chou et al., 2018) and China Zhu et al., 2012), although this cooling may potentially not be solely attributed to irrigation alone (van Oldenborgh et al., 2018). The excessive responses of LHF and T max to irrigation in the CMIP6.irr models may be related to imperfect implementations of the irrigation practices. Additional detection and attribution research are required to analyze how much of the extra cooling in CMIP6.irr is due to irrigation. A proper detection and attribution analysis would require running dedicated factorial experiments with and without transient irrigation extent over the historical period with multiple earth system models." Such type of analysis is beyond the scope of our study. Moreover, the gridded temperature products may underestimate the irrigation effect, as stations are typically not located in irrigated croplands and thus may miss the local cooling influence of irrigation.

Conclusions
Changes in land use and land cover can alter the biophysical properties of the land surface, including albedo, LAI, soil wetness, and roughness, which in turn influence the exchange of water and energy between the land surface and atmosphere. As already shown by Cook et al. (2015), the expansion of irrigated areas influences the long-term global terrestrial water balance via these biogeophysical effects. The present paper generalizes this seminal study by analyzing historical long-term simulations (1900-2014) of LHF, SM, Pr, T max , R net , SWnet, and LWnet using fully coupled CMIP6 climate models. Three models among the CMIP6 models account for irrigation activities (i.e., CMIP6.irr), while the remainder omit irrigation (CMIP6.noirr). The striking feature is the clear and visual distinction between the CMIP6.noirr and CMIP6.irr models in heavily irrigated regions. Based on time series anomalies and Sen's trend results, a distinct behavior of the three models was observed, with annual LHF and SM generally increasing over heavily irrigated regions. The trend in the corresponding annual maximum of monthly temperature over the irrigated areas showed a cooling trend. We also find that neglecting irrigation reduces the accuracy of LHF trend estimates over heavily irrigated regions, as the behavior of models with irrigation was consistent and confirmed by observational LHF products. In the intensively irrigated areas, the CMIP6.irr models capture the decrease in air temperature that started in 1960, observed by observational data sets. However, the cooling shown by CMIP6.irr after 2010 was not confirmed by the observational air temperature products. Furthermore, all ensemble members of the two CMIP6 models showed similar trends for all the variables, ruling out the possibility that the difference between CMIP6.irr and CMIP6.noirr is due to natural variability. Based on the recent CMIP6 experiments (with variables considered in this study) and their comparison to observational products, our results highlight the importance of irrigation as a historical climate forcing over intensively irrigated areas. None of the models considered in this study, however, currently differentiates between irrigation techniques, that is, they all apply one single irrigation technique across the globe (typically adding a calculated amount of extra water to the soil surface). This is an evident shortcoming of these models, which exhibit significant biases relative to observed/reported estimates of irrigation amounts. These models also overlook the reliance of irrigation on groundwater supply in many areas (Siebert et al., 2010), and the potential limitation of irrigation in cases of surface or groundwater shortage. Further improvements to the irrigation schemes and underlying land models are thus needed in the future generation of climate models, to better account for past and future evolutions of water resources, near-surface climate, and also food production over irrigated areas, expected to further expand within the next decades because of climatic and demographic pressures (Puy et al., 2020). Such an effort will need to face the challenge of representing realistic irrigation practices (e.g., timing, source of water, amount of irrigation forcing applied, and method of application) into climate models and comparing these parameterizations across modeling groups.
Large increases, compared to present, in global irrigated agriculture that could double by 2050 (Puy et al., 2020), and are projected in 20100 by two future scenarios (SSP3-7.0 and SSP5-8.5), with a range of values between 3.4 and 4.1 million km 2 (Hurtt et al., 2020). Irrigation areas drive irrigation volumes (Puy et al., 2021), but this relationship will not probably remain the same, as water can become limiting, either because irrigated areas have expanded too much or because climate change decreases the available water resources. Future research could consider CMIP6 projections to assess the impact of future changes in irrigation and its impact on climate variables. Once changes in irrigation are adopted by the majority of models as a standard historical anthropogenic climate forcing, its importance for historical changes in climate can be directly compared to other anthropogenic forcings (e.g., greenhouse gases and aerosols) and other LULCC anthropogenic activities (Sterling et al., 2013) that influence ET and therefore the near-surface climate (e.g., fertilizers, urbanization, etc.). Finally, an important issue is to properly quantify the relative weight of these factors of change in the past to try to project them into the future in a detection and attribution framework.