Deforestation reduces rainfall and agricultural revenues in the Brazilian Amazon

It has been suggested that rainfall in the Amazon decreases if forest loss exceeds some threshold, but the specific value of this threshold remains uncertain. Here, we investigate the relationship between historical deforestation and rainfall at different geographical scales across the Southern Brazilian Amazon (SBA). We also assess impacts of deforestation policy scenarios on the region’s agriculture. Forest loss of up to 55–60% within 28 km grid cells enhances rainfall, but further deforestation reduces rainfall precipitously. This threshold is lower at larger scales (45–50% at 56 km and 25–30% at 112 km grid cells), while rainfall decreases linearly within 224 km grid cells. Widespread deforestation results in a hydrological and economic negative-sum game, because lower rainfall and agricultural productivity at larger scales outdo local gains. Under a weak governance scenario, SBA may lose 56% of its forests by 2050. Reducing deforestation prevents agricultural losses in SBA up to US$ 1 billion annually.

Amazon. Beige to blue gradient color indicates the total annual rainfall (mm year -1) calculated by using a second-degree polynomial model. Note that the map clearly shows the signal of the South Atlantic Convergence Zone.
Supplementary Figure 4. 3-step procedure to calculate the anomalies. Diagram with method steps to calculate the annual rainfall anomaly for each grid cell. i, j, t are the subscripts representing space (i,j) and time (t) dimensions, P , , are the annual rainfall values in mm, P , are the estimated values of rainfall due to geographical location and elevation; P * , , is the difference between observed annual rainfall values and the estimated ones due to geographical location and elevation; P ̅ t are the annual averages of rainfall calculated throughout the study region; P , , ′ are the residual annual rainfall anomalies; , , are the subscripts representing latitude, longitude (in degrees), and elevation (in meters). Supplementary Figure 10. Spatial association between forest loss and estimated anomalies for 2019. Spatial variability in anomalies and the Cramer's V 10 and the Spearman Rank Order Correlation coefficient to measure the spatial association between forest loss and estimated anomalies for c and d cells that have reached the critical threshold of forest loss and a and b cells that have not crossed this limit. Brown to blue gradient color indicates precipitation anomalies (mm) for 2019. Red to pink gradient color indicates forest loss (%) within 28 km grid cells. Each point represents the location of a State Capital.
Supplementary Figure 11. Significant effect of Maximum cumulative water deficit on deforestation. One degree grid cells with a significant (α = 0.05) temporal effect of the Maximum Climatological Water Deficit on deforestation (MCWD; % yr -1 mm -1 ). Colors indicate change in deforestation for each mm reduction in MCWD as a percentage of the average annual local forest suppression. Note that cells with statically significant effect concentrate in the highly humid forests of Northwestern Amazon.     Table 7: Two-Sample t-test with unequal sample sizes and unequal variance of means of annual rainfall anomalies (Pi,j,t) for cell groups 1 and 2 for the 1999-2009 and 2010-2019 periods.

1.1.
Modelling the spatial variability of annual rainfall: The TRMM data show a swath of high values of annual rainfall extending from NW to SE ( Figure S1), which may relate to the South Atlantic Convergence Zone (SACZ).
Supplementary Figure 1. Mean observed annual rainfall in southern Brazilian Amazon. Beige to blue gradient colour indicates the total annual rainfall (mm year -1 ) from the Tropical Rainfall Measuring Mission between 1999 and 2019.
We computed a rainfall profile along a transect from NW to SE to evaluate this gradient by testing different non-linear regression models to select the one that provides a good fit to the data (Supplementary Figure 2). Although the third-degree polynomial model yielded a slightly higher fitting, the gains are not so expressive when compared with the results of a second-degree polynomial model. To avoid overfitting common to high degree polynomial models, we opted for the second-degree polynomial model. We also applied these non-linear regression models to all cells of the study region and analysed the distance between known data points and the ones predicted with the regression models (the residual standard error; Supplementary  Figure 3 shows the mean annual rainfall estimates (1999 to 2019) calculated by using a second-degree polynomial model. Note that the map clearly depicts the signal of the SACZ with a good agreement with the observed values (r² = 0.72; p < 0.005). Table S2 shows the model residues and Table S3 shows the summary statistics for the second-degree polynomial model.
Supplementary Figure 3. Mean annual rainfall estimates for the southern Brazilian Amazon. Beige to blue gradient color indicates the total annual rainfall (mm year -1) calculated by using a second-degree polynomial model. Note that the map clearly shows the signal of the South Atlantic Convergence Zone.
Supplementary We tested the assumption that these coefficients are constant over time by randomizing the annual data of the TRMM. We then calculated a correlation coefficient, i.e. Spearman Rank-Order Correlation Coefficients, between predicted values for the entire period of study and predicted values for three different time-periods. We concluded that, although there are small differences between the equations' coefficients, the results from the equation for the entire period of analysis strongly correlates 1 with the ones for the threeseparate time-periods (Supplementary Table 4). As a result, we find that Interannual changes in the spatial patterns are non-monotonic (as would be expected due to e.g.

ENSO)
, and as such, should not bias our estimated effects.
Supplementary To minimize omitted variable bias in our regression analysis, we firstly removed the effects of factors other than deforestation that may affect rainfall across both time and space. To remove the trends associated with geographic location and elevation, as well as part of the interannual variability associated with large-scale climate mechanisms, hence isolating the deforestation effects from those of natural variability, we calculated anomalies of annual rainfall (P'i,j,t) using a 3-step procedure, summarized in Supplementary Figure 4. 3-step procedure to calculate the anomalies. Diagram with method steps to calculate the annual rainfall anomaly for each grid cell. i, j, t are the subscripts representing space (i,j) and time (t) dimensions, P , , are the annual rainfall values in mm, P , are the estimated values of rainfall due to geographical location and elevation; P * , , is the difference between observed annual rainfall values and the estimated ones due to geographical location and elevation; P ̅ t are the annual averages of rainfall calculated throughout the study region; P , , ′ are the residual annual rainfall anomalies; , , are the subscripts representing latitude, longitude (in degrees), and elevation (in meters).
As aforementioned, the spatial pattern of annual rainfall depicts a swath of high values extending from NW to SE that may be associated with the South Atlantic Convergence Zone. This swath of high rainfall values traverses our study region. To remove this largescale geographical pattern, we used increasingly complex regression models on latitude, longitude and elevation calculated using the climatological averages from 1999 to 2019.
We found that a second-degree polynomial model on latitude, longitude and elevation is sufficient to describe this geographical climatological pattern (equation (S1); r² = 0.72; p > 10 -5 ). Using this model, we computed estimated values of rainfall due to this geographical pattern so that: Step 2: Next, we spatially detrended the values of annual rainfall (i.e., to remove the climatological trend related to geographic location and elevation). To do this, we calculated the difference between raw values of observed rainfall in each year and the estimated values due to geographical position obtained from equation (S1), so that.
Since the coefficients of equation (S1) are static over time, the results from equation (S2) represent the observed deviation (anomaly) from the climatological pattern due to the interannual variability.
Large-scale interannual climate phenomena may also affect the annual averages of rainfall values. El Niño Southern Oscillation (ENSO) is the most important interannual variability modulating rainfall in the region, which indeed shows a high correlation with rainfall in our study region (section S1.3; Supplementary Figure 6). Since our time-series is relatively short, there is a bias (albeit weak) toward higher Nino 3.4 indexes, possibly associated with a partial cycle of the Pacific Decadal Oscillation (PDO). As forest losses across the study region increase monotonically as a function of time, Nino 3.4 trend coinciding with forest loss may bias the analysis.
Step 3: To remove the signal of large-scale factors, such as that of the ENSO, we subtracted the mean annual rainfall for the whole study region from the outputs of equation (S2): The residual is assumed to be an 'anomaly' that is not explained by the geographic location, elevation or large-scale time-varying factors.

Verification of the detrending method
To show that the detrending method works, we present a map of difference (residues) between raw values of observed annual rainfall and the estimated values due to geographical location and elevation (Supplementary Figure 5). In general, the difference is between 120 mm and -120 mm, i.e., less than 10% of the annual rainfall.
Supplementary Figure 5. Residues between raw values of observed annual rainfall and the estimated values due to geographical location and elevation. Brown to blue gradient color indicates residues between raw values of observed annual rainfall from the Tropical Rainfall Measuring Mission data and the estimated values after equation (S1).
As mentioned earlier, we attempted to remove the spatial patterns that might be correlated with deforestation and would thus lead to bias in our regression analyses. We added the step 3 in which we subtracted the annual regional average rainfall from equation (2)  After this step, the residual is assumed to be a "rainfall anomaly" that is not explained by the geographic location, elevation and the interannual variability. As expected, after this procedure, the annual means situate around zero (Table S5).
Supplementary   The change in soybean yield in each deforestation scenario is denoted below: The percentage change ∆ % in soybean yield is as follows: 10 . 100 (S6) Change in revenue per hectare is calculated by multiplying the simulated yield and the price of a ton of soybean: Where ∆R ℎ is the change in soybean revenue per hectare in cell i,j (US$ ha -1 ) and Psoy is the price of a ton of soybean (US$ ton -1 ).

1.6.1.2.
Change in revenue per hectare for cattle beef The cattle beef production is derived using the regional stocking rates as of 2012: where C ref is the reference cattle beef production (arroba/ha/year). S ref is the reference stocking rate as of 2012 (heads ha -1 ); W is the national average weight per animal (540 kg head -1 ); r is the average ratio dead weight/live weight (r = 0.41); a is the conversion factor from arroba to kg (a = 15kg arroba -1 ); and ts is the average animal age at slaughter (ts = 2 years).
Change in pasture productivity P in each deforestation scenario is calculated using equations (S9) and (S10): is the pasture productivity in pixel i,j in the standard scenario F10 (ton ha -1 ); and A F is the total area deforested in scenario Fx (km 2 ). If the generic deforestation scenario d is different than one of the standard scenarios Fx (x = 10, 20, …, 40), then the result is interpolated using equation (S10): (S10) The percentage change ∆P % in pasture productivity is denoted as follows: F10 . 100 (S11) Equation (S12) then calculates the change in cattle beef production: The revenue per hectare is calculated by multiplying the cattle beef production per hectare per year and the price of the arroba: where E ℎ is the cattle beef revenue per hectare per year in cell i,j (US$ ha -1 yr -1 ) and P arroba is the price of an arroba of dead weight (US$ arroba -1 ).
Next, we projected productivity change until 2050 from estimates of soybean and pasture productivity losses due to lower rainfall caused by Amazon region-wide deforestation under WEG and SEG scenarios. To do so, we adjusted soybean productivity projections (3.7 ton hectare -1 ) 5 and pasture productivity projections (2.9 arroba hectare -1 ) 6 by the average productivity losses. We computed annual revenues in US$ per hectare using current soybean and cattle arroba prices (US$ 302.58 per ton and US$ 201.50 per arroba, respectively) and projected soybean and pasture productivity for each of the two forest loss scenarios (WEG and SEG) with and without decreases in yields. Under WEG and SEG scenarios, total annual revenues are calculated by weighting soy and pasture revenues using the simulated future croplands and pasturelands. Next, we calculate the Net Present Value (NPV; Eq. S14) of future revenues for the Southern Amazon using as discount rate the Special System of Clearance and Custody rate (Selic interest rate) of 3.75% (the Selic interest rate is determined by Brazil's Central Bank).
Where t is the time of the cash flow, i is the Selic interest rate and Rt is the net cash flow i.e. cash inflow -cash outflow, at time t. We computed the difference between the total NPV of revenues under the WEG versus the SEG scenario is the opportunity cost of the SEG scenario.
where i is the Selic discount rate, NPV is the net present value calculate trough equation (S14) and N is the simulation period (i.e., the planning horizon, in years). Similarly, the difference between total EAA of revenues under the WEG versus the SEG scenario is the annually opportunity cost of the SEG scenario. The fact that the reduction in rainfall could be smaller than the typical year-to-year variability in some of the grid cells does not mean that there is no impact on agricultural production. Since we have isolated the signal of deforestation, the effect of deforestation on rainfall is additional to the interannual variability, further increasing risks of crop failure especially in dry or drought years. Furthermore, crop yields and the adoption of double cropping systems vary significantly from year to year in the region, partially because of the interannual climate variability 8,9 . Therefore, it is reasonable to assume that even changes that are smaller than the interannual variability may be relevant to agriculture in the region.

Evidence for forest loss causation of rainfall reduction
We also evaluated the spatial variability in anomalies calculating the Cramer's V 10  The relation between forest loss and rainfall presented here is correlational. Recently study point out to that forest loss causes 4% of droughts, while droughts account for 0.13% of forest loss per millimetre of rainfall reduced in the Amazon biome 12 . To investigate the possible influence of a drier climate on deforestation in SAB, we superimposed data of Staal et al. 12 on the map of our study region. We found that only 12 cells of a total ≈ 86 non-null cells (14% of our study region) show a statistically significant effect that more severe forest loss tends to take place in drier (or drying) parts (Supplementary Figure 11; Figure was created using data from Staal et al., 2020 12 ). We also found for the few cells that have statistical significance (α = 0.05) that this effect is ≈ 0.09 % yr -1 mm -1 in our study region, i.e. near zero. In other words, only in these few areas forest loss may significantly increase (though very little according to the authors' own results) as climate becomes drier. Furthermore, the same authors point out that forest loss has intensified dry seasons especially in the southwestern Amazon (our study region), stressing that the causality of forest loss on rainfall in our study region is more important than the other way around. Additionally, empirical evidence that the length of the rainy season decreases ≈ 0.9 ± 0.34 day per each 10% of forest loss has been published recently 13 .
Supplementary Figure 11. Significant effect of Maximum cumulative water deficit on deforestation. One degree grid cells with a significant (α = 0.05) temporal effect of the Maximum Climatological Water Deficit on deforestation (MCWD; % yr -1 mm -1 ). Colors indicate change in deforestation for each mm reduction in MCWD as a percentage of the average annual local forest suppression. Note that cells with statically significant effect concentrate in the highly humid forests of Northwestern Amazon.

Spatially-explicit signal of the reduction of rainfall in the region
We tested whether the mean annual rainfall in Southern Amazon is different between two periods: 1999-2009 and 2010-2019, using Two-Sample t-test. Results are significant at 5% level of significance (Table S8). As a result, Supplementary Figure 13 depicts a clear signal of rainfall reduction, mostly notable where deforestation concentrates in Rondônia and Eastern Amazon.
Supplementary Table 8