Impacts of land cover changes and global warming on climate in Colombia during ENSO events

Colombia is highly vulnerable to climate change which may be intensified due to the climatic effects of regional deforestation. Here, we quantify the impact of historical (1900–2011) land cover changes (LCC) and of global warming during ENSO events (CC) on precipitation, temperature and surface energy balance components by running the Weather Research and Forecasting model WRF v3.9 at 10 km resolution. We find that historical anthropogenic CC causes a mean temperature increase of 0.77 ± 0.02 °C in Colombia, which is more pronounced in high altitudes. Precipitation is enhanced by 0.98 ± 0.30 mm/day (+ 9%), particularly over forested areas and reduced at the Pacific coast. LCC imply a reduction of precipitation particularly above the Andes (− 0.48 ± 0.10 mm/day) and Caribbean Coast (− 0.67 ± 0.12 mm/day), where LCC effects dampen CC effects by 24% and 72%, respectively. La Niña tends to intensify LCC and CC effects in the Andes but dampens them at the Coast, roughly by a factor of two compared to El Niño impacts in both regions. At the subregional level, LCC and CC can have impacts of similar magnitude on precipitation highlighting the need to precisely account for both drivers in hydroclimatic assessments. Contrary to almost all observations and similar simulations with climate models, WRF simulates a cooling bias after historical deforestation in Colombia, even with alternative WRF land surface models. We identify two main sources of biases in the default WRF parametrization to explain this inaccuracy: (1) surface shortwave radiation reflected after deforestation is overestimated; (2) associated evapotranspiration loss is underestimated. Improved model representation and validation of tropical vegetation properties are necessary to provide more robust and confident projections.


Introduction
Since preindustrial times, mean surface air temperature has risen globally by 1 °C due to anthropogenic emissions of greenhouse gases from fossil fuel burning and land use and land cover changes (LCC, IPCC 2018). Especially the tropics are recently experiencing accelerating deforestation rates (Hansen et al. 2013, Song et al. 2018. Deforestation influences climate by reducing carbon stocks and releasing carbon into the atmosphere (Pan et al. 2011), leading to global warming (biogeochemical feedback to deforestation) and by alterations of regional energy, momentum and water fluxes (biophysical feedbacks), which can outbalance biogeochemical effects and lead regionally to net cooling or additional warming (Claussen et al. 2001). Jia et al. (2019) suggest that on a global scale, biogeochemical effects dominate the interplay of both feedbacks. However, the net impacts of LCC highly vary in space and time and may have contrasting consequences depending on latitude (Davin and de Noblet-Ducoudré 2010). Modelling and observational studies suggest that boreal regions experience a cooling of surface air temperature, while lower latitudes show warming in response to local deforestation (Alkama and Cescatti 2016;Findell et al. 2017). Colombia, a megadiversity hotspot and one of the top 12 deforesting countries (FAO 2015), is particularly affected by such LCC impacts, which in tropical countries are expected to amplify the warming from CC. Every year, roughly 180,000 ha (IDEAM 2017(IDEAM , 2018, 2019) of 1 3 forests in Colombia are lost due to agricultural expansion, cattle ranching, urbanization and timber extraction. At the same time, the intensity of extreme weather events associated with the El Niño-Southern Oscillation (ENSO), e.g. droughts or floods, is increasing (Ávila et al. 2019;Hoyos et al. 2013). Thus, the question arises of how such climate events and LCC are interacting.
Understanding deforestation and global warming impacts on climatic variables in Colombia is crucial because this country concentrates very high levels of ecosystem services, terrain complexity, species richness and water provision. Modelling this region is challenging because global climate models cannot resolve small-scale terrain complexity inherent to the presence of the three Andean cordilleras (Posada-Marín et al. 2018;Espinoza et al. 2020). Regional Climate Models (RCMs) can represent local, small-scale phenomena more accurately than global General Circulation Models (GCMs) or Earth System Models (ESM). One of the most widely used RCMs is the Weather Research and Forecasting model (WRF) which has proven its capability to dynamically downscale global reanalysis data from coarse to high resolution in good agreement with observations (e.g. Posada-Marín et al. 2018;Caldwell et al. 2009;Lee and Berbery 2012). In complex terrain domains, WRF is considered to produce climatic features even better than the original reanalysis data sets due to more accurate representation of orography (e. g. Heikkilä et al. 2010;Sun et al. 2016;Soares et al. 2012;Gao et al. 2012;Jiménez-Esteve et al. 2018).
To date, no impact study for Colombia exists that evaluates the impacts of LCC on climate at high spatial resolution. Thus, we present here the first high-resolution, dynamical downscaling study investigating LCC and CC impacts in Colombia. We perform sensitivity simulations with WRF v3.9, forcing the model with historical land cover and climate input and compare the results against present-day conditions. Our aims are to answer the following questions: (1) How do historical LCC and CC influence surface air temperature and precipitation in different regions of Colombia, during ENSO events (2009 (2) What is the dominant forcing for regional climate changes in Colombia? (3) How are land cover changes and global warming impacts modulated by ENSO? (4) How is the surface energy balance affected by different scenarios of land cover change in Colombia?

Study site
Colombia is located in the northwest of South America and contains five major biogeographic regions with contrasting biophysical and land cover properties: Colombian Amazon, Caribbean, Pacific Coast, Orinoco plains (Llanos orientales) and Andes. The three branches (cordilleras) of the Andes, which cross the country in the center and west, confers Colombia a unique and complex orography with altitudes ranging between 5800 m and sea level, resulting in a wide variety of climatic and vegetation zones. Localscale phenomena such as inter-Andean valley circulations but also large-scale oscillations in the Atlantic and Pacific, interactions with the Amazon and Orinoco basin influence Colombia's climatic variability (Poveda et al. 2011;Espinoza et al. 2020). On the inter-annual scale however, the hydrological variability is dominated by ENSO, with El Niño causing less evapotranspiration followed by less recycled precipitation, consequently fewer river discharges, and less cloud cover which facilitates more radiation and eventually higher temperatures (Poveda et al. 2011). La Niña, the opposing ENSO-phase, has a reversed effect. ENSO impacts the western Andean branch first while the eastern branch is affected with a little time lag and less strong (Poveda et al. 2011). Generally speaking, Colombia experiences a drier and prolonged dry season during El Niño, and a wetter and prolonged wet season during La Niña.

Model settings
In order to simulate climate response to LCC and CC in Colombia, we use WRF version 3.9.1.1 (Skamarock et al. 2008). This non-hydrostatic, fully compressible and terrain-following sigma coordinate model (Powers et al. 2017) has been proven to perform well in highresolution regional modelling studies and is one of the world's most widely used numerical weather prediction models. The simulations are executed by the Advanced Research WRF (ARW) dynamical core of wrf.exe. It contains basic dynamical equations for advection, Coriolis, pressure gradient terms, buoyancy and diffusion. Options on microphysics, radiation schemes and planetary boundary layers are set in namelist.input. Lookup tables for soil characteristics and land cover and vegetation properties set important physical features like minimum and maximum albedo per land use category, surface roughness, rooting depth, which the model accesses during the model run.
The simulations are carried out on two regional domains encompassing a large external domain (D01) covering large parts of Central America, Northern South America, Amazon basin, the Caribbean Sea and Eastern Pacific Ocean, along with an internal (nested) domain covering continental Colombian territory (D02, see Fig. 1), excluding the Colombian islands of San Andres, Providencia and Santa Catalina off the coast of Nicaragua.
All simulations are performed from January 1st, 2009 to December 31st, 2011, excluding the first month from the analysis due to spin-up time, covering two major ENSO events with extreme socio-economic consequences in Colombia (El Niño 2009and La Niña 2010, Hoyos et al. (2013, Comisión Económica para América Latina y el Caribe (2012)) as well as some ENSO-neutral periods (i.e. March-July 2009, March-June 2010 and May-July 2011 For this period, we find an average Oceanic Niño Index (ONI) value (based on the 3-month running mean of SST anomalies in the Niño 3.4 region [5°N-5°S, 120°-170°W]) of − 0.34 (as computed by the Climate Prediction Center: https:// origin. cpc. ncep. noaa. gov/ produ cts/ analy sis_ monit oring/ ensos tuff/ ONI_ v5. php), meaning that our simulation period (2009)(2010)(2011) is slightly biased towards more representative of weak La Niña/Neutral periods. It appears that Summer 2010-Spring 2011 was characterized by an Eastern-Pacific (EP) La Niña event of Medium grade (although a strong minimum trimestrial ONI of − 1.6) while Summer 2009-Spring 2010 was marked by a Central-Pacific (CP) El Niño event of Medium grade (with an amplitude of 1.7 °C) (Ren et al. 2018).
We used two-way nesting to allow interactions between inner and outer domain which is the recommended method for simulations lasting longer than a few days ). Both domains have 35 vertical levels. D01 has a horizontal resolution of 30 km. A finer 10 km resolution for D02 was chosen as a trade-off between computational time and high resolution of fine-scale climate events. Initial and boundary conditions for the simulations were retrieved from the 6-hourly, 0.75° × 0.75° gridded ERA-Interim reanalysis dataset provided by the European Center for Medium-Range Weather Forecasts (ECMWF, Dee et al. 2011).
To avoid large deviations between simulation and driving fields (Bowden et al. 2012) we use spectral nudging on both domains during the whole simulation run, which has been shown to improve simulated results especially in high-resolution downscaling studies (e.g. Chotamonsak et al. 2012;Bowden et al. 2012;Heikkilä et al. 2010). Additionally, Paul et al. (2016) found a reduction in precipitation bias, specifically in rainy regions. Nudging, which comprises relaxation techniques, maintains on the one hand large-scale features from the input data driving field and on the other hand, simulates small-scale features. To let the model develop its own mesoscale and synoptic structures in the surface-near troposphere, we decide to nudge zonal and meridional wind (U and V), temperature (T) and geopotential height (PH) only above the planetary boundary layer additionally to switching off nudging in the first ten levels from the bottom of the model. Similar to Posada-Marín et al. (2018) and Paul et al. (2016), we use a relaxation time of approximately one hour, which corresponds to a value of 0.0003 for U, V, T and PH nudging coefficient (g uf , g t , g ph ). The top wave number, which is the number of waves contained in the domain, is set to three in x and y direction, with three being the maximum one that is nudged, to capture ERA-Interim features with wavelengths of approximately 1200 km and upwards (the WRF domain size is about 3000 km × 2940 km in zonal and meridional directions, respectively). Nudging is exerted every 6 h, which coincides with the frequency of our ERA-Interim reanalysis data.
We also rely on previous WRF studies in Colombia or in nearby countries that tested several microphysics, cumulus, longwave radiation, and planetary boundary layer schemes (Arregocés et al. 2021;Núñez 2014;García 2014;Posada-Marín et al. 2018;González-Rojí et al. 2022) to optimize our CTRL configuration. However, the sensitivity of our results to an exhaustive panel of WRF schemes is beyond the scope of this study. In consequence, further physical parametrization schemes considered include the microphysics scheme of the WRF single moment 6-class (Hong and Lim 2006), the new version of the Rapid Radiative Transfer Model (RRTMG) for shortwave and longwave radiation (Iacono et al. 2008), the Noah land surface model (Tewari et al. 2004), the Monin-Obukhov similarity scheme for the surface-layer (Jiménez et al. 2012) and the Yonsei University (YSU) planetary boundary layer ) as  (Kain 2004). We use adaptive time-stepping to avoid the model run crashing when calculating too high windspeed on steep slopes.
We stress here that those model settings (simulation months, schemes, boundary conditions, spin-up, nudging) are similar to the ones used in other recent studies on the detection or attribution of regional and local climatic events, and this standard nudging setting strongly limits the overall model internal variability across simulations (Teklay et al. 2019;Wang et al. 2020;Glotfelty et al. 2021;Lu et al. 2021). Meanwhile, across the simulation period, we find no significant drift in the differences between simulations (< 5%, not shown) and we perform sensitivity experiments in Sect. 4.4 to albedo parametrization and physical schemes to explain responses and potential biases.

Experimental setup
In total, we conduct five main simulations, each one representing a different land cover and/or climate input, for both domains: (1) CTRL simulation: forced by ERA-Interim data for 2009-2011, with standard WRF-built-in land cover map (2) WITHOUT_LCC simulation: forced by ERA-Interim data for 2009-2011, preindustrial land use/land cover (potential natural vegetation) (3) WITHOUT_CC simulation: forced by "pseudo anthropogenically-unforced" ERA-Interim data for 2009-2011, with standard WRFbuilt-in land cover map (4) WITHOUT_CC_LCC simulation: forced by "pseudo anthropogenically-unforced" ERA-Interim data for 2009-2011, preindustrial land use/land cover (5) 100% DEFORESTATION Simulation: forced by ERA-Interim data for 2009-2011 with 100% replacement of forest by cropland and pasture Those simulations are performed with the same physical parameterization schemes and parameters. Land use categories in the baseline simulation (simulation 1) are based on the USGS land-cover classification (24 categories, see Fig. 2a). For Colombia, this means that it is largely covered by natural vegetation: evergreen broadleaf forest (EBF, Amazonas, Chocó) and savanna (eastern Llanos). Urban and agricultural land is located mostly between the central and the eastern Andean branch, and in the north at the Caribbean coast. According to the built-in USGS land cover classification, only a few pixels are classified as urban and builtup land (Bogotá: 3 grid points; Cartagena, Cali and Barranquilla: 1 grid point, respectively). The most prominent agricultural land use type is dryland cropland and pasture. The Caribbean coast is dominated by a cropland/woodland mosaic.
To retrieve LCC effects on climate, in simulation 2, we changed all categories representing urban and agricultural land to potential natural vegetation (EBF and grassland) classes from inspecting coarse preindustrial land cover maps suggested by Levavasseur et al. (2012) and Hengl et al. and 4 (given is the landcover fraction changed in percent) and selection for regional analysis (2018, simulation 2). Each grid cell consists of fractions of more than one land cover type. Therefore, we decide to replace the vegetation properties (minimum and maximum values for albedo, leaf area index, emissivity, roughness length, as well as green vegetation fraction, rooting depth and stomatal resistance) in the classes "Urban and Built-Up land", "Dry cropland and pasture" as well as "Cropland/Woodland mosaic" with the properties of the class "Evergreen Broadleaf Forest". The properties of the classes "Irrigated cropland and pasture", "Mixed dryland/irrigated cropland and pasture" and "Cropland/Grassland Mosaic" are replaced by the vegetation properties of the class "Grassland" (see Fig. 2b). Simulation 3 is forced by near-preindustrial (~ 1900) ERA-Interim data. 'Pseudo anthropogenically-unforced (PAU)' ERA-Interim data consisted of reconstructing historical climate as if there would be no anthropogenic global warming influence. To achieve this, a three-step approach was employed: (a) Calculate 30-year monthly mean values for the periods 1900-1929 (~ near-preindustrial) and 1980-2009 (~ presentday) of variables (VAR) needed to force WRF e.g. air temperature, surface pressure, the two wind components and sea surface temperature. 30-year monthly mean values were calculated from the global reanalysis data set ERA-20C, equally provided by the ECMWF.
(b) Calculate the difference between preindustrial and present-day values of VAR.
(c) Create for each input variable new ERA-Interim PAU data at each time step: The underlying assumption is that the regional historical trend for those variables is mostly due to anthropogenic global warming forcing. Simulation 4 combines the adjustments made in simulation 2 and 3. The last simulation (simulation 5) represents a complete deforestation (100DEF) scenario which is forced by ERA-Interim data (as in simulation 1 and 2) with dryland cropland and pasture replacing present-day EBF in the whole domain.

Statistical analysis and visualization
To identify the relative contribution of LCC and CC, the areaweighted mean values of variables over the 3 year-simulation period are calculated at the national level as well as for two focus regions, namely the southern Colombian Andes (76.0°W -74.8°W, 5.2°N -1.8°N) between the central and eastern cordilleras, and the Caribbean coast (76.5°W -74.8°W, 11°N -7.6°N) in northern Colombia. These two regions were chosen because historically they were particularly affected by LCC (see Fig. 2a, b). Consequently, they presumably exhibit the highest temperature and/or precipitation changes from LCC in the domain.
Simulation 2 to simulation 5 were compared against the CTRL simulation. Hence, the effects were calculated as shown in Table 1.
We calculate the statistical significance of each simulation pair and associated confidence intervals at 95% using a Mann-Whitney-Wilcoxon non-parametric test, commonly used in regional climate studies, without presuming the shape of a variable distribution.
Further analysis focuses on the ENSO phases. Based on the Multivariate ENSO Index (MEI, PSL 2020) and the Oceanic Niño Index (ONI, CPC 2020), the El Niño study period is defined from June 2009 to March 2010, and La Niña from July 2010 to May 2011.

Evaluation of model performance
Since WRF has already been shown to simulate temperatures and precipitation satisfactorily in our model domain (Posada-Marín et al. (2018), we used a similar model configuration and visually evaluated model performance by comparing annual temperature and precipitation to the available observational 30-years average  provided by the Colombian Institute of Hydrology, Meteorology and Environmental Studies (IDEAM 2015). Underlying per-pixel data of the Fig. 3c, d were not publicly available to calculate biases, nor the 2009-2011 average (IDEAM, pers.comm.). However, visually, WRF is representing spatial patterns of annual temperatures satisfactorily (Fig. 3c). Only in the northern lowlands and between the central and eastern Andes WRF simulates peaks higher in altitude (i.e. colder than reanalysis) and valleys lower in altitude on average (i.e. warmer than reanalysis), as expected (see e.g. Figure 11 lower panel of Posada-Marin et al. 2018). Precipitation is generally overestimated across the whole domain but follows the observed spatial patterns.

Regional impacts of CC and LCC on temperature
In our simulations, global warming (CC) since the beginning of the last century increases regional temperature in Colombia (+ 0.77 ± 0.02 °C, Fig. 4b) while historical LCC and 100DEF lead to a significant weak cooling (− 0.01 ± 0.01 °C and − 0.16 ± 0.01 °C, respectively; Fig. 4a, c). Generally, these effects are significantly lower at the Caribbean coast (+ 0.58 ± 0.02 °C for CC and − 0.04 ± 0.01 °C for historical LCC) than in the Andean region (+ 0.65 ± 0.02 °C for CC and − 0.07 ± 0.01 °C for historical LCC). The maps show that the cooling effect from historical LCC is mainly located in the Andean region over potentially historically deforested areas, where forest was converted to cropland and pasture (Fig. 4a). CC effects are significantly stronger at higher altitudes (+ 0.17 °C/km above 500 m at the national level and + 0.12 °C/km across an Andean transect, r > 0.6 for both, in elevation-dependent warming vs. height diagram, see Fig. 5a, b).

Regional impacts of CC and LCC on precipitation
Historical LCC and 100DEF simulations both reduce mean precipitation in Colombia by -0.08 ± 0.22 mm/day (− 0.7%, not significant p > 0.05). In contrast, precipitation significantly increases due to CC by 0.98 ± 0.30 mm/ day (+ 9%). Combined effects are marginally smaller (+ 0.91 ± 0.31 mm/day) with negligible synergistic effects.
In terms of precipitation response, the Caribbean Coast is more affected by historical LCC and 100DEF than the Andean region. In the coastal region, LCC decreases precipitation by − 0.67 ± 0.12 mm/day (− 7%), while 100DEF leads to an increase of 0.44 ± 0.13 mm/day (+ 5%). In contrast, the CC impact is by a factor of two higher in the Andes than at the coast and at the national level. LCC and 100DEF both decrease precipitation in the Andes by − 0.48 ± 0.10 mm/day (− 4%) and − 0.18 ± 0.08 mm/day (− 1.5%), respectively (see Fig. 6a-d). Deforested areas clearly show precipitation reduction (Fig. 6a, c), while CC and combined effects show opposite trends, except for the Caribbean coast. Precipitation increase is also significantly higher at higher elevations, albeit less significantly than temperature (r < 0.3): by + 0.17 mm/day/km above 500 m at the national level and by + 0.32 mm/day/km across the Andean transect (Fig. 5c, d).
In the Andean region, LCC and CC effects are slightly enhanced in magnitude during La Niña (by − 0.01 ± 0.01 °C -i.e. + 14%-and by + 0.04 ± 0.01 °C -i.e. + 6%-respectively, compared to the full simulation period, Fig. 7c) and dampened during El Niño (− 0.02 ± 0.01 °C -i.e. − 28%-and − 0.08 ± 0.01 °C -i.e. − 12%-, respectively). The CC effect on precipitation depends on the ENSO phases: it is weakened by El Niño (+ 0.62 ± 0.31 mm/ day) but enhanced by La Niña (+ 1.08 ± 0.42 mm/ day, Fig. 7d). Combined effects respond in the same pattern, although slightly weaker than the CC effect (+ 0.55 ± 0.29 mm/day and + 1.0 ± 0.4 mm/day, for El Niño and La Niña respectively). The Andean region reacts in a similar pattern, just about twice as strong in magnitude (Fig. 7f). At the coast, an opposite reaction can be observed: CC and combined effects are stronger during El Niño and weaker during La Niña (Fig. 7e). Here, historical LCC reduces precipitation, especially during La Niña (− 0.84 ± 0.13 mm/day). If the complete domain would be deforested and replaced by dryland cropland and pasture, it would result in a general decrease in precipitation which would hit the Andes the  Fig. 7f). On the contrary, the 100% deforestation case increases precipitation at the coast across all studied time periods, with its highest impact during La Niña (+ 0.69 ± 0.19 mm/day, Fig. 7e).

Regional impacts of historical LCC and 100DEF on surface energy balance
Changes in the surface energy balance following historical LCC show a significant increase in ΔSW ↑ of 2.03 W/ m 2 (p < 0.05, Fig. 8a), resulting in a relative cooling effect. This represents a contribution of roughly 35% to the total magnitude of change (i.e. ΔSW ↑ divided by Δ sum of all absolute flux changes).
A similar change in magnitude but with an opposite sign for ΔL (− 1.92 W/m 2 , p < 0.05) has a warming effect on national level, thereby balancing the cooling from ΔSW ↑ . Δ(H + G) decrease and contribute also to warming, albeit the effect is smaller. Together, the heat fluxes account for 50% of the total magnitude of change (same as above: Δ non-radiative fluxes divided by Δ sum of all absolute flux changes).
Downward radiative fluxes are slightly decreased and contribute around 15% to the overall change. However, on a local scale, changes in fluxes due to LCC effects are imbalanced, especially at the coast, where the dominant land cover was changed from EBF to a cropland/woodland mosaic ( ΔSW ↑: + 12.43 W/m 2 vs. ΔL : − 9.96 W/m 2 ; p < 0.05). In the Andean region, where EBF was converted to dryland cropland and pasture, these two fluxes are similarly imbalanced ( ΔSW ↑: + 6.80 W/m 2 vs ΔL : − 4.26 W/ m 2 , p < 0.05), pointing to a dominance of albedo-driven cooling.
The 100DEF scenario with the replacement of all EBF by cropland and pasture in the model domain (Fig. 8b) reveals a more drastic picture compared to historical LCC: a stronger increase in reflected shortwave radiation from the surface and a smaller reduction in latent heat, which is strongest on country level. The effect exists also in the focal regions, though weaker because these areas are already largely deforested in the control simulation.
Albedo on average significantly rises by 0.008 and 0.04 following historical LCC and 100DEF in Colombia, respectively. At the coast, the albedo increases by 0.045 due to historical LCC and by 0.01 under 100DEF. The Andean region shows an albedo increase of 0.02, respectively, after LCC and complete deforestation.

Discussion
This study explores the effect of both land cover change and global warming particularly during ENSO events, as well as their combined effects on surface air temperature, precipitation and surface energy fluxes, over Colombia, a climate modeling coldspot.
The performance of the WRF model in simulating general temperature and precipitation patterns in Colombia is satisfactory for the purpose of this research ( Fig. 3; Sect. 3.1). An apparent overestimation of precipitation is found comparing our WRF CTRL simulation and IDEAsM patterns (the national observational dataset, see Fig. 3). This can be explained because ERA-Interim already slightly overestimates precipitation in Colombia (Posada-Marín et al. 2018), our large-scale climatic input data to WRF. In addition, Jin et al. (2010) and Teklay et al. (2019) explored different model configurations with varying land surface models (LSM) in a similar context and found that precipitation is overestimated regardless of the LSM used. Another reason for the discrepancy might be the simulation period of only three years of the recent decade.

CC effect on temperature and precipitation
Interestingly, our results on temperature changes substantially contrast with the official climate near-term projections from the Colombian IDEAM in its Third National Communication on Climate change in Colombia (TCNCC, IDEAM et al. 2017). Historical surface air temperature increases in our study compare well with national averages published by IDEAM (+ 0.77 °C vs. + 0.8 °C), however, spatial warming patterns are very different (Fig. 4b, d). Our simulations show the highest temperature increases occurring above Andean cordilleras and higher elevation areas (> 1 °C, dark red in Figs. 4b, d and 5a, b), lower temperature increases in the northern and western areas and intermediate increases in the south and eastern parts (e.g. Llanos and Amazon Forest). This is in contradiction with IDEAM results where higher elevation areas warm less (1-1.5 °C) and the Llanos or low elevation areas more but more homogeneously (> 1.5 °C, Fig. 4b, d). In other means, we show a positive elevationdependent warming (EDW) while the IDEAM showed the opposite result (negative EDW). Confusingly, the IDEAM states in this very same TCNCC, that the "most significant increase in mean, minimum and maximum temperatures will be in the Andean region, especially in high elevation areas" (IDEAM et al. 2017). In the literature, there is growing recent evidence that the rate of warming is amplified with elevation (Pepin et al. 2015;Wang et al. 2016;Hock et al. 2019). However, in the tropics, there is a particular problem with observational data from mountain regions, which is both sparse and inhomogeneous (Pepin et al. 2015), and EDW in the Tropical Andes is controversially discussed in the literature ( Concerning precipitation changes, global warming intensifies the global hydrological cycle (Jia et al. 2019), i.e. the warming of surface air temperature features more evaporation of water into the atmosphere, thus more water content in the air, since warmer air can hold more water vapor. However, according to an average of four emission scenarios, IDEAM et al. (2017) expect precipitation to decrease over the Caribbean coast and the Colombian Amazonas (by 10-40% by 2071-2100 compared to the reference period  and to increase in the Andean region. Our results suggest an overall precipitation increase in Colombia in response to historical global warming, except for the Pacific coast. Our findings are comparable to Skansi et al. (2013), who analyzed station-based time series data from 1950-2010 in South America. Our findings show a drying along the Pacific coast, which is confirmed at least for the southern part by Skansi et al. (2013) but in contradiction to Carmona and Poveda (2014). In WRF, the decrease in Western Pacific Colombian coast precipitation in response to global warming can be explained by less tropospheric winds coming from the Pacific which usually advect moisture, driven by less ocean-continent thermal contrast (not shown).
Furthermore, elevation-dependent rainfall changes under global warming are poorly studied (Li et al. 2017) and we are only aware of one study in the tropics (Urrutia and Vuille 2009, tropical Andes). This latter study does not show consensual patterns of elevation-dependent rainfall changes under global warming, which depends on the latitude or the slope orientation. The authors found that "In general the projected changes in precipitation are not as pronounced as for temperature" and that precipitation is slightly higher in the global warming scenario below 4000 m.a.s.l. and slightly less beyond that elevation (Urrutia and Vuille 2009). Our results show a slightly higher rainfall increase under historical global warming in Colombian Andean cordilleras vs. other low-elevation areas (Figs. 5b and 6c, d) which is a consistent pattern with IDEAM near-term projections (IDEAM et al. 2017) and some mid-latitude patterns (Tibetan Plateau, Li et al. 2017) but differs from Urrutia and Vuille (2009).
Finally, our results suggest that CC effects are more dominant on national climate in Colombia than LCC effects while LCC affect rainfall substantially more than temperature at regional and sub-regional scales, which is consistent with most recent studies (Hock et al. 2019;Quesada et al. 2017aQuesada et al. , 2017b. Moreover, at the subregional level (Coast or Andean sub-region of Colombia), LCC effects are of the same magnitude on precipitation as CC effects, which urges to systematically account for land-use and land-cover dynamics in any departmental and regional mitigation and adaption to climate change plans.

LCC effect on precipitation and temperature
A drying over tropical deforested areas is reported throughout the literature (Perugini et al. 2017;Coe et al. 2017;Spracklen and Garcia-Carreras 2015;Lejeune et al. 2015;Saavedra et al. 2020;Sierra et al. 2021;Eiras-Barca et al. 2020) which is in good agreement with our results, although the magnitude of precipitation change differs. Trees usually have a deeper rooting depth and are therefore able to access deep-lying groundwater even in dry periods. Removing tropical forests and replacing them by crops and pastures results in reduced evaporative leaves and canopies, foliar density, surface roughness, stronger winds and an increased albedo. Consequently, evapotranspiration and precipitation (due to fewer clouds) are expected to diminish in tropical agricultural areas. Furthermore, modeled (Sampaio et al. 2007;Nobre et al. 1991;Zemp et al. 2017) changes in precipitation show that drying is even more pronounced during the dry seasons and in low latitude areas.
Reviews on GCM and RCM modeling studies with a focus on Amazon deforestation report a precipitation decrease (Spracklen and Garcia-Carreras 2015;Lejeune et al. 2015;Sampaio et al. 2007). Lejeune et al. (2015) report a reduction of ~ − 0.7 mm/day (− 5 to − 10%) across 28 GCM and RCM experiments, in agreement with Spracklen and Garcia-Carreras (2015), who report a decrease of 16.5 ± 13% across 44 simulations of complete Amazon deforestation. This is similar to our findings of a 4% decrease in precipitation (− 184 mm/year or − 0.5 mm/day in case of total deforestation over the modeling domain) but more pronounced in magnitude.
Though the comparison with other regional WRF studies (Laux et al. 2017;Paul et al. 2016;Takahashi et al. 2017;Quesada et al. 2017a) shows discrepancies, due to the spatially specific influence of climatic phenomena, they agree that deforestation leads to decreasing precipitation over large-scale deforested areas.
Our simulations show that tropical deforestation causes a slight cooling at country level and more enhanced cooling over deforested grid-cells. Such cooling after deforestation simulated by default WRF parametrizations is at odds against all observations and climate simulations about the tropical deforestation impacts on temperature (Perugini et al. 2017;Lawrence and Vandecar 2015;Jia et al. 2019). The IPCC Special Report on Climate Change and Land (Jia et al. 2019) concludes that there is high confidence that large-scale deforestation leads to mean biophysical warming (+ 0.61 ± 0.48 °C) over the entire tropics and especially over deforested areas. After entire Amazon deforestation, GCMs and RCMs reviewed by Lejeune et al. (2015) simulated a temperature increase by 0.8 °C on average, whereas WRF in our study simulates a cooling of − 0.21 ± 0.02 °C over Colombian deforested grid-cells in the 100DEF scenario. Perugini et al. (2017) reviewed eleven modeling studies in the tropics and found a mean temperature increase of 0.68 °C (forest cropland/grassland) which is in agreement with three tropical observational studies they investigated (+ 0.41 °C on average). Temperature changes derived from satellite observations (+ 1.06 °C in Alkama and Cescatti (2016), + 1.34 °C in Duveiller et al. (2018)) or from forest change data (+ 0.38 °C in Prevedello et al. (2019) confirm the warming effect of deforestation in the tropics. However, caution is needed as observational studies do not capture nonlocal biophysical impacts from deforestation (Winckler et al. 2017).
Only a few modeling studies found a cooling as a result of tropical deforestation. For instance, Lejeune et al. (2015) report a mean increase in surface air temperatures of + 1.2 °C across 28 GCM and RCM studies conducted for the Amazon basin with only 3 studies showing cooling out of 28 in total. Robertson (2019), using the HadGEM2-ES ESM, obtained a cooling similar to our study (− 0.08 °C vs. − 0.01 °C) and showed with a decomposition of the surface energy balance that albedo and evapotranspiration responses were strongly biased (overestimated and underestimated, respectively). We suspected similar reasons being responsible for WRF incorrectly simulating cooling after tropical deforestation and thus conducted in the next section a decomposition of our results as well as sensitivity experiments modifying albedo parametrization and land surface models.

LCC and CC interactions with ENSO
Studies investigating the interactions between LCC and ENSO in tropical South America are rather rare. The consensus, however, is that LCC enhances the ENSO signal (Nobre et al. 2009;Beltrán-Przekurat et al. 2012;Bush et al. 2017;Zhang et al. 2009). Tölle et al. (2017), for example, used the RCM COSMO-CLM and investigated consequences of abrupt tropical deforestation during the ENSO phases in South East Asia and concluded that they "can amplify the impact of the natural mode ENSO". Another study by Chapman et al. (2020) on the island of Borneo (Indonesia) confirms that deforestation under El Niño conditions resulted in warmer and drier climate than under normal conditions. Syktus et al. (2007), whose study domain was Australia, came to the same conclusion. Bush et al. (2017) suggest that the ENSO signal has been even amplifying since humans started cultivating agricultural goods, which was buffered before by natural vegetation. The buffering effect has been confirmed by Meijide et al. (2018). However, our results suggest that it depends on the region, land cover type, ENSO phase and variable whether the LCC effect is amplified or not, i.e. there is no or little influence of LCC and CC on temperature during ENSO, whereas precipitation changes are more affected, especially during La Niña and at the Coast. Beltrán-Przekurat et al. (2012) concluded from their results that dry periods are enhanced by the occurrence of ENSO, but that the type of land cover conversion also plays an important role for the sign of change.

Decomposition of energy balance and model sensitivity to albedo parametrization and land surface models (LSM)
At a regional scale, the energy balance is deeply modified by forest losses. In the tropics, a deforestation-driven warming caused by decreased latent heat flux is expected to outbalance the albedo-driven cooling. When decomposing the energy balance into its single components and comparing it to observed values provided by Duveiller et al. (2018), WRF with the coupled Unified Noah LSM (Tewari et al. 2004), which is used here, overestimates the change in reflected shortwave radiation from the surface by a factor of 3.5 over converted grid-cells (EBF → crop/grass, see Table 2). At the same time, the magnitude in latent heat change reduction is underestimated by 41% in WRF. WRF simulates a latent heat flux reduction of 9.6 W/m 2 , while a reduction of 20 to 30 W/m 2 would be more realistic according to Duveiller et al. (2018), Nobre et al. (1991) and Snyder et al. (2004) when forests are converted to cropland and/or pasture. Furthermore, sensible heat is decreasing by -8.4 W/m 2 in WRF following deforestation, instead of increasing by + 5 W/ m 2 to + 20 W/m 2 as suggested in the literature (Duveiller et al. 2018;Nobre et al. 1991;Sampaio et al. 2007;Llopart et al. 2018;Takahashi et al. 2017;Coe et al. 2017;Oliveira et al. 2019). In summary, the most important reason for the deforestation-temperature bias in WRF seems to be the high overestimation of shortwave radiation reflected from the surface after deforestation, while biases in the latent/sensible heat partitioning tend to compensate each other. This is a coherent finding compared to Sierra et al. (2021) who mentioned the WRF overestimation of the reflected shortwave radiation and the strong underestimation of the sensible heat flux and evaporation/latent heat after cooling in response to forest loss in the tropical Amazon-Andes transition region. Indeed, Robertson (2019) concluded, that "to reduce this error the difference between the PFT's albedo parameters should be reduced". Following this recommendation, we adjusted the land use category's albedo parameters and ran new WRF simulations. The default albedo ranges from 0.17-0.23 for dryland cropland and pasture, from 0.19-0.23 for grassland and has a constant value of 0.12 for EBF. Observed values for albedo indicate that 0.13 to 0.19 for crops, 0.13 to 0.23 for grassland and 0.11 to 0.15 for tropical EBF are more appropriate (Cescatti et al. 2012, Hänchen 2017Prentice et al. 1992;Giambelluca et al. 1997). After running simulation 1 and 2 again with these adjusted values (e.g. using minimum albedo of 0.11 and a maximum albedo of 0.15 for EBF), ΔSW ↑ is now only slightly overestimated (33%). However, the underestimation of ΔL is amplified (by 83%, see Table 2).
Additional runs were conducted to explore the model's sensitivity to the LSM used. Table 2 shows that the cooling over deforested pixels is even more enhanced when applying the Rapid Update Cycle (RUC) LSM (− 0.13 °C, Benjamin et al. 2004) and the 5-Layer thermal diffusion scheme (− 0.13 °C, Dudhia 1996). Only the Pleim-Xiu LSM (Pleim and Xiu 1995) simulates a warming of 0.06 °C over deforested grid-cells. Furthermore, a decomposition of the alternative LSM's energy balance shows also changes in radiative and non-radiative fluxes that are in stark contrast to observed changes (Duveiller et al. 2018, see Table 2). Thus, we investigated several sources of deforestation-induced temperature biases in WRF, which were found in recent tropical studies performed with WRF in the default configuration (Glotfelty et al. 2021;Sierra et al. 2021).
Finally, we stress here that we are confident with WRF responses to CC (Sect. 4.2) as well as precipitation responses to LCC: (1) vegetation model properties are the same between CC simulations and climatic patterns are well represented by WRF in this region -see Sect. 3.1 and Posada-Marín et al. (2018); (2) precipitation responses after LCC and tropical deforestation (100DEF) are consistent with the observations in this tropical region (see Sect. 4.3), (3) the biases tend to cancel out: the underestimation of evapotranspiration loss after deforestation has an opposite effect with respect to the overestimated cooling for precipitation: the former tends to underestimate the magnitude of the reduced precipitation response (more evapotranspiration leads to more clouds and in turn more precipitation) while the latter tends to overestimate it (cool air tends to hold less water thus less precipitation, Mahmood et al. (2014); (4) across the different WRF land cover schemes tested, precipitation responses are not substantially affected and are systematically negative after deforestation (not shown).

Conclusion
Colombian climate regional modelling is challenging in several respects: it has a highly uneven and mountainous terrain highlighting the need for high-resolution modelling, it is a regional climate modelling coldspot (very few studies), it is highly vulnerable territory to climate change and biodiversity loss and thus, it calls for an urgent need of hydroclimatic assessments (Poveda et al. 2011). In this study, we quantified and disentangled the impact of LCC reduces precipitation everywhere over deforested areas, but especially in the Caribbean coast in the north. Conversely, CC increases precipitation, specifically in the Andean region. LCC effects are found to be strong drivers of sub-regional precipitation changes: dampening CC effects magnitude by 72% and 24% on average in Coastal Caribbean and Andean regions, respectively. Combining historical LCC and CC, temperature increases by 0.76 °C, especially in high altitudes. Conflicting with most other modeling studies and observations, WRF simulates a temperature decrease following historical LCC and 100DEF, which is intensified over deforested areas.
Our results suggest that CC effects are more dominant on regional climate in Colombia than LCC effects but both effects are comparable in magnitude over deforested areas or at the subregional level. La Niña tends to enhance the studied effects across the country, whereas El Niño weakens them, except at the coast where the pattern is mostly inverted. In the Andes, La Niña tends to intensify LCC and CC effects roughly to a factor of two compared to El Niño impacts in both cases.
Our simulations indicate a serious bias in WRF concerning an incorrect simulation of cooling after tropical deforestation, recently found as well in other tropical regions (Glotfelty et al. 2021;Sierra et al. 2021). We further investigate this apparent bias because of its great importance for climate studies using WRF with default parametrization in tropical contexts and probably beyond. We show here that WRF coupled with the Unified Noah LSM, RUC LSM and the 5-layer thermal diffusion LSM overestimates reflected surface shortwave radiation response and underestimates evapotranspiration-driven warming in the tropics which leads us to the assumption that vegetation properties are not appropriately parametrized for tropical regions in all investigated LSMs in WRF. We therefore recommend using default vegetation properties of WRF with caution for temperature assessments. Improved model representation of vegetation properties and more validation data in tropical areas are needed to provide more robust and confident climate projections at regional-to-local scale. We stress here that, although our study has limits concerning the limited sensitivity experiments across physical schemes or imprecise representation of land-cover characteristics, the model settings are tuned to limit internal model variability maximizing the signals response, following similar WRF climatic studies (see Sects. 2.2 and 2.4). Meanwhile, we provide statistical significance to ensure robust messages throughout this first-of-its-kind attribution study at high-resolution (10 km) analyzing global warming and land-cover changes impacts on climate in this region. Although our simulation period is limited (2009)(2010)(2011), our results on the whole period are likely to represent the simulated impact of historical land cover changes and global warming on Colombian climate during a climatological period, albeit slightly skewed towards La Niña.
To avoid obvious temperature biases in WRF simulations, in-depth analysis and further development of land cover classification and the representation of vegetation properties used by WRF is necessary. This includes a more realistic and up-to-date representation of tropical land cover data. USGS land use classification in WRF standard simulations dates back to the early 90s. Since then, tropical deforestation and land use rose tremendously. A thorough review of vegetation properties parametrization is necessary to assure a good representation of tropical land-atmosphere feedbacks. We implemented alternative albedo values for Colombia's most important land use categories and tested its sensitivity, which showed a substantial improvement of simulated reflected shortwave radiation. However, parametrizations of e.g. roughness length and stomatal resistance, important parameters which drive evapotranspiration, were not investigated, potentially explaining that simulated latent heat fluxes were even more biased. To understand the role of the ENSO phases better, longer simulation studies are needed, in order to have a more appropriate climatological period and an average of more ENSO events. In addition, a follow-up study could investigate climate extreme indices like maximum daily temperature and maximum daily precipitation to provide evidence for LCC and CC impacts on extreme weather events. Although land use data has been used for regional climate simulations in the Neotropics (e.g. Sierra et al. 2021;Saavedra et al. 2020) and although a more recent update land use data already exists (years 2000s, Eva et al. 2004), a newer high-resolution land-cover dataset for Colombia spanning a climatological period would represent a great input for more robust climate and hydrological assessment in Colombia. Finally, it should be noted that the official climatic projections by the TCNCC (IDEAM et al. 2017) contains methodological shortcomings, inconsistencies and need urgent reassessment by the national communities of meteorology and climatology (Arias et al. 2022).
Funding Open Access funding provided by Colombia Consortium. We deeply acknowledge the services of the High Performance Computing center hosted at Universidad del Rosario (Colombia) and Advanced Computer Laboratory for Research as well as the useful technical help provided (https:// www. urosa rio. edu. co/ Labor atorio-Compu tacion-Avanz ada-Inves tigac ion/ Infra estru ctura/"). A.Manciu acknowledges 1 3 TUM university and the Hanns-Seidel-Stiftung for the Mobility Fellowship as well as the Universidad del Rosario for supervision. B.Quesada acknowledges Climat AmSud program (project 21-CLI-MAT-13), Colombian-french researcher association COLIFRI in the framework of the French FSPI fund on Renewable Energy ecosystem in Orinoco Region project and the Starting Grant of Universidad del Rosario for funding this work.
Data availability For this research, we did not use any new data for simulating the studied effects. The ERA-Interim reanalysis dataset from Dee et al. (2011), with which we forced the model, and the ERA-20C dataset from Dee et al. (2016), that was used to create the forcing data without Global Warming influence, is available for download from the European Center for Medium-Range Weather Forecasts (ECMWF, available at https:// www. ecmwf. int/ en/ forec asts/ datas ets/ browse-reana lysis-datas ets) upon creation of a cost-free user account. The newly simulated and processed data in this research is available in the Center for Open Science repository under the following link: https:// osf. io/ pcrk5/? view_ only= 5f80e 18935 1b484 5bd9d f72b2 8ed87 fe.

Conflict of interest
The authors declare no conficts of interest or competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.