An Assessment of Land–Atmosphere Interactions over South America Using Satellites, Reanalysis, and Two Global Climate Models

: In South America, land–atmosphere interactions have an important impact on climate, particularly the regional hydrological cycle, but detailed evaluation of these processes in global climate models has been limited. Focusing on the satellite-era period of 2003–14, we assess land–atmosphere interactions on annual to seasonal time scales over South America in satellite products, a novel reanalysis (ERA5-Land), and two global climate models: the Brazilian Global Atmospheric Model version 1.2 (BAM-1.2) and the U.K. Hadley Centre Global Environment Model version 3 (HadGEM3). We identify key features of South American land–atmosphere interactions represented in satellite and model datasets, including seasonal variation in coupling strength, large-scale spatial variation in the sensitivity of evapotranspiration to surface moisture, and a dipole in evaporative regime across the continent. Differences between products are also identiﬁed,withERA5-Land,HadGEM3, andBAM-1.2showingoppositeinteractions tosatellitesoverpartsoftheAmazon and the Cerrado and stronger land–atmosphere coupling along the North Atlantic coast. Where models and satellites disagree on the strength and direction of land–atmosphere interactions, precipitation biases and misrepresentation of processes controlling surface soil moisture are implicated as likely drivers. These results show where improvement of model processes could reduce uncertainty in the modeled climate response to land-use change, and highlight where model biases could unrealistically amplify drying or wetting trends in future climate projections. Finally, HadGEM3 and BAM-1.2 are consistent with the median response of an ensemble of nine CMIP6 models, showing they are broadly representative of the latest generation of climate models.


Introduction
Global climate models offer a way to predict how increasing anthropogenic emissions and land-use change will impact climate. However, before models can be reliably used for future projections, it is necessary to validate their performance in the present day, including their ability to represent key physical processes in the climate system (Flato et al. 2013). Such process-based model evaluation has become increasingly important as models have gained in complexity, with novel evaluation approaches over the past decade helping to drive improvements in model development, and to assess the credibility of future climate projections (Eyring et al. 2016a;Duveiller et al. 2018;Eyring et al. 2019Eyring et al. , 2020Fasullo 2020).
In South America, interactions between the land surface and the atmosphere are particularly important for climate, and thus need to be accurately represented in climate models. Studies integrating remote sensing and reanalysis datasets have highlighted the importance of soil moisture (SM), with variation in SM shown to drive spatiotemporal variation in evapotranspiration (ET), surface temperature (T), and precipitation (P) (Spennemann and Saulo 2015;Bedoya-Soto et al. 2018). SM impacts local climate by affecting surface albedo and the partitioning of heat loss to the atmosphere into sensible and latent heat fluxes (Bowen ratio) (Avissar 1995;Eltahir 1998;Koster et al. 2004;Seneviratne et al. 2010;Dirmeyer 2011). SM may also influence climate remotely, with P over some areas of central and southern South America dependent on water that Denotes content that is immediately available upon publication as open access.
Variation in vegetation cover is also important for regulating land-atmosphere interactions in South America (Spracklen et al. 2018). Extensive land-use change, both inside and outside Amazonia, has already impacted climate on local and regional scales (Baidya Roy and Avissar 2002;Silvério et al. 2015;Salazar et al. 2015Salazar et al. , 2016Baker and Spracklen 2019). Deforestation has been linked to lengthening dry seasons in the southern Amazon and in the seasonally dry Cerrado biome (Costa and Pires 2010;Fu et al. 2013). Many models struggle to represent the complex biophysical changes that accompany land-use change, disagreeing on the direction of ET, T, and P responses over South America (Duveiller et al. 2018;Boysen et al. 2020). These studies demonstrate that accurate representation of landatmosphere interactions in models is necessary to quantify the climate response to future land-use change.
Models can contribute to mechanistic understanding of land-atmosphere interactions, as model fields may be available for variables that cannot be observed directly, and model experiments can be used to isolate the importance of different variables and processes. Studies based on regional climate models have identified the Amazon, eastern Brazil, and southeastern South America, including La Plata basin, as key regions where interactions between land and atmosphere have an important impact on the local climate (Sörensson and Menéndez 2011;Llopart et al. 2014;Ruscica et al. 2015;Menéndez et al. 2016;Spennemann et al. 2018;Menéndez et al. 2019). Over eastern Brazil and La Plata Basin, where T is highly variable, a study based on two regional circulation models showed that models with fully interactive land-atmosphere feedbacks simulated greater T seasonality and higher interannual T variation, than models where climatological SM was prescribed . Similarly, a comparison of models run with and without land-atmosphere coupling over the Amazon showed that models with prescribed SM underestimated the dry season T response to El Niño over the eastern basin (Levine et al. 2019). In southeastern South America, land-atmosphere coupling strength varies seasonally, peaking in austral summer when SM and ET variability are sufficiently high to influence the atmospheric response (Sörensson and Menéndez 2011;Ruscica et al. 2015Ruscica et al. , 2016. Detailed evaluation of variables in the surface-to-atmosphere moisture feedback pathway showed that SM influences P by modulating the surface heat flux, thus coupling SM to the vertical gradient in moist static energy (MSE), where higher (lower) MSE is usually associated with lower (higher) atmospheric stability, favoring (suppressing) P (Eltahir 1998;Ruscica et al. 2015). These studies illustrate how models can enhance our understanding of land-atmosphere interactions, though they depend on models providing an accurate representation of the relevant processes, which may not be guaranteed.
Evaluation of land-atmosphere interactions over South America in global climate models has so far been limited, although evidence suggests representation is variable at best. On short time scales, models show a tendency of enhanced afternoon rainfall over areas with high SM, which is opposite to the observed response (Taylor et al. 2012). Furthermore, a metaanalysis of studies modeling the climate impacts of complete Amazon deforestation found simulated precipitation changes ranged from 15% to 255%, illustrating high uncertainty in the atmospheric response to changes in the land surface (Spracklen and Garcia-Carreras 2015). Another study evaluating landatmosphere interactions over the Amazon showed 20 out of 38 Coupled Model Intercomparison Project (CMIP) phase 5 (CMIP5) models analyzed, and 10 out of 26 CMIP6 models, simulated a water-limited evaporative regime over the Amazon, rather than a radiation-limited regime as indicated by observations (Baker et al. 2020, manuscript submitted to Environ. Res. Lett.). The authors showed that the direction of interactions influenced whether future changes in climate were likely to be amplified or moderated. These examples illustrate the importance of evaluating model representation of land-atmosphere interactions, and identifying model processes to target for improvement.
A key challenge facing researchers wishing to study landatmosphere interactions over South America, and evaluate their representation in climate models, is the limited availability of high-quality observational datasets, particularly over the Amazon. Ground-based climate data are scarce in the Amazon and therefore datasets based on interpolated station data (e.g., Harris et al. 2014) are less likely to be reliable in this region. Satellite products can provide useful estimates of some variables, such as T, SM, and P. However, SM estimates are unavailable over most of the Amazon, as closed-canopy vegetation obscures the land surface, thus interfering with satellite retrievals of SM (Dorigo et al. 2017). ET is even more challenging to quantify over large spatial scales as retrieval from space is particularly complex, in situ data are rare, and comparison studies have revealed important discrepancies between products over the Amazon Baker et al. 2021). Where satellite data are unavailable or unreliable, reanalysis products can potentially provide useful estimates of climate variables. Reanalysis fields are not pure observations, but the result of assimilating observations into a numerical model to generate a best estimate of the climate state on a homogenous grid. Reanalyses have been widely used in model evaluation studies (e.g., Flato et al. 2013, and references therein), though they themselves require validation, where possible, to ensure reliability.
There is a clear need to evaluate the representation of landatmosphere interactions in climate models in order to determine where model improvements are needed, and ultimately to reduce the uncertainty in future climate projections. In this study, we selected two global climate models representing different stages of model development for a detailed assessment of land-atmosphere interactions over South America: the U.K. Hadley Centre Global Environment Model version 3 (HadGEM3), which has been developed over several decades and has been included in multiple CMIP phases (Murphy 1995;Johns et al. 1997;Pope et al. 2000;Collins et al. 2001;Johns et al. 2006;Martin et al. 2011;Williams et al. 2018;Andrews et al. 2020), and the Brazilian Global Atmospheric Model version 1.2 (BAM-1.2), which was developed comparatively recently (Figueroa et al. 2016;Coelho et al. 2021). Furthermore, we wanted to understand how representation of land-atmosphere interactions in these two models compared with wider ensemble of CMIP6 models, to assess the generalizability of our findings. Limited in situ observations and uncertainties in satellite retrievals and reanalysis datasets make such evaluations challenging, suggesting we need to understand the behavior of land-atmosphere processes across all types of datasets before models can be critically assessed. Therefore, the aims of this study were 1) to quantify land-atmosphere interactions over South America using satellite products, a novel reanalysis, and the global climate models HadGEM3 and BAM-1.2; 2) to compare results from all products, and identify similarities and differences between them; and 3) to compare results from HadGEM3 and BAM-1.2 to results from an ensemble of nine CMIP6 models.

Methods
We compiled a suite of land-atmosphere coupling metrics based on previous literature to evaluate the surface soil moisture-evapotranspiration-precipitation (SM-ET-P) pathway in satellite products, a reanalysis and two global climate models. We focused on evaluating land-atmosphere interactions over South America, but the metrics applied here can be adapted to evaluate other processes, over any region of interest across the world. A link to the repository containing the fully commented diagnostic scripts (written in the programming language Python) is provided at the end of the paper.

a. Land-atmosphere coupling metrics
The metrics used here measure the strength of coupling between the land and the atmosphere, where ''coupling'' refers to the sensitivity of the atmospheric state to variation in the land surface (land-atmosphere coupling), or vice versa (atmosphere-land coupling) (Seneviratne et al. 2010;Levine et al. 2016). We specifically aimed to diagnose and understand the surface-to-atmosphere moisture transfer pathway (i.e., SM-ET-P). Given the challenges associated with measuring SM and ET over large parts of South America, we selected five metrics that could be used to quantify and visualize this pathway using complementary approaches and data requirements. The first three metrics specifically relate to the SM-ET segment of the pathway, offering alternative methods to quantify SM-ET coupling. These include the terrestrial coupling index (TCI), which measures the sensitivity of ET to variation in SM (Guo et al. 2006;Dirmeyer 2011); the temperature-evapotranspiration metric (T_ET), which uses relationships between T and ET to infer SM-ET coupling (Seneviratne et al. 2006(Seneviratne et al. , 2010; and Zeng's gamma (G), a quantity that uses P as a proxy for SM to estimate SM-ET relationships (Zeng et al. 2010). In addition to these we used the Betts' approach, which calculates and displays associations between any two variables averaged over a specified domain of interest (Betts 2004), and the two-legged metric, which traces the full SM-ET-P pathway from the surface to the atmosphere in a mechanistic way (Guo et al. 2006;Dirmeyer 2006;Dirmeyer et al. 2014). A full description of each metric is provided in the online supplemental material.

b. Satellite and reanalysis products
The satellite and reanalysis products used in this study are listed in Table 1. Satellite products were obtained for the common time period of 2003-14 and all datasets were scaled to monthly resolution. Satellite retrievals of SM were obtained from the European Space Agency Climate Change Initiative (ESA CCI) Combined product (Liu et al. 2012), which merges information from active and passive satellite sensors. Satellite SM products typically offer information on the moisture content of the top 5 cm of the soil profile (Ray and Jacobs 2007), and thus provide a measure of surface SM only. We used satellite ET from the Moderate Resolution Imaging Spectroradiometer (MODIS) MOD16 ET product (Mu et al. 2007(Mu et al. , 2011. Precipitation data came from the Tropical Rainfall Measuring Mission (TRMM) Multi-Satellite Precipitation Analysis (TMPA) 3B43 monthly version 7 product (Huffman et al. 2007). Monthly surface T data were obtained from the Level 3 AIRX3STM version 6 product, based on retrievals from the Atmospheric Infrared Sounder (AIRS) instrument on board the Aqua satellite (Aumann et al. 2003;Tian et al. 2017). We took the mean of measurements from ascending [local overpass time 1330 local time (LT)] and descending (local overpass time 0130 LT) orbits (variable names SurfSkinTemp_A and SurfSkinTemp_D), which was previously shown to provide a reliable estimate of surface T (Susskind et al. 2019). Shortwave incoming solar radiation (RDN) data were retrieved from the 0.258 3 0.258 CLARA-A1 dataset (Karlsson et al. 2013). Previous work has shown inconsistencies between CLARA-A1 and other satellite radiation products (e.g., Loew et al. 2016;Sörensson and Ruscica 2018), and thus comparisons between models and satellite radiation data should be interpreted with some caution. However, as radiation data were not used for any of the metric calculations, radiation uncertainties are expected to have a relatively small impact on the interpretation of results in this study.
We obtained monthly output from ERA5-Land, a stateof-the-art global reanalysis from the European Centre for Medium-Range Weather Forecasts (Hersbach et al. 2020). ERA5-Land is a downscaled ERA5 product reported to have improved accuracy for land applications. To produce ERA5-Land, remote sensing and in situ observations were assimilated into the Integrated Forecasting System Cy4lr2 model to quantify the climate state on a four-dimensional grid. The land surface scheme is the Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (HTESSEL; Balsamo et al. 2015). As a reanalysis, ERA5-Land theoretically occupies an intermediate position between observations and pure simulation, and may provide insights over regions where observational data are unavailable. We downloaded ERA5-Land ET, P, T, and RDN as single-layer datasets. For surface SM we used data from the uppermost layer of the soil column (0-7 cm) to allow close comparison with satellite SM.

c. Models
We used Atmospheric Model Intercomparison Project (AMIP) simulations following the CMIP6 protocol. CMIP6 AMIP runs cover the historical period only , and are prescribed with observed sea surface temperatures, seaice concentrations, and all volcanic, solar, and anthropogenic forcings, including atmospheric CO 2 concentrations (Eyring et al. 2016b). Our analysis focused on two models, the U.K. Hadley Centre Global Environment Model version 3 (HadGEM3), and the Brazilian Global Atmospheric Model version 1.2 (BAM-1.2), but we later compare these results with AMIP runs from nine additional CMIP6 models, where the necessary data were available. The model variables needed to replicate the analyses shown in this study are listed in Table S1 in the online supplemental material and include ET, P, surface T, RDN, and surface-layer (0-10 cm) SM. Detailed information on the processing of model output is provided in the supplemental material.
AMIP runs from HadGEM3-GC31-MM (HadGEM3), the medium-resolution, global coupled model from the U.K. Met Office (Met Office Hadley Centre 2019; Williams et al. 2018;Andrews et al. 2020), were downloaded from the CMIP6 Earth System Grid Federation (ESGF) archives (https://esgf-index1.ceda.ac.uk/search/cmip6-ceda/) at monthly resolution (surface SM data only were downloaded at daily resolution and converted to monthly means, due to monthly output being unavailable for this variable). MM (medium resolution in atmosphere and ocean) simulations (N216) have a horizontal resolution equivalent to approximately 60 km in the midlatitudes (Roberts et al. 2019). HadGEM3 uses the Unified Model global atmosphere (GA7) and the JULES global land (GL7.1) schemes (Table S2). Further details on the configuration of HadGEM3 is provided by Walters et al. (2019).
In addition, we obtained monthly AMIP-type runs from BAM-1.2 (Coelho et al. 2021 (Figueroa et al. 2016;Guimarães et al. 2020). In this study, BAM-1.2 was run with triangular 126wave truncation (TQ126), corresponding to a horizontal grid of approximately 1.08 3 1.08 at the equator, 42 sigma vertical levels (32 levels in the troposphere and 10 in the stratosphere), and the model top was at 2 hPa. BAM-1.2 uses the land surface scheme IBIS-CPTEC (Kubota 2012). A full description of BAM-1.2 and its physical components is provided by Coelho et al. (2021). A summary of the characteristics of ERA5-Land, HadGEM3, and BAM-1.2 is presented in Table S2.
To compare land-atmosphere interactions in the U.K. and Brazil climate models with the latest generation of global models, we downloaded monthly AMIP simulations from a further nine CMIP6 models (Table S3) for the period 2003-14. We selected models that provided output for surface SM, ET, P, and T, excluding models from the same modeling centers as HadGEM3 or BAM-1.2.

d. Data analysis
We focused on assessing land-atmosphere moisture interactions over South America and four key subregions (Fig. 1). These included the Amazon and La Plata basins, which represent the two largest watersheds in South America, and the Brazilian Cerrado and Caatinga biomes, seasonally dry ecoregions of woody savanna and semiarid vegetation, respectively. The shapefile for the Amazon was obtained from the Observation Service SO HYBAM (http://www.ore-hybam.org), La Plata from Centro de Investigaciones del Mar y la Atmósfera (CIMA, http://www.cima.fcen.uba.ar/ClarisLPB/), and shapefiles for the Cerrado and Caatinga biomes were downloaded from Map for Environment (https://mapforenvironment.org/ layer/info/751/#4.3/-13.19/-48.45). We analyzed land-atmosphere interactions in satellite products, ERA5-Land and global climate model output using the five metrics outlined in section 2a. All datasets were regridded to the same resolution as BAM-1.2, using an areaweighted regridding approach. For HadGEM3 and BAM-1.2, which each had four ensemble members, we computed landatmosphere interactions for each member individually, and then calculated the mean result for each metric. Landatmosphere interactions varied slightly between members, and metrics computed using ensemble mean data were not necessarily representative of the mean of the results from individual members (Figs. S1 and S2). This is because computing the ensemble mean prior to applying the landatmosphere coupling metrics averages out spatiotemporal variability in the surface, flux and atmospheric variables, and thus land-atmosphere interactions estimated from the ensemble mean tend to be slightly muted. Therefore, all of the metric results presented here represent the mean of metrics computed from individual members, unless otherwise stated.
To gain a broad overview of coupling at multiple time scales, metrics were calculated at the annual scale using data from all months (i.e., a time series of 12 3 12 5 144 data points), and at the seasonal scale using data from 3-monthly periods [December-February (DJF), March-May (MAM), June-August (JJA), and September-November (SON), i.e., a time series of 3 3 12 5 36 data points), for the 12-yr period 2003-14 (the common period of available data). Since most of the climate data analyzed here were not normally distributed (Fig. S3), Spearman's rank correlation coefficients were computed for all metrics, with the exception of the Betts' approach, which used Pearson's correlation coefficients derived from a scatterplot of the two investigated variables (though correlations from Spearman's rank were found to be very similar). Finally, while it would have been preferable to analyze landatmosphere interactions over a longer time period, we were limited by the availability of appropriate satellite reference data, and the fact that CMIP6 AMIP simulations finish in the year 2014.
Finally, the nine CMIP6 models (Table S3) were evaluated at the annual time scale to provide a wider context for the rest of the results. The TCI, T_ET, and Zeng's G metrics were calculated for each model member in turn (the number of members per model is listed in Table S3). The mean of the results from individual members was calculated for each model, and the median metric value over each of the four study regions was extracted. The distribution of median values across the nine CMIP6 models was compared with median values from satellite products, ERA5-Land, and the four members of HadGEM3 and BAM-1.2.

a. Controls on surface moisture fluxes
A key step in the hydrological pathway is the transfer of moisture from the land to the atmosphere via ET. Following the terminology of Seneviratne et al. (2010), we refer to areas where variation in ET is strongly related to surface SM as having strong land-atmosphere coupling. When land-atmosphere interactions are dominated by an atmospheric influence, such as a radiation control on ET, the coupling is in the opposite direction (atmosphere-land coupling). To evaluate representation of land-atmosphere interactions in HadGEM3 and BAM-1.2 at the annual mean time scale, we synthesized results from three land-atmosphere coupling metrics calculated using twelve years  of monthly data (TCI, T_ET, and Zeng's G), and the four climate variables used in their computation (P, ET, surface SM, and surface T), over the four regions in South America shown in Fig. 1, and compared the results against satellite and reanalysis products (Fig. 2, Fig. S4, Table S4). In general, we consider the satellite products as the ''reference'' in our comparison, though we acknowledge that they also have an associated degree of uncertainty, which should not be overlooked.
Satellite SM data were unavailable over three quarters of the Amazon (Fig. S4u, Table S4), though the T_ET and Zeng's G metrics computed using satellite products suggest Amazon land-atmosphere interactions are dominated by an atmospheric control on the land surface (i.e., atmosphere-land coupling, Figs. 2e,i). Positive T_ET values and negative Zeng's G values indicate that ET fluxes in this region are controlled by incoming shortwave radiation. This is consistent with the fact that in very wet areas, such as the Amazon, P is sufficiently high that the land surface is nearly always moist, therefore, the primary limiting factor for ET is the amount of available radiation. Metrics calculated using ERA5-Land, HadGEM3, and BAM-1.2 showed a broader distribution of land-atmosphere interactions in the Amazon, with TCI and T_ET results showing variation in sign across the region (Figs. 2e,i). BAM-1.2 in particular tended to simulate a stronger land surface control on coupling over the Amazon than shown by other datasets (Figs. 2e,i). This finding can be understood by the fact that BAM-1.2 simulates lower surface SM over the Amazon than ERA5-Land or HadGEM3 (Fig. 2u), and therefore BAM-1.2 Amazon ET is more likely to be water limited. Amazon P distributions are similar in all datasets (Fig. 2m), therefore the lower Amazon SM in BAM-1.2 relative to ERA5-Land and HadGEM3 is likely caused by differences in the processes controlling SM.
In the Cerrado, Caatinga, and La Plata, the TCI and T_ET metrics showed that ERA5-Land, HadGEM3, and BAM-1.2 largely captured a land surface control on ET, as indicated by satellite products (i.e., metric results have the same sign, Figs. 2b-d, f-h). Land-atmosphere interactions were strongest over the Caatinga (Figs. 2c,g,k). This is the driest of the regions analyzed, with mean annual P of just 825 mm yr 21 (Fig. S4m), and therefore water availability has an important influence on surface moisture fluxes. BAM-1.2 simulated weaker coupling here than shown in other products, likely because annual P was more than 50% higher than in satellite estimates (Fig. 2o). For the Caatinga and La Plata, the distributions of land-atmosphere interactions in ERA5-Land were closest to the results from satellites. However, HadGEM3 performed slightly better than ERA5-Land over the Cerrado (Figs. 2b,j), which is likely due to ERA5-Land showing higher surface SM than satellites, causing weaker land-atmosphere coupling.
Across all regions analyzed, the TCI, T_ET, and Zeng's G metrics, which all provide a measure of the land surface influence on surface moisture fluxes, showed fairly good agreement in the sign of coupling at the annual scale (Figs. 2a-l). However, there were some differences among metrics in the direction and magnitude of model biases. Zeng's G in particular gave slightly different results from the TCI and T_ET metrics. Given that the three metrics have different input variables (see section 2a), we would not have expected the results from each to be identical. Indeed, climate variables from ERA5-Land and the two climate models over the four South American domains showed variable correspondence with satellite (u)-(x) surface soil moisture; and (y)-(ab) surface temperature in satellite products (white boxes), ERA5-Land (green boxes), and AMIP simulations from HadGEM3 (blue boxes) and BAM-1.2 (red boxes) using data from all months in the year. For the metrics, only pixels with statistically significant coupling ( p , 0.05) were included (Table S4). Boxplots show the quartiles (box) and upper and lower extremes (whiskers) for all pixels over the Amazon, Cerrado, Caatinga, and La Plata regions shown in Fig. 1. Arrows to the right of the top three rows indicate the direction of increasing land-atmosphere coupling strength. Note that the y axes differ between panels to optimize data visualization. The y axes of (e)-(h) have been reversed, to reflect the fact that strong land-atmosphere coupling is indicated by negative values in the T_ET metric. Satellite P came from TRMM, ET from MODIS, SM from ESA-CCI, and T from AIRS (Table 1). The same annual coupling and climate variable data are shown as maps in Fig. S4.
products . Difference between satellite and reanalysis SM did not always mirror those for P, with ERA5-Land tending to represent P relatively well (Figs. 2m-p), but showing higher SM than measured by satellites (Figs. 2u-x). Since Zeng's G uses P as a proxy for SM (supplemental methods), Zeng's G results for ERA5-Land should therefore be interpreted with some caution (Fig. 2j). ERA5-Land and the global models also tended to show higher annual ET than satellites, which may be associated with simulated T that was lower than observed (Figs. 2q-t, y-ab).
Since South America has a strong seasonal hydrological cycle, we evaluated model representation of land-atmosphere interactions using data from 3-monthly periods, to observe how model performance varied seasonally (Fig. 3). The three metrics computed from satellite products showed general agreement in the seasonality of land-atmosphere interactions across regions (Figs. 3a-l), though there were some discrepancies over the Amazon and Caatinga (Figs. 3c,g,k). Land-atmosphere interactions have previously been shown to vary spatially, with the strongest land-atmosphere coupling in wet/dry transition regions (e.g., Koster et al. 2004). Our analysis of satellite products showed seasonal variation in water availability also influenced coupling magnitude (Fig. 3). For example, in the Amazon where atmosphere-land coupling was dominant at the annual scale (Figs. 2e,i), coupling was strongest in DJF when P was highest (Figs. 3e,i,m, Fig. S5a) and radiation at its lowest (Fig. S6a). The direction of interactions was opposite over the Cerrado and La Plata, with variation in the land surface controlling ET. These areas also showed strongest coupling during the wettest months (DJF, Fig. 3,  columns 2 and 4). These results illustrate that due to high variability of P in DJF (Fig. S7), and therefore also radiation and SM, areas with both an atmospheric forcing on the land (e.g., the Amazon), and land surface forcing on the atmosphere Zeng's G metrics calculated using data from 3-monthly periods, plus climatological 3-monthly mean (m)-(p) precipitation, (q)-(t) evapotranspiration, (u)-(x) surface soil moisture, and (y)-(ab) surface temperature from satellite products (black), ERA5-Land reanalysis (green), and AMIP simulations from HadGEM3 (blue) and BAM-1.2 (red) using data averaged over the Amazon, Cerrado, Caatinga, and La Plata regions shown in Fig. 1. For the metrics, only pixels with statistically significant coupling (p , 0.05) were included. Arrows to the right of the top three rows indicate the direction of increasing coupling strength. Note that the y axes differ between panels to optimize data visualization. The y axes of (e)-(h) have been reversed to reflect the fact that strong landatmosphere coupling is indicated by negative values in the T_ET metric.
(Cerrado and La Plata) show strongest coupling at this time of the year. There were differences in seasonal land-atmosphere interactions estimated from ERA5-Land, HadGEM3, and BAM-1.2 over the four study regions (Fig. 3). HadGEM3 and BAM-1.2 captured seasonality in land-atmosphere interactions reasonably well over the Cerrado and La Plata, and to some extent over the Amazon (Figs. 3a-l). Seasonal cycles in P, ET, and SM for the Cerrado and La Plata were mostly well represented in the two models (Figs. 3m-x). For the Amazon, seasonal variation in P was well characterized, with model SM and ET following the same seasonal cycle (Figs. 3m,q,u). However, no Amazon SM satellite data were available for comparison (Fig. 3u), and satellite ET data showed little agreement with models over the Amazon at the seasonal scale (Fig. 3q). ERA5-Land performed least well over the Cerrado, where landatmosphere interactions in DJF showed the opposite sign to satellite estimates (Figs. 3b,f), possibly due to surface SM being overestimated by more than a third in this period (Fig. 3v).

b. Spatial variation in model performance
Figures 4 and 5 show spatial variability in the TCI and T_ET metrics, respectively. Results from the Zeng's G metric showed similar spatial patterns to TCI and T_ET, and are presented in the supplemental material (Fig. S8). Dark red shading in Figs. 4 and 5 indicate ''hot spots'' where surface moisture fluxes are most sensitive to variation in soil water. The regions with the strongest land-atmosphere interactions, as indicated by the satellite TCI results, were northeast Brazil over the seasonally dry Caatinga and Cerrado ecoregions, and northern Argentina, over La Plata basin (Fig. 4, column 1). Since the metric is weighted by variability in surface state, the land influence on the atmosphere was strongest in DJF (see large region of dark red shading in Fig. 4a), in agreement with the results presented in Fig. 3. With the exception of the Cerrado, ERA5-Land, HadGEM3, and BAM-1.2 mostly captured these hotspots (Fig. 4), indicating spatial variation in the ET response to surface SM was well represented over these areas. The maps also reveal where there is spatial variation in land-atmosphere interactions within the defined study regions, thus influencing the regional means presented in Fig. 3, for example, over the Cerrado in MAM (Figs. 4e-h) and the Amazon and La Plata in JJA (Figs. 5i-l).
The Amazon is dominated by an atmospheric control on the land surface throughout the year (blue shading in Fig. 5, column 1). This feature can also be detected in ERA5-Land, BAM-1.2, and HadGEM3, though these datasets tend to overestimate the strength of the signal, while underestimating the spatial extent. This is particularly evident in JJA and SON, when the models show large areas with a land surface control (red shading, Fig. 5, columns 2-4). Output from the Zeng's G metric, which is based on relationships between P and ET, also indicates an atmospheric control on interactions over most of the Amazon throughout the year (Fig. S8). The P-ET relationships have previously been used to distinguish between regions where ET is primarily limited by soil moisture (areas with a positive correlation between anomalies of P and ET) or available energy at the surface (areas with a negative correlation; Baker et al. 2020, manuscript submitted to Environ. Res. Lett.). The fact that Amazon atmosphere-land interactions are strongest during DJF (Fig. 5a) supports this interpretation, as solar radiation is lower during the wetter months of the year (Figs. S5 and S6).
The reanalysis and models performed less well along the North Atlantic coast, particularly in DJF and MAM, simulating strong land-atmosphere coupling that was not suggested by satellite data (Figs. 5b-d, f-h). This indicates that the variables controlling ET, such as SM, P, vegetation functioning and/or cloudiness, may be misrepresented over these regions, leading to misrepresentation of land-atmosphere interactions. Indeed, ERA5-Land, BAM-1.2, and HadGEM3 all showed lower P than the satellite product in this area (Fig. S5), which may explain why moisture availability is shown to be the dominant control on surface moisture fluxes over this region.

c. Seasonal variation in SM-ET relationships over the Caatinga
We used the Betts' approach to investigate the relationship between anomalies of surface SM and ET over the Caatinga by season (Fig. 6), where the models and reanalysis seemed to capture the direction of land-atmosphere interactions relatively well (Figs. 4 and 5). For HadGEM3 and BAM-1.2, results from a single model member are shown in Fig. 6 (r1i1p1 and member 1, respectively), though comparable results were found across all members (Table S5).
Although the period of analysis was only short (2003-14), we identified strong and consistent positive relationships between satellite surface SM and ET throughout the year (Pearson's r 5 0.83-0.87), indicating surface water fluxes are strongly dependent on terrestrial water availability in this region (Fig. 6, column 1). ERA5-Land, HadGEM3, and BAM-1.2 captured positive relationships between surface SM and ET in all seasons (Pearson's r 5 0.80-0.96, Fig. 6, columns 2-4), though there were some differences in the distributions of these variables from those shown in satellite products. For example, in the wet-to-dry transition season (MAM) HadGEM3 showed a wider range in surface SM than shown in satellites (Fig. 6g), implying this model may overestimate variability in the surface state in these months, while BAM-1.2 showed a narrower range in surface SM and ET over the same time period (Fig. 6h). In the dry season (JJA), the models and ERA5-Land all showed narrower distributions of surface SM than shown in satellites. This could explain why coupling in these products is weaker than in satellites in these months (Figs. 3c,g,k). These results show that despite capturing the sign and magnitude of SM-ET relationships, representation of variables impacting land-atmosphere interactions over the Caatinga in ERA5-Land, HadGEM3, and BAM-1.2 could still be improved.

d. Tracing surface-to-atmosphere moisture pathways
When considering land-atmosphere moisture transfer pathways, it can be helpful to distinguish between processes that operate at the land-atmosphere interface and processes that occur in the atmospheric boundary layer. For example, the coupling between the land surface and surface fluxes, and between surface fluxes and the atmosphere, can be quantified separately (Guo et al. 2006;Dirmeyer 2011;Dirmeyer et al. 2014). This can make it easier to understand interactions and feedbacks in a mechanistic way, and to properly evaluate whether models are representing physical processes accurately. We focused on the DJF period, when seasonal analysis showed land-atmosphere interactions were at their strongest over the Amazon, Cerrado, and La Plata (Fig. 3), and therefore, when it is particularly important that models perform well.
The two-legged metric for the land-to-atmosphere moisture transfer pathway (SM-ET-P) in DJF is presented in Fig. 7, and its proxy pathway (T-ET-P) is shown in the supplemental material (Fig. S9). ERA5-Land, HadGEM3, and BAM-1.2 captured the positive association between surface SM and P over most of South America, showing similar spatial patterns of land-atmosphere interactions to those in satellites (Figs. 7i-l). However, on inspection of the surface  and atmospheric (Figs. 7d, legs of the metric, it is possible to detect regions where the reanalysis and models simulated the correct feedback between surface SM and P, but for the wrong reasons. For example, ERA5-Land, HadGEM3, and BAM showed positive correlations between surface SM and ET and between ET and P along the North Atlantic coast (red shading, Figs. 7b-d, f-h), suggestive of a strong surface control on land-atmosphere interactions. However, satellite products showed inverse SM-ET (positive T_ET in Fig. S9a) and ET-P relationships over this region (blue shading,Figs. 7a,d), indicating that the direction of control is not from the surface to the atmosphere, but rather the atmosphere is exerting an influence on the surface. This is consistent with this region having an energy-limited evaporative regime (Guan et al. 2015). Over the Cerrado, ERA5-Land showed the opposite pattern, with negative relationships in the surface and atmospheric legs combining to give a positive SM-P association (Fig. 7, column 2), whereas the results from satellite products show the surface (SM-ET) and atmospheric (ET-P) legs of the pathway are both positive in this region (Fig. 7, column 1). These results highlight how the two-legged metric can provide a useful way of analyzing representation of land-atmosphere feedback pathways, and highlight regions where model processes could be improved. e. Land-atmosphere interactions in CMIP6 models Land-atmosphere interactions in nine CMIP6 models were analyzed at the annual mean time scale and compared with results from the datasets already analyzed in this study (Fig. 8). For each region analyzed, median estimates of landatmosphere coupling from satellites, ERA5-Land, HadGEM3, and BAM fell within or close to the distribution of median values from the CMIP6 models (gray boxes, Figs. 8a-l). In most cases, median coupling values from members of HadGEM3 or BAM-1.2 overlapped with the interquartile range of median values from CMIP6 models, showing that the two models that were assessed in detail in this study are consistent with the median response shown by the wider CMIP6 ensemble. For example, over La Plata, although results from the T_ET metric showed most CMIP6 models tended to simulate stronger land-atmosphere coupling than satellites, ERA5-Land, and HadGEM3, the CMIP6 models were comparable to results from BAM-1.2 (Fig. 8h). For BAM-1.2 and the wider CMIP6 ensemble, the overestimation of land-atmosphere interactions in La Plata could be attributed to low P over the region (Fig. 8p). In another example, CMIP6 TCI and T_ET results over the Caatinga showed relatively strong land-atmosphere coupling, in good correspondence with HadGEM3 (Figs. 8c,g), though coupling in BAM-1.2 was weaker in this region due to P and surface SM being overestimated (Figs. 8o,w). Over the Amazon, six out of nine CMIP6 models were able to capture the atmospheric influence on land-atmosphere interactions indicated by satellites, as evidenced by median T_ET and Zeng's G values that were positive and negative, respectively (Figs. 8e,i). Amazon P, ET, and T from CMIP6 models also showed good agreement with satellite products (Figs. 8m,q,y). However, Amazon surface SM in the CMIP6 ensemble was notably lower than in ERA5-Land, HadGEM3, and BAM-1.2. Outside of the Amazon, nearly all of the CMIP6 models analyzed showed coupling results of the same sign as satellites (e.g., Figs. 8b-d, f-h), though representation of climate variables was variable. A more in-depth of analysis of land-atmosphere interactions over South America in CMIP6 models should be the focus of another study, however, the results presented here suggest that CMIP6 models capture key features in land-atmosphere coupling over the continent.

Discussion
We conducted a detailed assessment of land-atmosphere interactions over South America using satellite products, a novel reanalysis (ERA5-Land) and simulations from two global climate models, HadGEM3 and BAM-1.2. The climatology and variability of these models have previously been evaluated at the global scale (Figueroa et al. 2016;Williams et al. 2018;Kuhlbrodt et al. 2018;Coelho et al. 2021), and landatmosphere interactions over South America were analyzed in similar versions of these models at subseasonal time scales (Chevuturi et al. 2021, manuscript submitted to Climate Resilience Sustainability). Here, we focused on assessing landatmosphere interactions over longer, climate-relevant time scales, with the aim of comparing results from different FIG. 7. Two-legged metric for the DJF surface-to-atmosphere moisture-transfer pathway. Relationships between (a)-(d) a surface variable and a surface flux variable, (e)-(h) a flux variable and an atmospheric variable, and (i)-(l) the full coupling pathway. In this example, surface, flux, and atmospheric variables are surface SM, ET, and P, respectively. The metric was calculated using (first column) satellite products, (second column) ERA5-Land reanalysis, and AMIP simulations from (third column) HadGEM3, and (fourth column) BAM-1.2, for austral summer (DJF). Satellite data are from the ESA-CCI (SM), MODIS (ET), and TRMM (P; Table 1). Stippling indicates where coupling is statistically significant at p , 0.05. For the models, stippling indicates where at least three out of four members showed a statistically significant correlation (p , 0.05). Black lines indicate the four subregions of South America shown in Fig. 1. Gray shading shows where satellite results were unavailable. The proxy pathway, T-ET-P, is shown in Fig. S9.  (Table S3) for the Amazon, Cerrado, Caatinga, and La Plata regions shown in Fig. 1. Spatial medians for satellite products (black stars) and ERA5-Land reanalysis (green triangles) are shown, plus the spatial median values from each of the four members of HadGEM3 (blue circles) and BAM-1.2 (red diamonds). The satellite and ERA5-Land values are the same median values as those presented in Fig. 2. The BAM-1.2 and HadGEM3 values are the same as the median values for individual members presented in Figs. S1 and S2. Note that the y axes differ between panels to optimize data visualization. The y axes of (e)-(h) have been reversed, to reflect the fact that strong land-atmosphere coupling is indicated by negative values in the T_ET metric.
sources, identifying strengths and weaknesses of different products, and comparing our results to the wider ensemble of CMIP6 models. Models must be able to reproduce the direction and strength of interactions between the land and the atmosphere to accurately simulate climate. In this study, results from satellites showed a dipole in the direction of land-atmosphere interactions across South America. Interactions in the Amazon were dominated by an atmospheric control on the land surface (i.e., an energy-limited evaporative regime) throughout the year. This feature was detected in ERA5-Land, BAM-1.2, and HadGEM3, and two-thirds of CMIP6 models analyzed. A recent analysis showed half of CMIP5 models were unable to accurately represent land-atmosphere interactions over the Amazon (Baker et al. 2020, manuscript submitted to Environ. Res. Lett.), so the ability of HadGEM3 and BAM-1.2 to capture this large-scale signal at the annual scale is a strength of these two models. However, both models tended to overestimate the importance of the land surface state (SM) on ET over parts of the Amazon in JJA and SON. The overestimation of land-atmosphere coupling, seemingly driven by underprediction of surface SM, has implications for future climate projections over the Amazon, as drying or wetting trends could be unrealistically exacerbated (Baker et al. 2020, manuscript submitted to Environ. Res. Lett.). Furthermore, future deforestation and climate change in the Amazon are expected to have most impact on dry season climate (Boisier et al. 2015;Guimberteau et al. 2017), which is when differences between models and satellites were largest. Deforestation has already been linked to increases in dry season length over recent decades (Fu et al. 2013;Marengo et al. 2018). Our analysis suggests that in parts of the Amazon where land-atmosphere interactions are misrepresented in the dry season, HadGEM3 and BAM-1.2 may not reliably simulate climate responses to future deforestation. In contrast to most other CMIP5 and CMIP6 models, HadGEM2-ES and the U.K. Earth System Model, successors to HadGEM3, simulate a net cooling in response to Amazon deforestation, partly driven by the ET response to deforestation being too weak (Robertson 2019;Boysen et al. 2020).
Over the Cerrado, the Caatinga, and La Plata basin, landatmosphere interactions were controlled by variation in surface state. These regions have a moisture-limited evaporative regime, since water availability determines the magnitude of ET fluxes (Budyko 1974;Seneviratne et al. 2010). BAM-1.2 and HadGEM3 generally captured the correct sign of landatmosphere interactions over these areas, which have at least five months of the year with less than 100 mm of P, and are examples of ''transition zones'' between wet and dry climate regimes. The CMIP6 models analyzed in this study were also able to capture the direction of land-atmosphere interactions over these areas. Previous work has shown that landatmosphere interactions are often strongest over regions with transitional climate (Koster et al. 2004;Seneviratne et al. 2010;Levine et al. 2016). Accurate simulation of landatmosphere interactions is especially important over areas where coupling is strong. One study, based on four regional climate models, found overestimation of coupling in La Plata contributed to P and T biases, due to positive feedbacks between the surface and atmosphere (Carril et al. 2012). On the other hand, underestimating coupling could cause climate models to underpredict climate variability , thus impacting the reliability of future projections.
Our results confirmed the importance of surface SM for accurate representation of land-atmosphere interactions at annual and seasonal scales. High surface SM over the Cerrado in the wet season resulted in ERA5-Land simulating opposite coupling to that shown by satellite products, despite ERA5-Land being highly constrained by observational data (Hersbach et al. 2020), and being reported to represent SM better than previous reanalyses over most of the Northern Hemisphere (Li et al. 2020). ERA5-Land, HadGEM3, and BAM-1.2 all showed lower surface SM variability over the Caatinga in the dry season than shown by satellite data, contributing to weaker land-atmosphere coupling. SM is key for influencing spatiotemporal patterns of P and ET over South America (Bedoya-Soto et al. 2018), and the land surface ''memory'' of SM anomalies can influence the climate from a few days to several months Ruscica et al. 2014;Levine et al. 2019). Simulation of SM has been shown to be highly dependent on the model land surface scheme (Koster et al. 2009), and our results confirm this, with the reanalysis and models all showing different patterns of SM behavior. SM is influenced by many processes besides P, including factors which affect plant water consumption (e.g., vegetation fraction, leaf area index, and rooting depth), partitioning of heat fluxes (e.g., albedo, roughness length), and properties that determine soil-water storage and runoff (e.g., soil porosity, depth). Targeting these processes for model development could help to improve representation of land-atmosphere interactions, and thus regional climate projections.
Seasonal analysis showed land-atmosphere interactions were strongest during the wet austral summer (DJF) in areas with a land surface forcing on the atmosphere (e.g., Cerrado and La Plata) and an atmospheric forcing on the land (e.g., the Amazon). Previous studies have reported coupling maxima during DJF in southeast South America (Ruscica et al. 2015(Ruscica et al. , 2016. Conversely, Wei and Dirmeyer (2012) observed that climatologically wet regions have a greater land surface influence in drier times of the year, though they did not analyze interactions over the Amazon (due to a lack of SM observations) and were primarily focused on understanding the surface-to-atmosphere forcing, rather than atmospheric forcing on the land surface.
Spatial analysis highlighted where coupling in the models differed from coupling estimated from satellite products. Along the North Atlantic coast of South America, ERA5-Land, HadGEM3, and, to a lesser extent, BAM-1.2 simulated strong land surface forcing on the atmosphere that was not apparent in satellites. Precipitation and SM biases were most likely responsible for the misrepresentation of coupling over this area. A recent analysis suggested that this region could see large absolute reductions in aboveground carbon in response to 28C of global warming (Sullivan et al. 2020), highlighting the importance of reliable climate projections for understanding future changes in the carbon cycle here. Furthermore, since atmospheric water vapor enters the continent from the Atlantic Ocean, having a bias in landatmosphere coupling along the Atlantic coast could impact simulation of downwind precipitation, as the biased atmospheric moisture signal propagates inland via transpiring forests (Staal et al. 2018). In the wet season, when the difference between coupling in satellites and models in this region was greatest, moisture from the Amazon is transported via the South American low-level jet along the Andes and southward to La Plata basin van der Ent et al. 2010;Zemp et al. 2014;Drumond et al. 2014). Model errors in coupling along the Atlantic coast could therefore have detrimental consequences for climate simulation over the whole region.
There are some caveats to the analysis presented in this study. Our evaluation focused on the period of overlap between AMIP simulations and satellite-derived climate products (i.e., 2003-14). This is a relatively short timeframe for analysis, especially at the monthly scale, however, we were still able to detect robust land-atmosphere interactions that were statistically significant. Reanalysis datasets offer climate insights over longer time periods and are often used for model evaluation, however, as indicated in this study, reanalyses may not always provide a reliable representation of reality. We therefore decided to use satellite products as a reference, despite the short period of available data, because satellites can provide measurements over remote areas such as the Amazon, where instrumental data are scarce (Harris et al. 2014). However, there were doubts over the reliability of some of the reference satellite datasets used. Previous studies have shown satellite ET products show discrepancies over the Amazon, and struggle to represent the seasonal cycle in Amazon ET (Maeda et al. 2017;Sörensson and Ruscica 2018;Baker et al. 2021). Furthermore, our results showed that satellite estimates of P and SM, which might be expected to covary, were not always consistent, highlighting another area of uncertainty. Doubts over the reliability of reference datasets present a challenge for climate model evaluation, though quantifying land-atmosphere interactions using several independent metrics based on different data sources, as applied here, may help toward overcoming this issue.
A drawback of the simple correlation-based evaluation metrics described in this study should also be noted. Significant relationships between land surface, flux and atmospheric variables do not provide evidence of causality, and could be caused by covariation with another climate variable. Compensating errors in models can result in relationships that appear to match observations, although the processes driving them in models do not reflect reality (Dirmeyer et al. 2018). The two-legged metric, which was used to trace moisture feedback pathways, was able to detect unreliable land-atmosphere coupling arising from such error compensation, highlighting parts of South America where ERA5-Land, HadGEM3, and BAM-1.2 captured relationships between surface SM and P but for the wrong reasons. Another approach to demonstrate causality is to conduct model experiments, for example, through comparing simulations with and without land-atmosphere coupling (Sörensson and Menéndez 2011;Ruscica et al. 2015;Spennemann et al. 2018;Menéndez et al. 2019;Giles et al. 2020). However, it may be desirable to evaluate land-atmosphere interactions in preexisting model simulations, such as output from global model intercomparison projects (Eyring et al. 2016b), and therefore metrics such as those presented in this study offer a useful alternative approach.

Summary and outlook
We assessed representation of land-atmosphere interactions over South America in satellite products, the ERA5-Land reanalysis and two global climate models, BAM-1.2 and HadGEM3. The two models captured seasonal and spatial variation in coupling strength over South America relatively well, simulating a strong land surface influence on ET over regions of savanna and seasonally dry forest in the south of the continent, and an atmospheric control on ET over the wet Amazon. Some issues with model performance were identified over the Amazon, and along the North Atlantic coast, which could result in climate biases elsewhere. BAM-1.2 and HadGEM3 were broadly consistent with an ensemble of nine CMIP6 models, which captured key features of South American land-atmosphere interactions at the annual scale.
Climate change studies suggest spatial coupling patterns in South America are likely to change over the next century (Dirmeyer et al. 2013). Areas where P is expected to increase, such as La Plata, will tend to see a weakening of land-atmosphere coupling as moisture fluxes become less sensitive to variation in surface SM, while areas that are projected to become drier, such as the Amazon, may show a shift toward a stronger land surface influence on coupling (Llopart et al. 2014;Menéndez et al. 2016;Ruscica et al. 2016;Zaninelli et al. 2019). Crucially, CMIP6 models disagree in their local T and P responses to tropical deforestation, due to different representations of land cover and land cover change (Boysen et al. 2020). Overall, our analysis shows that while global models show promise in their ability to represent landatmosphere interactions over South America, improvements are required for more reliable investigation of how future changes in climate and land use will impact regional hydrological cycling. tl51554219774!4!!&q5modis%20mod16&ok5modis%20mod16& ac5true); TRMM P (https://disc2.gesdisc.eosdis.nasa.gov/data/ TRMM_L3/TRMM_3B43.7/); AIRS T (https://acdisc.gesdisc. eosdis.nasa.gov/data/Aqua_AIRS_Level3/AIRX3STM.006/); CLARA-A1 RDN (https://wui.cmsaf.eu/safira/action/viewDoiDetails? acronym5CLARA_AVHRR_V001); ERA5-Land P, ET, T, SM, and RDN (https://climate.copernicus.eu/climate-reanalysis). The CMIP6 model datasets analyzed in the study are available at https://esgf-index1.ceda.ac.uk/search/cmip6-ceda/. The landatmosphere coupling diagnostic scripts applied in this study can be obtained from https://github.com/jcabaker/land_atm_coupling.