Final Response to Reviewers: Dominant effect of increasing forest biomass on evapotranspiration: Interpretations of movement in Budyko Space

During the last six decades, forest biomass has increased in Sweden mainly due to forest management, with a possible increasing effect on evapotranspiration. However, increasing global CO2 concentrations may also trigger physiological water saving responses in plants that induce an opposite effect. Additionally, changes in other forest attributes may also affect evapotranspiration. In this study, we aimed to detect the dominating effect(s) of forest change on evapotranspiration by studying changes in the ratio of actual evapotranspiration to precipitation, known as the evaporative 20 ratio, during the period 1961–2012. We first used the Budyko framework of water and energy availability at the basin scale to study the hydroclimatic movements in Budyko space of 65 temperate and boreal basins during this period. We found that movements in Budyko space could not be explained by climatic changes in precipitation and potential evapotranspiration in 60% of these basins, suggesting the existence of other dominant drivers of hydroclimatic change in these basins. In both the temperate and boreal basin groups studied, a negative climatic effect on the evaporative ratio was counteracted by a positive 25 residual effect. The positive residual effect occurred along with increasing standing forest biomass in the temperate and boreal basin groups, increasing forest cover in the temperate basin group and no apparent changes in forest species composition in any group. From these three forest attributes, standing forest biomass was the one that could explain most of the variance of the residual effect in both basin groups. These results further suggest that the water saving response to increasing CO2 in these forests is either negligible or overridden by the opposite effect of the increasing forest biomass. Thus, we conclude that 30 increasing standing forest biomass is the dominant driver of long-term and large-scale evapotranspiration changes in Swedish forests.

In general, the hypothesis and analysis were well thought out and executed. R1: Thank you for this contextualization of our work and for appreciating its contributions. We have also included the new important references in the text. Q2: Broadly, I think the clarity of the paper could be improved. And, more specifically, I have some difficulty understanding the analysis linking evapotranspiration changes to forest expansion rather than CO2 fertilization. "Possible physiological water saving response to rising CO2" The authors state in their title that the ET increase from forest expansion overrides ET decrease from rising CO2. The authors then present evidence that (1) the aridity index decreased over time due to an increase in precipitation; (2) this decrease in aridity index lead to a decrease in E/P, as expected from the Budyko curve; and (3) there was an overall increase in E/P. This increase in E/P was then attributed to changes in forest standing biomass. There was no evidence or estimate of the CO2 fertilization effect and, therefore, I find it difficult to conclude that this effect was present and indeed over-compensated. To address this issue, I would suggest the authors reduce the focus on CO2 fertilization in the title and abstract.
R2: Thanks for this suggestion. It is correct that we did not directly detect or estimate the possible CO2induced water saving response, but rather conclude from the positive direction of change in the nonclimatic component of the evaporative index that it be must small or absent. We removed the water saving response from the title. In addition, we have followed the Reviewer's advice and reformulated the abstract to mention that the water saving response is one of many forest attributes that may influence the evaporative ratio and evapotranspiration: "During the last six decades, forest biomass has increased in Sweden mainly due to forest management, with a possible increasing effect on evapotranspiration. However, increasing global CO2 concentrations may also trigger physiological water saving responses in plants that induce an opposite effect. Additionally, changes in other forest attributes may also affect evapotranspiration". We also now explain that an interpretation of our results is that "These results further suggest that the water saving response to increasing CO2 in these forests is either negligible or overridden by the opposite effect of the increasing forest biomass." Q3: My intuition is that the CO2 fertilization effect was relatively weak in this ecosystem. The atmospheric CO2 concentration increased approximately 85 ppm (315 too 400 ppm) over the study period, 1961-2012. In a meta-analysis of the FACE experiments (Ainsworth and Rogers 2007), trees showed one of the lowest responses of stomatal conductance to elevated CO2 (20% decrease). In these experiments, CO2 was increased from 366 ppm to 567 ppm, an increase 2.35 times that experienced from 1961 to 2012. As a first guess, one might expect less than 10% decrease in stomatal conductance whereas the authors show that forest standing biomass increased 25% and 55% in boreal and temperate watersheds ( Figure 8).
Further, in several places, the authors make the same argument based on results of species-specific responses. See page 2, lines 9-20 and page 7, lines 14-23. To paraphrase, the watersheds studied are dominated by coniferous species and CO2 water saving response has not been observed in these species. In conclusion, there is little evidence to expect a CO2 water saving response in the studied watersheds. R3: We now state more clearly in the Introduction that it is indeed possible that there exists a water saving response by citing the work of Kellomäki and Wang (1996) with the coniferous species Scots pine and by Rey and Jarvis, (1998) with the deciduous species Silver birch. Both experiments have found water saving responses (P2 L29-38).
Q4: Therefore, I do not think it is appropriate to set up the paper with this hypothesis that is later not supported with the data. R4: Thanks for this suggestion. We modified the main objective/hypothesis of the study by removing any specific mention to the stomatal water saving response: "Our main objective was to determine in Swedish forests the dominant forest effect(s) driving changes in the ratio of actual evapotranspiration to precipitation, known as the evaporative ratio, from a basin scale approach" (P2, L40-41). We also improved the justification in the Introduction section of any discussion in the manuscript about the water saving response (P2, L5-12).
Q5: On the other hand, I do think it is appropriate to address this weak CO2 effect as one reason why an increase in E/P was observed, as the authors have done on page 7. R5: Thanks. We also now state that the water saving response may be even negligible as implied from our results, explaining in the discussion how the results agree with other previous tree experiments in the region showing a negligible stomatal water saving response to increased CO2 (New Section 4.2). See also responses to comments above.
In my opinion, the authors have made a convincing argument for residual changes in basin-level ET that goes the extent dictated by climate. They have postulated that, because the observed ET has increased despite a decrease in aridity index (when Budyko's curve, under stationary conditions, would suggest otherwise based solely on climatic effects), there must exist some non-climate related mechanisms that offset this increase. R1: We thank the Reviewer for the good summary of our work and for appreciating our "convincing arguments" and analyses.

Q2:
To make this point, however, I think that Figure 4 is redundant with Figure 6. Figure 4's use of "wind roses" does not add additional support for the authors' main point. While they claim that these wind roses are "a simple way to synthetize general tendencies of movement," the general direction of these movements are well summarized by the histograms presented in Figure 6, so to me these are two different graphical representation for very similar sets of information.
R2: Thank for the suggestion. The reviewer was right in the fact that the roses (Fig. 4) and the boxplots ( Fig. 6) showed similar data. Because of this, we have expanded the description on the way that these roses are analyzed (Figure caption of figure 4) and added more of the important results that can be extracted from such figures (P6 L18-24). These roses show the simultaneously all changes in the aridity index and evaporative ratio for all basins, which becomes an advantage when processing a hydroclimatic change analysis for a large number of basins. The roses are also necessary in order to reach an important conclusion mentioned in the Abstract: "We found that movements in Budyko space could not be explained by climatic changes in precipitation and potential evapotranspiration in 60% of these basins, suggesting other dominant drivers of hydroclimatic change in these basins". Furthermore, we have now modified Figure 6 to include an uncertainty assessment that cannot be included in the roses of Figure 4; this also accentuates the difference between both figures.
As such, both figures are now vital for the red line of the manuscript. We also expanded on the usefulness of each figure by describing them in more in detail and by extending the analysis of the results of the roses (P6 L9-28). The main result of the roses, "…that changes in precipitation and potential evapotranspiration could not explain observed changes in the evaporative ratio in 60% of these basins in both basin groups" is now mentioned fully in the abstract, results, discussion and conclusion to better justify the use of the roses.
Q3: In addition, "spectra of movements in Budyko space," used repeatedly in Section 3.1, need to be rephrased. Since "spectra" has a very specific meaning in time series analysis, I would suggest the authors avoid this term in reference to movement in the Budyko coordinates.
R3: We removed any reference to "spectra" and referred to the wind-rose type diagrams as "roses" instead.
Q4: I think also that the weakness of the manuscript as it stands lies in formulating the argument for the second point, e.g., in attributing the observed increase in ET to a specific, hypothesized mechanism.
R4: Thanks for this suggestion. We have eliminated from the manuscript any reference or mention of attribution of a "particular hypothesized mechanism". We now state that standing forest biomass is the forest attribute that can best explain the variance of the residual of the evaporative ratio (E/P)r (P7 L19, L25, L30 and P9 L4) in comparison to forest cover and forest composition and even the fraction of precipitation falling as snow. Please also now refer to the new methodology used to reach this conclusion (See Section 2.5 of the Methods). We have now removed any sentence mentioning that it is the "only" specific mechanism.
Q5: In Figure 8, boreal and temperature forests showed opposite changes in this deciduous proportion, though how this might contribute to the overall increase in ET in both forest types is not discussed.
R5: Regarding the changes in deciduous proportion in the forests, we now state more explicitly that the changes in forest composition in the temperate and boreal basin groups are smaller than the sampling standard error, so this study cannot actually confirm a statistically significant change in forest compositing in any of the two groups between the periods 1961-1986and 1987 Figure 7). Additionally, forest composition cannot statistically explain the variance of the residual of the evaporative ratio in any of the two basin groups along the period 1961-2012 ( Figure 9 and Table S1 of the Supplementary materials.

Q6:
In addition, the relationship between forest attributes and Ψr is described in Section 3.2 using very vague terms like "in agreement with" and "followed that of." I would suggest applying more statistical analysis (and plot out the correlation between Ψr and each of the forest attributes) in this section to more quantitatively describe these relationships.
R4: Thanks for this valuable suggestion. We removed the former figure from which we had performed the qualitative analysis. In its place, we now include a statistical analysis focusing on the regressions between the residual of the evaporative ratio (E/P)r and each of the forest attributes (A, V, Qv) in order to determine which of the forest attributes can better explain the variance of the residual of the evaporative ratio (i.e., which are the dominant effects on the residual). We now presents these results in Figure 8 and Table S1 of the Supplementary materials and make a thorough quantitative description of these results in P7 L10-18.
Q5: I also remain unconvinced of the authors' use of the cumulative [Ψr] in comparison to the forest attributes (Figure 9), and the application and choice of a 5-year moving window for [Ψr]. Both of these usages require further justification. R5: We have removed the previous confusing parameter "cumulative [Ψr]" and replaced it with the much simpler one, which is the residual of the evaporative ratio (E/P)r=E/P-(E/P)c (P5 L36). In this way, we can use the annual values of (E/P)r to perform the regression analysis just described in response 4.
The idea of the original five-year mean moving window was to reduce the effect of water storage change for the calculations of evapotranspiration in each basin. What we do now is to divide the period 1961-2012 into three-year and five-year periods and also perform the regression analysis on the means of the hydroclimatic and forest data of these periods. Regressions are also done using annual data, and results for one-, three-and five-year data are shown in Supplementary Table S1. We explain this in detail in Section 2.5-Regression analysis of the residual effect and forest attributes.

Q6:
If the authors can address these concerns, this paper will make a good contribution to the study of water partitioning at high latitudes.

Response to Reviewer Nr. 3
We thank Reviewer Nr. 3 for explaining in detail the concerns regarding the robustness of our results as for suggesting options to increase such robustness. We have now: 1) Incorporated the uncertainty of P and E0 into the analysis, 2) Performed a boot strapping procedure to the forest attribute data to account for inherent uncertainty, and 3) Developed a statistical and regression analysis to determine in a more quantitative manner the parameters that can best explain the variance of the residual of the evaporative ratio ad its change. We describe these in the following answer to the Reviewer: Q1: This paper uses the Budyko framework to study the effect of changes in evaporative ratios at a number of boreal and temperate catchments in Sweden. The study looks at changes in the location of each catchment in Budyko space during two consecutive 25-year periods in the early 21st century and second half of he 20th century, and separates the changes into climatic and non-climatic effects. The significant non-climatic effect is then attributed to forest expansion. However, I have a few methodological concerns (detailed below) that leave me concerned about the robustness of the results. I also find the analysis of the results to be fairly limitedthe temperate vs. boreal differences are barely discussed for example. R1: We thank the reviewer for the recommendations. In order to highlight the differences between the boreal and temperate basin groups, we have included information on species for each biome group (temperate and boreal) in the new panel Fig 1c and described the differences accordingly (P5 L14-20). We have now quantified the increase in standing forest biomass of both the deciduous and coniferous fraction of the volume (Fig. 7b). Furthermore, we have also discussed in more detail the differences of the results between both basin groups in the analysis of each figure. For Figure 3, we mention the differences between temperate and boreal forests in P6 L13-15. For Figure 4, P6 L20-24 and the figue caption. For Figure 5, in P6 L29-34. In Figure 6, there is no important difference between boreal and temperate forests. In Figure 7, P7 L4-9 and in Figure 8 in P7 L10-18. As the Reviewer can now see, we now deal with the difference between both basin groups through the entire manuscript.
Furthermore, we now state more explicitly that the changes in forest composition in the temperate and boreal basin groups are smaller than the sampling standard error, so this study cannot actually confirm a statistically significant change in forest compositing in any of the two groups between the periods 1961-1986and 1987. Additionally, forest composition cannot statistically explain the variance of the residual of the evaporative ratio in any of the two basin groups ( Figure 8 and Table S1).

Q2:
The only real result presented is a qualitative statement of relative dominance that confirms previous studies nor is the amount of variability in climatic and vegetation drivers within each biome (despite data on this clearly being used before aggregation in this study).
R2: Thanks for this critical comment. We have modified the manuscript to include a statistical analysis focusing on the regressions between the residual of the evaporative ratio and each of the forest attributes (A, V, QV) in order to determine which of the forest attributes can better explain the variance of the residual of the evaporative ratio (i.e., which are the dominant effects on the residual). We now presents these results in Figure 8 and Table S1 of the Supplementary materials and make a thorough quantitative description of these results in P7 L17-26. Concerning "the amount of variability in climatic and vegetation drivers within each biome", we respond to the reviewer in the next response R3.
Q3: My methodological concerns are as follows: 1) It is argued that forest inventory data cannot be used because they represent too large of an area (e.g. a county that may be larger than the watershed of study within it). In response, the authors aggregate the data even further, to cover an even larger area! How do we know that forest changes and climatic changes are consistent across all of the temperate and all of the boreal areas? The authors should assess the spatial variability of both forest inventory… R3: Thank you. We have extended the explanation on the NFL forest dataset in the manuscript in order to answer the reviewer's questions. We now specify in a more explicit manner that a spatial aggregation into large basin groups is necessary in order to reduce the standard sampling error of the forest attribute data. Please refer to text from in P5 L5-12. Furthermore, we now explain that a more refined spatial data for forest attributes is not available in the Swedish Forest Inventory, and that all the sample plots falling within the area of the basins here studied have been used in this study in order to account for the spatial variability of forests within the basins. The inherent uncertainty of the forest data is shown in terms of the sampling standard error of biomass, area and forest composition in Figure 7. This error is then the uncertainty related to the forest data. It is worth noting that as suggested by the reviewer, we now perform a bootstrapping procedure in order to account for the forest data uncertainty in the statistical regression assessment used to determine in a quantitative manner the parameters that can best explain the variance of the residual of the evaporative ratio and its change. Please refer to Section 2.5-Regression analysis of the residual effect and forest attributes.
Q4: …and rainfall data in each biome to ensure this is a reasonable approach.
R4: Now, in order to study the spatial uncertainty of precipitation (P), we have now added to the analysis the daily P from the Climatic Research Unit-gridded P product CRU TS3.23 (Harris et al., 2014) and an additional P product of the Luftwebb (http://luftwebb.smhi.se/) portal of SMHI. The latter is a biascorrected gridded dataset of precipitation for Sweden during the period 1961-2014 with a 4 x 4 km horizontal resolution and based in data collected from over 87 precipitation stations around the country. A further assessment of P uncertainty is not possible since we are now using all available data for precipitation in the analysis. Of course, we agree with the author in the fact that precipitation may still have a fair amount of spatial variability between stations, but the effects of such variability cannot be assessed due to spatial data limitations. We have also used three data sources for potential evapotranspiration E0, so that in the current analysis each basin has now nine possible combinations of annual P and E0 data area now used in an uncertainty analysis. Please refer to Figure 6 and its corresponding methods (Section 2.4-Calculation of the residual effect on the evaporative ratio).

Q5
: 2) Similarly, LAI is calculated by using a constant leaf mass per area and biomass data from biomeaggregated NFI data (I think the exact treatment of the NFI data is not clearly explained in Sec. 2.4). The authors then argue that LAI and areal forest cover is constant even as biomass increases by 23%. This would imply a huge trend in stem and branch biomass without any changes in other forest properties, which seems somewhat unlikely. R5: We must clarify that our results did not show that the leaf area index of deciduous and temperate forest cover was constant in time. What was constant in time was the ratio of deciduous LAI to total LAI (LAIQ), which was used here as a proxy for forest composition.
Q6: Have the authors checked whether there are changes in forest composition over that time?
R6: Yes we have. The LAIQ was issued precisely for that, as a proxy of forest composition. We had found that the change in LAIQ was not significant since it was not larger than the sampling standard error propagated from the measurement errors involved to calculate it.
In this new version of the manuscript, we have removed the LAI analysis from the manuscript as it only permitted a restricted analysis for the period 1978-2012, which is the time of availability of leaf mass data per species. Instead, we now use the ratio of standing forest volume of deciduous trees to total standing forest biomass (Qv), as a proxy for forest composition since it can be calculated for the entire period 1961-2012. The corresponding propagated uncertainty of Qv is shown in Figure 6. As with LAIQ, the figure shows that the forest composition expressed as Qv did not changed significantly from 1961-1986 to 1987-2012 since the changes in both basin groups were much smaller than the propagated sampling standard error of Qv. Additionally, the new regression analysis showed that Qv could not statistically (p<0.05) explain the variance of the residual of the evaporative ratio in any of the two basin groups ( Figure 8 and Table S1).
Q7: What is the uncertainty induced in the LAI calculation based on assuming constant LMA for two species, and no other species contributions, however small? Furthermore, the statement that LAI is constant on page 6 line 32, directly contradicts the statement that LAI is changing on page 7 line 9.
R7: Please see our responses R5 and R6. LMA is supposed to be constant for each species. It is the leaf mass per unit area as it is a typical characteristic of the species.

Q8
: 3) Even for a catchment with unchanging vegetation conditions, there can be quite a lot of scatter on where a specific catchment's point falls relative to a theoretical Budyko curve due to interannual variability and imperfects in the Budyko framework. While a 50-year average may reduce noise to some degree, the entire climatic vs. non-climatic calculation is potentially highly sensitive to the exact value of n used. Some bootstrapping and uncertainty propagation for n would be helpful for demonstrating that the results are robust.
R8: Thanks for mentioning this point. Based on this suggestion of the Reviewer and that of Q4, we have now performed an uncertainty analysis of P and E0. This analysis gives us now nine new possible combinations of P and E0 that in turn by propagation imply nine possible values of n (Eq. 2) for each basin (and more for each basin group composed of several basins). The further corresponding propagated uncertainty ranges of the calculated arithmetic and area-weighted means of ΔE/P, Δ(E/P)c and residual Δ(E/P)r for each basin group are now shown in Figure 6 as blue and pink vertical ranges of uncertainty, respectively. The results show again the robustness of the counteracting signal between Δ(E/P)c and residual Δ(E/P)r in both basin groups despite the uncertainty of P, E0 and n. WE again thank the reviewer for this important suggestion. Q9: 4) As both the introduction and discussion mention, changes in the fraction of precipitation falling as snow could have a significant effect on the evaporative ratio (Berghuijs et al., 2014) that is not captured in the present analysis. A study of similar effects in China studying the effects of such a change is dismissed for making unrealistic assumptions, but that does not mean that the change itself could not be a factor here. The authors should at a minimum check if there are trends in the fraction of precipitation falling as snowfall. This is particularly troubling since Figure 5 shows a significant change in the seasonal cycle of rainfall in temperate areas.
R9: We agree on the importance of accounting for changes in the fraction of precipitation falling as snow (fs). This is why we now include fs in the new regression analysis used to identify the attributes/parameters that can better explain the variance of Δ(E/P)r (P6 L5-9). The regression analysis found that standing forest biomass V was again the parameter that could best explain (E/P)r across basin groups. Although not as evident as the effect of V, the almost significant p-value (p=0.06) and relatively high R 2 (R 2 =0.27) in the boreal basin group when using three-year mean values for the regression suggests a possible important influence of fs on (E/P)r. We mention this in P7 L20-22 of the Results and devote a section of the Discussion section to hypothesize on such influence (Section 4.3. Intra-annual hydroclimatic changes and their effect on the evaporative ratio). Regarding the trends in fs, we had indeed found that fs had decreased from the first period to the second period as we again mention in this section.

Q10:
The presentation of this paper could be significantly improved in several areas: 1) The specific Ep dataset used in Figures 3-6 is never stated.
R10: The data sets of P, T, Tmin and Tmax used to calculate E0 and the four models/products used for the calculation were and are once again mentioned and explained in P3 L22-25.

Q11
: 2) I find Figure 4 quite hard to follow. Why are the colors not the same across the 4 sub-plots? This would be easier to read. If the colors represent the radius of each paddle, why are different paddles reaching the same radius colored different (e.g. 4a). In addition, how is the r chosen for each paddle, given that it presumably represents multiple catchments? R11: Thanks to the Reviewer for making us realize that we had not explained well how to read these roses. The roses show movement direction and displacement in Budyko space in the same way that a wind rose shows information on wind direction and speed. We have now expanded on the interpretation of the Budyko roses and provided an example on how to read them in the figure caption. The radius of the roses is in fact not relevant for their interpretation. The colors represent the variable r calculated by Eq. 4 and the radius.

Q12
: 3) Figure 6 suggests differences in the climatic vs non-climatic effects magnitudes between boreal and template. Possible reasons for these differences should be mentioned in the Discussion section, since this is one of the main ways in which your analysis allows detailed study. For example, are there differences in composition?
R12: Although there are differences between the results for the two basin groups, the differences of ΔE/P, Δ(E/P)c and residual Δ(E/P)r between them were rather small. The counteracting effect of a positive Δ(E/P)r and a negative Δ(E/P)c and hence the increasing (E/P)c was present in both basin groups. We now show based on the uncertainty analysis suggested by this reviewer, that this result is independent of the assumptions of P and E0.
Concerning the differences in forest composition, we agree with the Reviewer in the fact that we had to expand on the role of different forest composition in both biome groups. Please refer to R1 to this reviewer. However, changes in forest composition were smaller than the sampling standard error in any of the two basin groups (Figure 7). Q13: 4) Can the authors comment on whether possible changes in air quality may play a role? R13: Yes, since air pollution such as ground-level ozone and acid rain affects the canopy of the forest and the LAI, evapotranspiration and the evaporative ratio may be lower than under unaffected conditions (i.e., a negative change). Such is the case of the study by Renner et al. (2013) mentioned in Line 23 of Page 2. However, since the residual of change in the evaporative ratio is in our study positive, we can assume that any possible effect of air quality (which should most likely have been negative, decreasing LAI and annual E) in this case is not as important as that of increasing biomass.
Q14: Other minor comments: R14: We removed the word "former" Q15: Page 3, line 15: Would be helpful to explain 1986 is the midpoint of your data period R15: We mentioned now that the two 26-year sub-periods are equal in length.
Q16: Page 7, line 16: This is not really a conflict with global studies. Even if global average trends are a certain way, showing that a specific location does not follow them is not a contradiction but indeed just a sign of spatial variability -CO2 effects can still dominate elsewhere and therefore for the global average cycle. However, see also Swann et al, PNAS 2016 for additional discussion on this topic. R16: Thanks. We removed the statement saying that our results were in contradiction with global studies and included the reference of Swan et al. 2016.
Q17: Page 7, line 33: That "most of [drainage] was implemented before the present study period" conflicts with you statement that there is a peak in forest drainage implementation in the late 1970's and 80's (line 31) R17: We rewrote this sentence. Most of the drainage occurred before our study period, but there was also a later smaller peak within our period, in the late 1970's and 80's.

Response to Reviewer Nr.4
We thank Reviewer Nr. 4 for highlighting that our results are important and should be communicated and for proposing valuable suggestions to improve the manuscript. We have addressed below the Reviewers remarks, questions and suggestions as follows: Anonymous Referee #4 Q1: The manuscript by Jaramillo et al., 2017 analyses long term changes in ET/P and PET/P of Swedish catchments. The topic is of general interest and well suited for the journal. Many catchments show increasing ET/P even though the aridity index is decreasing due to slightly higher precipitation rates. The data is compared with forest inventory data which shows a significant increase in forest biomass. The data suggests that the overall increase in biomass is the dominant driver in increasing ET and thus ET/P. The authors argue that this "overrides physiological water saving responses". I do agree, however, the improved water use efficiency due to higher CO2 levels might still be an important effect and could decrease ET, but clearly only for the same amount of biomass. For biomass aggregated results are presented, but there is no estimate of the physiological water saving response. R1: It is correct that we did not directly detect or estimate the possible CO2-induced water saving response, but rather conclude from the positive direction of change in the non-climatic component of the evaporative index that it be must small or absent. We have removed the water saving response from the title and reformulated the abstract to mention that the water saving response is one of many forest attributes that may influence the evaporative ratio and evapotranspiration. We also now explain throughout the manuscript that our results suggest, but not specifically show, that "the water saving response to increasing CO2 in these forests is either negligible or overridden by the opposite effect of the increasing forest biomass". We have also now reformulated the hypothesis of the study by excluding the water saving response as follows: "Our main objective was to determine in Swedish forests the dominant driver(s) of changes in the ratio of actual evapotranspiration to precipitation, known as the evaporative ratio, from a basin scale approach" (P2, L40-41).
Q2: Furthermore, the time series data does not provide statistical correlations between biomass and ET/P. Thus the study can not provide quantitative links between biomass or physiological water saving response and ET/P. However, both topics are suggested by the title and hypotheses. Therefore, I recommend to adapt the red line of the manuscript or improve the analysis. Nevertheless, the observation that increases in forest biomass are potentially linked with increasing ET/P is important and should be communicated. R2: We thank the reviewer for the suggestions. As we said in R1, we have reduced the emphasis on the water saving response. Based on this Reviewer's suggestion, we have also now improved the analysis on the attribution of the effects of the forest attributes on the residual of the evaporative ratio. We now include a quantitative approach that uses a statistical regression analysis between the residual of the evaporative ratio and each of the forest attributes (A, V, QV) in order to determine which of the forest attributes can better explain the variance of the residual of the evaporative ratio. We now present these results in Figure 8 and Table S1 of the Supplementary materials and make a thorough quantitative description of these results in P7 L10-18. In general, we found that as mentioned in the abstract "standing forest biomass was the attribute that could best explain the variance of the residual of the evaporative ratio in both biomes.
Thanks once more for all the recommendations.
Abstract. During the last six decades, forest biomass has increased in Sweden mainly due to forest management, with a possible increasing effect on evapotranspiration. However, increasing global CO2 concentrations may also trigger physiological water saving responses in plants that induce an opposite effect. Additionally, changes in other forest attributes may also affect evapotranspiration. In this study, we aimed to detect the dominating effect(s) of forest change on evapotranspiration by studying changes in the ratio of actual evapotranspiration to precipitation, known as the evaporative 20 ratio, during the period 1961-2012. We first used the Budyko framework of water and energy availability at the basin scale to study the hydroclimatic movements in Budyko space of 65 temperate and boreal basins during this period. We found that movements in Budyko space could not be explained by climatic changes in precipitation and potential evapotranspiration in 60% of these basins, suggesting the existence of other dominant drivers of hydroclimatic change in these basins. In both the temperate and boreal basin groups studied, a negative climatic effect on the evaporative ratio was counteracted by a positive 25 residual effect. The positive residual effect occurred along with increasing standing forest biomass in the temperate and boreal basin groups, increasing forest cover in the temperate basin group and no apparent changes in forest species composition in any group. From these three forest attributes, standing forest biomass was the one that could explain most of the variance of the residual effect in both basin groups. These results further suggest that the water saving response to increasing CO2 in these forests is either negligible or overridden by the opposite effect of the increasing forest biomass. Thus, we conclude that 30 increasing standing forest biomass is the dominant driver of long-term and large-scale evapotranspiration changes in Swedish forests.

Introduction
Boreal and temperate forests provide important ecosystem services at local, regional and global scales. These services include water regulation, soil stabilisation, biodiversity preservation, and the provisioning of timber, fibre, fuel, food and cultural 35 values for humans (Chopra, 2005). Changes to the attributes of these forests, such as forest coverage, have had important implications for the functioning of the Earth's climate and water and carbon systems (Abbott et al., 2016;Sterling et al., 2013). The return flow of water vapor from the earth's surface to the atmosphere, known as evapotranspiration, is a key hydroclimatic variable linking these systems. Major regulators of forest evapotranspiration include forest biomass (Feng et al., 2016), leaf stomatal conductance (Lin et al., 2015), canopy leaf area index (LAI) (Mu et al., 2013), tree hydraulic traits (Gao et al., 2014) 40 and stand surface roughness (Donohue et al., 2007). In the context of temperate and boreal Sweden, changes in forest attributes are likely to play an important role because of the higher proportion of available energy partitioned into latent heat than other ecosystems in these regions, such as grassland, wetlands and tundra (Baldocchi et al., 2000;Kasurinen et al., 2014;van der Velde et al., 2013).
However, the dominant effects of forest change on evapotranspiration are still a matter of debate. Zeng et al. (2016) states that 5 the increase in total canopy leaf area linked to the recent "greening" of the Earth is responsible for more than 50% of the increase in terrestrial evapotranspiration during the last thirty years. On the contrary, Betts et al. (2007) argue for a dominant global decrease of evapotranspiration induced by a water saving response to increasing carbon dioxide (CO2) concentrations. Recent modelling studies have found that ignoring the existence of this water saving response may overpredict future terrestrial evapotranspiration and drought (Milly and Dunne, 2016;Prudhomme et al., 2014;Swann et al., 2016). Piao et al. (2007) argues 10 instead that this physiological effect is cancelled by simultaneous CO2-induced increases in plant growth and total canopy leaf area, and that changes in climate and land use are the dominant drivers of evapotranspiration changes.
Tree or plot experiments in temperate and boreal forests focusing on the effect on evapotranspiration and water yield relate generally to forest management. These experiments are generally of short duration due to the difficulty and the costs of 15 performing long-term controlled experiments in the field. As forest clear-cutting has been found to reduce evapotranspiration and increase runoff, reforestation or regrowth has resulted in an opposite effect (Andréassian, 2004;Sørensen et al., 2009). For long-term and large-scale studies of this effect, basin-scale hydrological assessments are otherwise required (Andréassian, 2004). Long-term availability of precipitation and runoff data can be used to estimate evapotranspiration changes by water balance over long periods and to study the combined effect of climatic and forest change on forest water use (Hasper et al., 20 2016). However, the Budyko framework (Budyko, 1974) which links water and energy availability to basin-scale hydrology, is required to further separate this combined effect (Zhang et al., 2001). The Budyko framework has already been applied to understand impacts on forest evapotranspiration in Germany by air quality (Renner et al., 2013) and in North America to differentiate the responses of forested basins to climate and human drivers (Creed et al., 2014;Jones et al., 2012;Wang and Hejazi, 2011;Williams et al., 2012). It has also been applied in China and at the global scale to study the hydrological effects 25 of reforestation programs and the forest controls on water partitioning (Fang et al., 2001;Huang et al., 2003;Li et al., 2016;Xu et al., 2014;Zhang et al., 2008b;Zhou et al., 2015).
Observation-based studies of tree water use responses to elevated CO2 are also limited to tree-or plot-scale experiments of a few years (with the exception of a 17-year study by Tor-ngern et al. (2015)). Although most of these short-term tree field 30 experiments have found a stomatal water-saving response that reduces stomatal conductance and transpiration (Ainsworth and Rogers, 2007;Assmann, 1999;Medlyn et al., 2001), the response is often small or absent in northern tree species. In general, stomatal conductance tends to be less sensitive to high CO2 concentrations in gymnosperms (e.g. conifers) than in angiosperms (Hasper et al., 2017;Medlyn et al., 1999). For example, the conifer species Scots pine has shown no response or small reduction in stomatal conductance under elevated CO2 (Kellomäki and Wang, 1996;Sigurdsson et al., 2002, respectively) and 35 experiments with Norway spruce have resulted in no leaf water saving response (Hasper et al., 2016;Sigurdsson et al., 2013). Instead, the deciduous species Silver birch tends to show leaf water saving responses (Rey and Jarvis, 1998).
Our main objective was to determine in Swedish forests the dominant forest effect(s) driving changes in the ratio of actual evapotranspiration to precipitation, known as the evaporative ratio, from a basin scale approach. This, since 1900, forests in 40 Sweden have seen important changes in forest structure and composition (Antonson and Jansson, 2011). For this purpose, we used a former application of the Budyko framework Destouni, 2014, 2015;van der Velde et al., 2014) to study movement of basins in the space comprising the aridity index and the evaporative ratio, known as the Budyko space. We separated the climatic effect on evapotranspiration that relates to changes in the aridity index from a residual effect that relates to other drivers of change. We further explored the relationships between this residual effect and the forest attributes of standing 45 forest biomass, forest cover area or species composition. We chose the period 1961-2012 for analysing the change due to the large availability of runoff and forest attribute data across Sweden during this period.

Hydroclimatic data
We gathered data on the daily runoff (R) for 65 unregulated Swedish basins (Fig. 1) monitored by the Swedish Hydrologic and Meteorological Institute (SMHI) and compiled by Arheimer and Lindström (2015). The daily R data for the 65 basin outlets and corresponding basin boundaries were obtained from SMHI's hydrologic server (https://vattenwebb.smhi.se/station/) for 5 the period 1961-2012. The daily R data was aggregated to annual values and we considered only complete years in the analysis, defined as years with at least 98% data capture. We obtained the corresponding daily precipitation (P) data for the basins selected from 68 climatic stations from SMHI's online server. In order to deal with precipitation uncertainty, we also calculated daily P from the Climatic Research Unit-gridded P product CRU TS3.23 (Harris et al., 2014) and an additional P product of the Luftwebb portal (http://luftwebb.smhi.se/) of SMHI. The latter is a bias-corrected gridded dataset of precipitation for 10 Sweden during the period 1961-2014 with a 4 x 4 km horizontal resolution and based in data collected from over 87 precipitation stations around the country (Johansson, 2000;Johansson and Chen, 2003). We made a spatial interpolation by Thiessen polygon to obtain three different estimates of mean P for each basin and used the annual P and R data to calculate actual evapotranspiration (E) as where dS/dt is the change in water storage within the basin. We performed our assessments of hydroclimatic change during the period 1961-2012 by calculating the inter-annual means of the two 26-year sub-periods 1961-1986 and 1987-2012 and defined all hydroclimatic and forest changes as the difference between these means. Using such long-term periods becomes an advantage since when averaged over long periods dS/dt should be considerably smaller than P, R and E, permitting the basin-scale assumption of essentially zero long-term water storage change. These assumptions are required in order to calculate 20 E by Eq. (1).
We also obtained mean, minimum and maximum daily temperature data (T, Tmin, Tmax) from the climatic stations with P data for calculations of potential evapotranspiration (E0). Estimates of annual E0 in the 68 stations were obtained by three models: 1) The Langbein model in terms of T (Langbein, 1949), 2) the Hargreaves model in terms of Tmin and Tmax (Hargreaves et al., 1985) and 3) the gridded E0 product from the Climatic Research Unit CRU TS3.23 (Harris et al., 2014). For the first two 25 models, the E0 station estimates were interpolated spatially to each basin area by a trivariate spline which allows a spatially varying relationship between E0 and elevation (Tait and Woods, 2007). We obtained a mean E0 from the three estimates to perform an uncertainty assessment and to minimize the potential problems of relying on a single model.

Budyko framework
We used the Budyko framework (Budyko, 1974) to characterize the observed partitioning of P into R and E during the period 30 1961-2012. The Budyko framework is based on the relationship of the partitioning of water and energy on land and states that evapotranspiration is limited by the supply of water (i.e., P) and energy (i.e., E0). This relationship (Fig. 2a) is often represented in the space (i.e., Budyko space) created between the ratio of actual evapotranspiration to precipitation known as the evaporative ratio (E/P) and the ratio of potential evapotranspiration to precipitation known as the aridity index (E0/P).
The relationship between E/P and E0/P is represented by Budyko-type curves expressing the first in terms of the latter; E/P= 35 F(E0/P). The E/P is expected to increase linearly with E0/P at low values of E0/P, while gradually reducing its increasing rate at higher E0/P values. Here, we use the "Budyko-type" formulation of Yang et al. (2008) which is a climatic model of the evaporative ratio, (E/P)c, in terms of E0/P and a parameter representing the effect of the contribution of catchment characteristics, such as vegetation, soils and topography for each basin (n). This formulation has the same functional form as that of Choudhury (1999) and has been synthetized by Zhang et al. (2015) as 40 For this study, a mean n value was calculated for each basin by solving Eq. (2) with the mean E0/P and E/P values of the initial period 1961-1986, as done by Wang and Hejazi (2011). A basin that changes hydroclimatic conditions from a subperiod 1 (t1) to a subperiod 2 (t2) can be represented in Budyko space by a point moving from initial conditions (t1: E0/P1, E/P1) (Fig. 2). If the movement is only due to change in the aridity index (E0/P)), the movement will occur along the corresponding Budykotype curve to a new point (t2*: E0/P2*, E/P2*), implying a change in the evaporative ratio, termed henceforth as the climatic 5 effect ((E/P)c). However, a basin will most certainly move in time due to combination of change in the aridity index and other drivers of change Jaramillo and Destouni, 2014). Under such realistic circumstances, the basin will move to a new location not falling on the initial Budyko-type curve (t2: E0/P 2, E/P 2) implying a corresponding observed total change in the evaporative ratio ((E/P)).
Regarding some possible forest-related effects on evapotranspiration, increasing forest biomass should increase E due to an 10 increase in root depth and total canopy leaf area, which, under constant conditions of precipitation, would imply an upward movement in Budyko space (Fig. 2b). In a similar way, a tree water saving response to increasing atmospheric CO2 concentrations reduces stomata conductance and should decrease E, implying a downward movement in Budyko space under constant P.
Following van der Velde et al. (2014), we represented movement due to changes in the aridity index as a vector (v*) with 15 direction of movement (θ) and magnitude (r) calculated as where b is a constant and θ is in degrees (0º<θ<360º) starting clockwise from the upper vertical (b = 90° when E0/P) >0 and 20 b = 270° when E0/P)<0). In a similar way, we represented the total movement vector (v) based on R and P observations by replacing (E/P)c with (E/P)in Eqs. 3 and 4.
We used the approach of Jaramillo and Destouni (2014) to synthetize climate-driven and total movements for the 65 basins as typical wind roses of direction and magnitude. The roses illustrate in the same plot the combination of changes in the aridity index and evaporative ratios for all basins. It is worth noting that since the magnitude and direction of change in Budyko space 25 depends on the specific hydroclimatic conditions of each basin Gudmundsson et al., 2016;Jaramillo and Destouni, 2014) and on the definition of the space in which such movement occurs, such wind roses oversimplify the variability of movements in Budyko space. However, using these wind roses is a simple way to synthetize general tendencies of movement in large sets of basins and enables a first-order identification of the importance of drivers of change different from the aridity index. For instance, directions of movement that go beyond the range of slopes of any Budyko-type curve (i.e., 45º<θ<90º and 30 225º<θ<270º) imply that drivers different from the aridity index dominate the observed changes in the evaporative ratio. This when the n parameter in Eq.
(2) is held constant in time.

Forest attributes
According to recent data from Swedish National Inventory (NFI; http://www.slu.se/en/Collaborative-Centres-and-Projects/the-swedish-national-forest-inventory/), productive forestland now accounts for 57% of the Swedish land surface after 35 constantly expanding throughout the 20 th century (Nilsson et al., 2016). The productive forest standing volume has increased by 85% since the first NFI took place in 1923 (KSLA, 2015). In Sweden, forest management is responsible for most changes in the abundance of coniferous and deciduous species (Elmhagen et al., 2015;Laudon et al., 2011). The forest data from the NFI here used quantifies the surface area of all types of land use (i.e., productive forestland, mires, mountainous regions, agricultural land, rock surfaces and urban areas); however, the most comprehensive data is collected for forestland. The NFI 40 utilizes a stratified systematic sample based upon clustered sample plots designed to deliver statistics at the county level and as of today accounts for changes in methodology across time. The standing forest volume is differentiated into several categories (i.e., species, diameter and age composition, forest management stage). The species differentiated are Norway spruce (Picea abies), Scots pine (Pinus sylvestris L.), Silver birch (Betula pendula) and other deciduous broadleaf species. The strata of the NFI are the Swedish counties; the sample plots have been distributed within each stratum. A single sample distribution is completed every five years, however as each year (representing a fifth of the sample) is evenly distributed over the country, any consecutive five year period can be used. A detailed description of the inventory design and methods was provided by Fridman et al. (2014).
The main forest attributes studied here were forest cover area and standing volume of productive forest, extracted for the period 5 1961-2012 from the plots falling within the 65 basins. The sample plots are not restricted to provide only county wise estimates every sample plot has an upscaling factor that can be used to create estimates of forest attributes for a given area after grouping the sample plots falling in it. However, a larger area will result in more sample plots being used and a lower standard error. Many of the original 65 basins contained too few sample plots to provide meaningful estimates and therefore the basins were aggregated into larger groups. We thus aggregated the 65 basins into two main basin groups (i.e., temperate and boreal) 10 according to a terrestrial ecoregion classification (Olson et al., 2001) to show the change in forest statistics within the 65 basins and reduce the sampling standard error.
For the analysis, forest cover and standing forest volume (i.e. biomass) were divided by the area of each of the two basin groups to obtain relative values of area (A in km 2 /km 2 ) and volume (V in km 3 /km 2 ), respectively. According to the NFI data, in both basin groups more than 50% of the area is covered by forests (Fig. 1b), and these forests consist mainly of coniferous 15 species (Fig. 1c). Although both the boreal and temperate basin groups have a similar proportion of coniferous and deciduous species, species distribution varies between them in terms of standing volume. In the boreal group, Scots pine accounts for 50% of all trees where as in the temperate group Norway spruce is the dominating species (55% of all trees). Furthermore, as a proxy for species composition, we studied the ratio of the deciduous standing forest volume, Vd, to the total standing volume (Qv=Vd/V). We assumed a conservative sampling standard error of A (2.7%), V (5%) following (Fridman et al., 2014) and a 20 corresponding propagated error for Qv (10%) (Taylor, 1996).

Calculation of the residual effect on the evaporative ratio
In order to identify possible forest-related effects on the evaporative ratio, we separated the non-climatic residual effect from the climatic effect. The residual effect on the evaporative ratio from 1961-1986 to 1987-2012 ((E/P)r) was calculated as ∆( / ) r = ∆( / ) − ∆( / ) c (6) 25 The combination of three P products with the three E0 products resulted in nine possible combinations of P and E0 for each basin. This is important since the evaporative ratio estimates depend on the sources of P and E0 and the models being used (Greve et al., 2015;Wang et al., 2015). For each basin, we built the corresponding uncertainty ranges for n, (E/P)c, E/P, (E/P) and the components (E/P)c and (E/P)r. The arithmetic and the area-weighted (i.e., based on basin areas) means of (E/P), (E/P)c and (E/P)r for each basin group, temperate and boreal, were calculated by averaging the values of all basins within 30 that group that represent the largest area, that is, excluding any nested basins. In similar way, the uncertainty ranges of each basin were also aggregated into basin group uncertainty ranges.

Regression analysis of the residual effect and forest attributes
In order to determine which of the forest attributes (A, V, QV) could better explain the changes in the residual effect on the evaporative ratio, we performed for each basin group linear regressions between all annual values (52 values) of the residual 35 of the evaporative ratio, (E/P)r=E/P-(E/P)c, and each of the attributes A, V, Qv. A bootstrapping procedure was performed to account for the sampling standard errors of the forest attribute data in the regression analysis. We constructed a normal distribution of twenty data points for each annual value of each forest attribute with the given mean and sampling standard error of each attribute data point. The coefficient of determination (R 2 ) and the statistical significance of the linear regression (p<0.05) were calculated for each of 1000 random samples of the annual time series 1961-2012. We finally obtained the mean 40 of the 1000 combinations of R 2 and p-values to determine the attributes that could best explain the residual effect.
The assumption of steady state conditions (dS/dt0) of the Budyko framework required to maintain the constant water supply limit may also become a potential source of uncertainty for calculations of E/P)r (Chen et al., 2013;Greve et al., 2016;Wang 6 and Tang, 2014;Zhang et al., 2008a). It is possible that a time interval of one year may not guarantee stationary conditions (i.e., dS/dt≠0) in these basins, affecting the calculations of E (Eq. 1) (Bring et al., 2015;Budyko, 1974;Moussa and Lhomme, 2016). We therefore applied the same bootstrapping and regression analysis to the 17 values and 10 mean values resulting from averaging the annual data of the period 1961-2012 into three-and five-year means, respectively.
Additionally, recent findings have shown that changes in the ratio of precipitation in the form of snow to total precipitation (fs) 5 might affect E/P (Berghuijs et al., 2014). We therefore additionally calculated fs per year and for each basin based on the mentioned collected daily P and T data; we assumed that in days with mean temperatures below 1 ͦ C precipitation fell as snow and above 1 ͦ C as rain, following Berghuijs et al. (2014). Finally, we spatially aggregated these values to the basin-group level and incorporated fs into the regression analysis previously described.

Movement in Budyko space and separation of effects
Most of the 65 basins presented energy-limited conditions since their aridity index E0/P fell below one (Fig. 3). This means that evapotranspiration in these basins was more limited by energy than by water availability. The mean E0/P and evaporative ratios E/P of the boreal basins were generally lower than those of the temperate basins due to their northerly and cool location. Some basins plotted low in Budyko space, i.e., evaporative ratio E/P near zero, possibly due to underestimates of precipitation 15 due to precipitation undercatch in rain gauges because of falling snow and/or wind, unrealistic measurements of runoff or significant groundwater flux across the basin boundary.
The roses of movements in Budyko space (Fig. 4) show the direction and magnitude of movement of each all basins from the period 1961-1986 to the period 1987-2012 in the same way that a typical wind rose shows wind direction and wind speed. The roses showed that from the period 1961-1986 to the period 1987-2012, all basins in temperate and boreal biomes 20 experienced a decrease in the aridity index (180º<θ<360º; Fig. 4a, c). The roses also evidenced an increase in the evaporative ratio in most basins in both the temperate (60%) and boreal (61%) groups, moving them upwards in Budyko space (Fig. 4a, c), with the directions of movement in both biomes being much similar. Furthermore, the movements of highest magnitude occurred in upward directions, with increasing E/P. However, in the absence of other drivers of change, a decrease in the E0/P can only result in a decrease in the E/P (Fig 4b, d). Hence, the occurrence of the upward movements (270º<θ<360º) suggests 25 the influence of other drivers of that movement that are from E0/P. As such, changes in the aridity index are not the only or the most important drivers in these basins.
The general decreases in the E0/P appeared to be the result of an increase in P between the two time periods that was considerably larger than an increase in E0 (Fig. 5). The increase in P occurred mostly around winter (January and February) 30 and late springsummer (May -August) and was larger in the temperate than in the boreal group of basins. For instance, the maximum increase of P in the boreal group occurred in June and reached 25 mm/yr. On the contrary, although the increase in E0 was larger in the temperate basin group, it was rather small with the highest value occurring in April in both biomes (7 mm/yr across biomes).

The residual effect on the evaporative ratio
The separation of the climatic and residual effects on the evaporative ratio at the basin group scale shows that the evaporative ratios in the temperate and boreal basin groups have experienced a decreasing climatic effect driven by less arid conditions, (E/P)c<0, between the periods 1961-1986and 1987. They have also experienced an increasing residual effect due to other drivers of change, (E/P)r>0. Note that the distribution of change for all basins shows that these counteracting 40 effect applies to the median, arithmetic and area-weighted means and interquartile ranges of both basin groups. The uncertainty assessment of the arithmetic and area-weighted means of (E/P)c and (E/P)r demonstrates that regardless of the data sources and models of P, T, E0 used, (E/P)c and (E/P)r always show these counteracting effects in both basins groups.

Forest change and effects on the evaporative ratio (E/P)r
The positive residual effect (E/P)rin both biomes might be the result of increasing E due to increasing standing forest volume or surface cover. Indeed, between the two 26-year periods, the forest cover area A and the standing forest biomass V increased 5 beyond the standard sampling error in the group of temperate basins (Fig. 7). In the group of boreal basins, V also increased more than the standard sampling error but A remained stable. There was no statistically significant change in the overall forest composition in any of the basin groups; the change in Qv was considerably smaller than the propagated sampling standard error.
The regression analysis between the annual (E/P)r and V in the period 1961-2012 presented a significant (p<0.05) positive 10 relationship in the boreal basin group when using any of the one-year, three-year and five-year mean options used to divide the period 1961-2012 ( Fig. 8 and Table S1 of the Supplementary materials). Hence, V was the attribute that could best explain the residual effect on the evaporative ratio in this basin group. The A, Qv and fs could not explain significantly (p>0.05) the variance of (E/P)r in any case, although fs was close to be significant (p=0.06) based on data with a three-year time step. The V was also the only attribute that could significantly (p<0.05) explain the variance of (E/P)r in the temperate basin group, but 15 only when using annual data. The relationship of (E/P)r and V was almost significant (p=0.07) for three-year mean data and that of (E/P)r and A when using annual data (p=0.08). Hence, among the four parameters studied, V was the attribute that could best explain the residual effect on the evaporative ratio also in the temperate group.

Dominant forest-related effects on the evaporative ratio
Our results show that the increase in the residual effect on the evaporative ratio E/P)r during the period 1961-2012 was best explained by an increase in standing forest biomass. Under a hypothetical situation with no changes in climate and stomatal conductance, increasing standing forest biomass should increase transpiration if the canopy leaf area index LAI and root biomass increased along. At least in the case of LAI, an increase in Sweden has been indicated by remote sensing and modelled 25 data (Zhu et al., 2016), and LAI and standing forest biomass are typically strongly correlated in Scandinavian forests (Heiskanen, 2006). Studies in many forest ecosystems across the world have found with reforestation a similar increase in evapotranspiration and/or decrease in runoff at the plot scale (Andréassian, 2004;Bosch and Hewlett, 1982;Sørensen et al., 2009) and basin scale (Fang et al., 2001;Huang et al., 2003;Li et al., 2016;Xu et al., 2014;Zhang et al., 2008b).

No evidence of a water saving response-effect to rising atmospheric CO2 concentrations 30
A possible water-saving response caused by reduced stomatal conductance under rising atmospheric CO2 concentrations (Ainsworth and Long, 2005) should have instead decreased E/P)r. However, our results show that in both basin groups E/P)r rather increased. Hence, the increase in E/P)r found in the two basin groups between the two periods implies that this water saving response is either negligible or overridden by the positive effect from increasing standing forest biomass. The trees in both the temperate and boreal basin groups are dominated (>80%) by coniferous tree species, and stomatal conductance tends 35 to be less responsive to increased CO2 concentrations in such species than in deciduous species (Hasper et al., 2017;Medlyn et al., 1999). This has indeed been shown in controlled CO2 experiments with trees of the two species dominating our basins, Scots pine and Norway spruce (Kellomäki and Wang, 1996;Sigurdsson et al., 2002;Hasper et al., 2016). Furthermore, there was no statistically significant change in forest composition in the two basin groups towards more deciduous trees with more CO2 responsive stomata that could have affected E/P)r. 40

Intra-annual hydroclimatic changes and their effect on the evaporative ratio
Intra-annual or seasonal changes of precipitation and temperature with no corresponding change in the annual absolute values of these variables may also affect E/P)r (Chen et al., 2013;Milly, 1993;Zanardo et al., 2012). However, the general seasonality pattern of potential evapotranspiration and precipitation in Northern latitudes is persistent in time. A decreasing fraction of precipitation falling as snow has been also related to an increase in E/P)r (Berghuijs et al., 2014). Indeed, our calculations 5 show that the fraction of precipitation falling as snow decreased from the period 1961-1986 to the period 1987-2012 from 0.20 to 0.14 in the temperate basin group and from 0.45 to 0.43 in the boreal basin group. In the case of the boreal basin group, this fraction may partly explain a small variance of E/P)r when performing the regression with three-year intervals. Studies in China have found changes in the evaporative ratio with decreasing snow fraction after assuming that in such latitudes, spring snowmelt does not infiltrate the soil and mainly flows on land as runoff, not being used by plants for their biological processes 10 (Zhang et al., 2015). Nevertheless, this extreme assumption does not apply to the conditions of Swedish forests. Laudon et al. (2004Laudon et al. ( , 2007 found that a majority of the snowmelt water infiltrates the soil in boreal Sweden and is available for plants in spring. It is therefore possible that the apparent effect of the decrease in spring snow on E/P)r in boreal basins is rather a consequence of a positive effect of longer growing season (and thus higher annual transpiration) in years with lower snow fraction. 15

Uncertainty on other possible drivers of change on the evaporative ratio
Apart from the inherent uncertainty in the hydroclimatic and forest-attribute data, are there any other factors that may explain the increase in E/P)r in both basin groups? The areal extent of forest fires and wind-throws in Sweden is small, and such effects should thus have minor influences on large-scale evapotranspiration. Drainage practices used in forestry and agriculture may have played a role in the partitioning of water in Swedish landscapes, although the magnitude of this effect is uncertain 20 and differs between different drainage techniques (Wesström et al., 2003). Farmers implemented open ditches in forests throughout the 20 th century, mainly during the 1930s, before the present study period, and the extent of functioning drainage in Sweden today is about 5% of the total land area (Berglund, Ö, et al., 2009).
Flow regulation and hydropower development have been shown to increase the evaporative ratio at local to global scales Jaramillo and Destouni, 2015;Levi et al., 2015). However, all the basins used in this study are 25 unregulated (Arheimer and Lindström, 2015). Technological improvements over time in rain gauges to reduce losses of water due to wind, evaporation and snow undercatch may also reflect as a false increase in precipitation and corresponding increase evapotranspiration. Nevertheless, in Sweden, these improvements occurred before the period of the present study; windshields were introduced in Sweden in 1853, and the last installations were reported in 1935. Even the replacement of the old cans made out of galvanized iron for new aluminum ones occurred around 1958 (Sofokleous, 2016). 30

Conclusion
We used the Budyko framework to study movements in Budyko space of 65 unregulated Swedish basins, aiming to detect the dominating effect(s) of forest. In all basins, the aridity index decreased due to an increase in precipitation larger than a corresponding increase in potential evapotranspiration. In 60% of the basins, this decrease was accompanied by an increase in the evaporative ratio, which as the Budyko framework of water and energy availability implies, cannot be explained by changes 35 in precipitation or potential evapotranspiration. In both the temperate and boreal basin groups studied, a positive residual effect on the evaporative ratio counteracted the negative climatic effect, maintaining the evaporative ratio relatively stable over time.
The positive residual effect occurred along with increasing standing forest biomass in both the temperate and boreal basin groups, increasing forest cover in the temperate basin group, and with no apparent changes in forest species composition. In general, standing forest biomass was the forest attribute that best explained the variance of the residual effect on the evaporative 40 ratio across basin groups. These results further suggest that the water saving response to increasing CO2 in these forests or any other negative effect is either negligible or overridden by the opposite effect of the increasing forest biomass. We conclude that forest expansion is the dominant driver of long-term and large-scale evapotranspiration changes in Swedish forests.  Evaporative ratio (E/P) vs. aridity index (E0/P) and the parameters describing movement in such space from period 1 (t1) to period 2 (t1) by a vector (v). The t1 and t2 represent mean conditions of E/Pand E0/Pduring the periods 1961-1986 and 1987-2012, respectively. The change in evaporative ratio (E/P is divided into an effect from changes in the aridity index, E/P)c,termed the climatic effect and represented by the vector (v*), and an effect from other drivers, termed the residual effect, E/P)r. The elliptic curve describes the Budyko-shaped curve of Chodbhurys's relationship between E/Pand E0/P. (b)

5
The expected directions of movements in Budyko space with fixed E0/P for increasing forest standing biomass and water saving response to increasing atmospheric CO2 concentrations under constant precipitation illustrated in Budyko space, in terms of the aridity index (E0/P; x-axis) and evaporative ratio (E/P; y-axis) for temperate (purple) and boreal (green) basins.
5 Figure 4: Hydroclimatic movement in Budyko space for the 65 basins due to changes in aridity index (E0/P) and evaporative ratio (E/P) between the two comparative periods 1961-1986 and 1987-2012. Wind roses of individual basin movements of the boreal basingroup according to: (a) the combined effect (green) of all drivers of change by calculating E/P from runoff observations (i.e., E/P by Eq. (2)) and (b) the effect (grey) of only change in the aridity index by calculating E/P from modelled data (i.e., (E/P)c by Eq. (4)), with fixed n and varying E0/P). (c-d) Similar wind roses for temperate basins (purple and grey, respectively). The range of directions of movement 5 (0<θ<360°) is divided into 15° interval-paddles that group all basins moving in each direction interval, with directions θ starting from the upper vertical and clockwise. We chose this 15°-value arbitrarily, to provide the sufficient detail of directions. Intensity of colour intervals represent the magnitude of the movements (r) in Budyko space in such given direction θ. As example, 37% of all boreal basins (Fig. 4a) moved in the range of directions 270°<θ<295°. Now, the intensity of the color describes the range of magnitudes (r: -) of those movements in Budyko space. As example, of those boreal basins moving in the range of directions previously described (37% of the boreal basins), 14% 10 have moved with magnitudes between 0 and 0.05 (light green), 14% with magnitudes between 0.05 and 0.10 (medium green) and 9% with magnitudes between 0.10 and 0.30 (dark green). On the contrary, no boreal basin moved in this range of directions (270°<θ<295°) when using the Chodbhury Equation Fig. 4b).
15 Figure 5: Intra-annual change in precipitation (P) and potential evapotranspiration (E0). Changes in the mean of P and E0 between the periods 1961-1986 and 1987-2012 in the boreal (a-b) and temperate (c-d) basin groups for each month (January (1) to December (12)).  1961-1986 to 1987-2012 in the evaporative ratio, ΔE/P, and its climatic, Δ(E/P)c, and residual, Δ(E/P)r, effect-components in the (a) boreal and (b) temperate groups of basins. Boxplot statistics include arithmetic mean (blue triangles), area-weighted mean (pink circles), median (thick horizontal black line), interquartile range (IQR) (boxes), whiskers (confidence interval of ±1.58*IQR√ where N is the number of basins in each biome group) and outliers (small black circles). Corresponding uncertainty ranges of the calculated arithmetic and area-weighted means of ΔE/P, Δ(E/P)c and Δ(E/P)r for 5 each basin group are shown as blue and pink vertical ranges of uncertainty, respectively. These uncertainty ranges arise from the use of three different precipitation and potential evapotranspiration data products.  1961-1986 and 1987-2012, in: a) forest cover (A), b) standing forest volume (V) and c) forest composition represented as the ratio of the standing deciduous-forest volume to total standing forest volume (QvThe whiskers represent the uncertainty range of A, V and Qv corresponding to the sampling standard errors for A (2.7%), V (5%) and propagated Qv (10%). For (b), the dark blue represents V of deciduous trees and light blue of coniferous trees.
5 Figure 8: Performance of forest attributes as predictors of the residual of the evaporative ratio. Coefficients of determination (R 2 ) of the linear regressions between the residual of the evaporative ratio, (E/P)r, and the forest attributes of forest biomass (V), forest cover (A), forest composition (Qv) and fraction of precipitation falling as snow (fs), in the boreal and temperate basin groups. The Pearson coefficient of determination (R 2 ) is shown for the regressions using annual (N=52), three-year (N=17) and five-year means (N=10) during the period 1961-2012, with N being the number of values available for each regression with each interval selection. The significance (p<0.05) of the 5 R 2 of each linear regression is shown with a black dot. See Table S1 in Supplementary materials for values.