Manuscript Number: EGUSPHERE-2025-4820

Revision notes: We sincerely thank the Editor for recognizing the potential value of our manuscript and for the time and effort devoted to handling the review process, including coordinating the assessment and inviting expert reviewers. We also sincerely thank both reviewers for their careful reading of our manuscript and for providing constructive and insightful comments. These comments have helped us identify several aspects that require further clarification and improvement. In the following response, the reviewers’ comments are shown in black, and our point-by-point responses and planned revisions are shown in blue. At this final-response stage, we provide a detailed explanation of how we will revise the manuscript according to the Editor’s and reviewers’ suggestions. We will carefully improve the Methods, Results, Discussion, figures, captions, and uncertainty analysis to make the manuscript more rigorous, transparent, and internally consistent. We respectfully hope to be given the opportunity to revise the manuscript and incorporate the suggested modifications.

Responses to the reviewer’s comments

Response to Reviewer #1

We sincerely thank the reviewer for carefully reading our manuscript and for providing constructive and detailed comments. We fully recognize that several aspects of the manuscript require further clarification, particularly the consistency of the nonlinear threshold values, the methodological assumptions of the moving-window approach, the unit and interpretation of dE/dVPD, the description of interannual VPD variability, the scale and interpretation of the SEM analysis, and the uncertainty associated with soil moisture and evapotranspiration products. In the responses below, we provide a point-by-point plan for how we will revise the manuscript according to the reviewer’s suggestions. We will correct the inconsistent terminology and units, clarify the spatial scale and assumptions of the analysis, revise unsupported or over-causal interpretations, and strengthen the uncertainty discussion. We sincerely appreciate the reviewer’s comments and respectfully hope to be given the opportunity to revise the manuscript and incorporate these suggested modifications.

Comment 1: Lines 27–29: The climatic-gradient thresholds are reported inconsistently. In the Abstract, the arid-zone threshold is given as 1.31 kPa, whereas in Section 4.1 and the Conclusion (Line 468) it is given as 1.67–1.68 kPa. The same discrepancy applies to the temperate zone (0.68 kPa in the Abstract; 1.67–1.68 kPa referred to in Section 4.3, Line 407). The authors must reconcile these values, clarify whether the Abstract reports the GAM-derived or piecewise-derived threshold, and ensure consistency throughout.

Response: Thank you for pointing out this inconsistency. We agree that the threshold values were not reported in a clear and consistent way in the original manuscript. In the revision, we will use one clear aridity-based framework throughout the manuscript. We will divide global land areas into four aridity zones: arid, semi-arid, semi-humid, and humid regions. We will no longer mix different climate grouping methods when reporting threshold values. This change will make the threshold comparison clearer. We will also clarify how the threshold values are obtained. We will not assume in advance that the one-breakpoint piecewise model is the only suitable model. Instead, we will compare three models: the one-breakpoint piecewise model, the two-breakpoint piecewise model, and the GAM. We will use model performance metrics, including model fit and error, to judge which model better describes the VPD–E relationship in each aridity zone. The reported thresholds will then be based on the model that gives the most suitable description of the nonlinear response. Under this revised framework, we will report one consistent set of threshold values throughout the Abstract, Results, Discussion, and Conclusion. The threshold values will be 1.90 kPa for the arid zone, 1.46 kPa for the semi-arid zone, 0.49 kPa for the semi-humid zone, and 0.47 kPa for the humid zone. We will remove the earlier inconsistent values, including 1.31 kPa, 0.68 kPa, and 1.67–1.68 kPa, from the relevant parts of the manuscript. We will also revise the wording around these values. We will describe them as estimates of the main transition in the observed VPD–E relationship. We will avoid presenting them as exact ecological limits. This wording is more appropriate because the threshold estimates depend on the model structure and the input data.

These changes will make the threshold values and their methodological basis consistent across the Abstract, Methods, Results, Discussion, and Conclusion.

Comment 2: Section 2.2.2: The "space-for-time" moving-window approach to compute dE/dVPD is borrowed from land-use-change literature, where the logic is well-established for discrete cover transitions. Applying it here to derive VPD sensitivity assumes that spatial VPD gradients at 5 × 5 km are driven solely by biophysical feedbacks, rather than by mesoscale atmospheric dynamics. Please provide justification for this assumption for the global domain, especially in complex terrain or coastline-adjacent pixels. Furthermore, the authors should discuss the robustness of this assumption and, if possible, offer a cross-validation against pixel-level temporal regression (dE/dVPD from time series at each pixel).

Response: We agree that the space-for-time moving-window approach needs a clearer methodological explanation. We also agree that this method should not be interpreted as direct causal evidence.

In the revision, we will first clarify the spatial scale of this method. The analysis will be conducted on a common 0.1° grid. The moving window will be defined as a 5 × 5 grid-cell neighborhood around each target pixel, rather than a 5 × 5 km window. This description will avoid confusion about the spatial resolution of the input datasets.

We will also add local screening rules within the moving window. Candidate pixels will be retained only when they share the same dominant MODIS land-cover type as the target pixel. The difference in the fraction of the dominant land-cover type will be required to be less than 10%. Candidate pixels with an elevation difference greater than 100 m from the target pixel will also be excluded. These rules will reduce the influence of strong land-cover differences and topographic gradients.

We agree that the moving-window approach cannot fully remove the influence of mesoscale atmospheric processes. Local VPD gradients may still be affected by topography, air movement, and land-sea contrast. This issue may be more important in mountainous and coastal regions. We will add this limitation to the manuscript and will state that the estimated dE/dVPD should be interpreted as an apparent local sensitivity, not as direct causal evidence.

Given the limitations of the current final-response stage, we will not overstate the validation here. In the revision, we will either provide a supplementary temporal-regression comparison where feasible, or explicitly describe the absence of such a comparison as a limitation and a priority for future testing. In addition, we will support the reliability of the input data by comparing the main ERA5-Land VPD and GLEAM E datasets with independent VPD and E products. This dataset-level comparison will help assess whether the main results depend strongly on a single data source.

These changes will make the method clearer and will make the interpretation more cautious.

Comment 3: Line 164: The moving-window size is stated as "5 × 5 km based on previous studies," but no citation is provided, and the dataset native spatial resolutions (ERA5 ~31 km, GLEAM ~25 km, MODIS MCD12C1 ~5.5 km) are substantially coarser than 5 km. This creates a fundamental mismatch: in practice, a 5 × 5 km window over ERA5 or GLEAM data contains a single pixel. The authors need to explicitly state the actual grid resolution at which the analysis was conducted and reconcile it with the stated window size.

Response: Thank you for pointing out this problem. We agree that the phrase “5 × 5 km” was misleading. This wording does not match the spatial resolution of the datasets used in this study.

In the revision, we will correct this description. We will state that the analysis is conducted on a common 0.1° grid. This grid is used for both the ERA5-Land VPD data and the GLEAM E data in our analysis. Therefore, the moving window is not a 5 × 5 km window. It is a 5 × 5 grid-cell window. This means that each target pixel is compared with candidate pixels within a 0.5° × 0.5° neighborhood. We will describe this clearly in Section 2.2.2. This correction will remove the mismatch between the stated window size and the actual data resolution.

We will also explain how candidate pixels are selected within this window. To reduce the effect of surface heterogeneity, we will keep only neighboring pixels with the same dominant MODIS land-cover type as the target pixel. We will require the difference in dominant land-cover fraction to be less than 10%. We will also remove pixels with an elevation difference greater than 100 m from the target pixel.

These changes will make the spatial resolution, moving-window size, and sample selection rules clearer. They will also make the method consistent with the actual resolution of the input datasets.

Comment 4: Line 270: The global mean sensitivity of 293.27 ± 62.28 mm·hPa⁻¹·yr⁻¹ is a central quantitative result, yet the units are unconventional. Sensitivity (dE/dVPD) is dimensionally [mm / hPa], not [mm·hPa⁻¹·yr⁻¹]. The inclusion of yr⁻¹ suggests this may be derived from annual trend magnitudes rather than a pure partial derivative. Unless I have misunderstood the definition of what quantity this represents, e.g., is it a sensitivity coefficient, a trend-normalised ratio, or a regression slope?

Response: Thank you for this careful comment. We agree that the original unit of dE/dVPD was incorrect. In the revision, we will clarify that dE/dVPD is a sensitivity coefficient. It describes the change in E per unit change in VPD. Its correct unit is mm hPa⁻¹. It is not an annual trend rate, and it is not a trend-normalized ratio. We will therefore remove yr⁻¹ from this unit throughout the manuscript. We will also clarify how this sensitivity is estimated. In this study, dE/dVPD is calculated as the slope between E differences and VPD differences within the moving-window framework. The Theil–Sen slope estimator is used to reduce the influence of outliers. We will report the global mean sensitivity as 293.27 ± 62.28 mm hPa⁻¹, rather than 293.27 ± 62.28 mm hPa⁻¹ yr⁻¹. We will also revise the related text in the Results, figure captions, and Conclusion to avoid confusion between sensitivity and temporal trend.

Comment 5: Lines 296-303: The statement that "VPD rose abruptly in 1998 (0.034 hPa yr⁻¹)" reports a trend rate, not an abrupt change, and the phrasing conflates trend magnitude with event-driven anomaly. Additionally, Figure 3b shows VPD peaking at 7.839 hPa in 2010, the connection between the 1998 El Niño and the 2010 peak is not straightforward and warrants more rigorous treatment, potentially with ENSO-index partial correlation.

Response: Thank you for this careful and constructive comment. We agree with the reviewer that our original wording was not rigorous enough. The phrase “VPD rose abruptly in 1998 (0.034 hPa yr^-1)” was inappropriate, because 0.034 hPa yr^-1 is a trend rate rather than the magnitude of an abrupt event. This wording may cause confusion between a long-term trend and a short-term climate anomaly. We also agree that the VPD peak in 2010 should not be directly linked to the 1998 El Niño without a dedicated attribution analysis. Such a link would require additional tests, such as an ENSO-index partial correlation or other attribution methods. These analyses are beyond the current scope of this manuscript. Therefore, we will remove the statement about the abrupt increase in 1998 and will avoid linking the 2010 VPD peak to the 1998 El Niño.

In the revision, we will rewrite the relevant part of Section 3.3 in a more cautious way. We will only describe the observed long-term trend and interannual variability of VPD and E during 1981-2020. We will not make an event-based explanation unless it is directly supported by the analysis. We will also adjust the focus of this section. Instead of discussing ENSO-related interpretation, this section will focus on the relationships between VPD and the two E components, namely E_t and E_b. We will use partial correlation analysis to examine the relationships between VPD and E_t/E_b after controlling for precipitation and soil moisture. This change will make the interpretation more consistent with the evidence shown in Fig. 3. It will also reduce the risk of overinterpreting the time-series patterns. We sincerely appreciate this comment, because it helps us separate trend description from event attribution and makes the discussion more reliable.

Comment 6: Lines 362-364: The SEM model is stated to explain 77% of the variance in E, and this is presented as a single global metric. However, SEM was presumably fitted globally on spatially aggregated, or grid-cell mean data. The authors should clarify the unit of analysis (pixel-level, regional mean, or annual global mean) and include model fit statistics (Fisher’s C, AIC) in the main text rather than only referencing Figure 5.

Response: We agree that the original description of the SEM analysis was not clear enough. The statement that the SEM explained 77% of the variance in E may be misunderstood as a pixel-level or spatially explicit global result. This is not what we intended. In the revision, we will clearly state that the SEM is fitted using annual global-mean variables. The unit of analysis is therefore the annual global mean, not individual pixels or regional averages. We will also revise the wording around R² = 0.77. We will explain that this value refers only to the SEM fitted at the annual global-mean scale.

We will also add the main model fit statistics directly to the main text. We agree that these statistics should not only appear in Figure 5. In Section 4.2, we will report Fisher’s C, P , and AIC. We will also explain that these values are used to evaluate the overall fit of the SEM. This change will make the SEM description more transparent. It will also help readers understand the scale, meaning, and limits of the reported 77% explained variance. We sincerely appreciate this comment, because it helps us avoid overstating the SEM result and makes the interpretation more precise.

Comment 7: Lines 375-378: The finding that LAI increases markedly under high VPD (standardised coefficient = 0.90) but has only a weak direct effect on Et (standardised coefficient = 0.07) is counterintuitive and deserves deeper mechanistic discussion. A strong VPD–LAI relationship with near-zero LAI–Et effect suggests that either collinearity in the SEM is absorbing the LAI pathway into the VPD direct effect, or that the LAI variable used (GIMMS V1.2 monthly) does not resolve the sub-seasonal dynamics at which VPD–stomatal coupling operates. Please, make sure to address this explicitly.

Response: The combination of a strong VPD–LAI path and a weak direct LAI–E path may look counterintuitive. We will revise Section 4.2 to explain this point more carefully. We will make clear that the weak direct LAI effect does not mean that vegetation is unimportant for E. Instead, it may show that the vegetation signal is partly shared with VPD and SM in the SEM. In this case, part of the LAI-related effect may be absorbed by the direct VPD pathway or by the SM-related pathway. We will therefore describe this result as a statistical pathway under multi-factor coupling, not as direct evidence that high VPD promotes vegetation growth.

We will also add a clear limitation about the LAI data and the temporal scale of the analysis. The SEM uses annual global-mean variables, while the VPD–stomatal response can occur at much shorter time scales. The GIMMS LAI data may not capture these short-term vegetation and stomatal changes. This mismatch may weaken the direct LAI–E pathway in the SEM. In the revised text, we will state that VPD regulates E mainly through atmospheric moisture demand and soil moisture limitation at the analysis scale used here. We will also state that the LAI pathway should be interpreted with caution because of possible collinearity and scale mismatch. We sincerely appreciate this comment, because it helps us avoid an over-strong mechanistic interpretation of the SEM result.

Comment 8: Lines 406-410: In the discussion of the SM–atmospheric water vapor feedback under high VPD in arid regions, the mechanism is described qualitatively without connecting it to the quantitative results established earlier (e.g., the SEM path coefficients or the dE/dVPD spatial distribution). The discussion would be very interesting if the authors could map the regions where VPD already exceeds the 1.67-1.68 kPa threshold (as identified in Section 4.1) onto the negative-sensitivity regions shown in Figure 2a to demonstrate that the feedback mechanism is already active.

Response: We agree that the previous discussion was too qualitative. The link between the proposed soil moisture–atmospheric water vapor feedback and our quantitative results was not clear enough. In the revision, we will rewrite this part of Section 4.3. We will connect the discussion more directly to the dE/dVPD spatial pattern, the aridity-zone thresholds, and the SEM path coefficients. In particular, we will discuss the strong negative VPD–SM pathway in the SEM and the negative dE/dVPD regions shown in Fig. 2a. This change will help show how high atmospheric dryness may weaken the positive response of E to VPD when water supply becomes limiting. We will also adjust the threshold discussion to match the revised threshold framework. The manuscript will use aridity-index-based zones rather than the earlier threshold values mentioned in the original text. We will therefore discuss the arid and semi-arid thresholds as 1.90 kPa and 1.46 kPa, respectively. If the revised analysis confirms a clear spatial overlap, we will add an overlay analysis between areas exceeding the relevant VPD threshold and areas with negative dE/dVPD. This analysis will be used to test whether high-threshold regions are already linked with suppressed E under rising VPD. If the overlap is not strong or not spatially robust, we will not overstate this mechanism. Instead, we will describe it as a plausible feedback pathway supported by the SEM and spatial sensitivity results. This revision will make the discussion more quantitative and more cautious.

Comment 9: Lines 431-436: The acknowledgement that ERA5 overestimates SM in arid regions (Kokkalis et al., 2024) is important, but its consequence is underplayed. For example, given that the SEM assigns the VPD–SM path a standardized coefficient of -0.84, the strongest single path in the entire model, a systematic positive bias in SM would directly attenuate this coefficient, inflating the apparent direct effect of VPD on E. This has a directionally predictable effect on the ms's central conclusions and should be discussed with explicit reference to the SEM coefficients.

Response: We agree that the implication of ERA5-Land soil moisture uncertainty was underdeveloped in the original manuscript. This issue is especially relevant to our SEM interpretation, because the VPD–SM pathway has the strongest standardized coefficient in the model (−0.84). In the revision, we will explicitly discuss this point rather than only mentioning ERA5-Land SM bias in general terms.

We will clarify that possible overestimation of SM in arid regions may affect the estimated balance between the direct and indirect effects of VPD on E. If the SM product overestimates moisture availability, or if it weakens the representation of dry anomalies and soil moisture variability, the negative VPD–SM coupling may be underestimated. In that case, part of the SM-mediated indirect effect of VPD on E could be shifted into the apparent direct VPD–E pathway. This would make the direct effect of VPD appear stronger than it would be if soil moisture stress were fully captured. We will therefore state that the SEM result should be interpreted as a statistical pathway based on the available datasets, and that the relative magnitudes of the direct VPD effect and the SM-mediated indirect effect remain sensitive to SM uncertainty.

We will also strengthen the uncertainty discussion for E and its components. Current global ET products differ in long-term trends and in the partitioning of E into transpiration, soil evaporation, and interception evaporation. These differences may influence both the spatial dE/dVPD pattern and the SEM path coefficients. We will make clear that these uncertainties do not remove the evidence that VPD is closely linked to global E variability, but they do affect how confidently we can separate the direct atmospheric-demand pathway from the indirect soil-moisture pathway. We sincerely appreciate this comment, because it helps us present the SEM results with more appropriate caution and avoids overstating the central mechanism.

Comment 10: Lines 444-449: The authors acknowledge that piecewise regression assumes a single breakpoint. Yet in the Conclusion (Lines 466-469), a full table of single-threshold values per climate zone is presented as a primary finding. The authors should either present evidence that the VPD–E response in each zone is indeed well described by a single breakpoint (e.g., by showing that a two-breakpoint model does not significantly improve fit) or explicitly qualify the reported thresholds as approximations of the dominant transition, avoiding overstatement of their precision.

Response: We agree that the threshold values should not be presented as overly precise ecological limits. The original wording may have given the impression that the VPD–E response was assumed to have only one breakpoint in each zone. This was not our intention. In the revision, we will clarify that the thresholds are identified after comparing one-breakpoint piecewise regression, two-breakpoint piecewise regression, and GAM. The purpose of this comparison is to evaluate whether a dominant transition point can adequately describe the observed nonlinear VPD–E response, rather than to impose a single-breakpoint structure a priori.

We will revise Section 4.1 and the Conclusion to make this point explicit. We will state that the reported thresholds represent the dominant transition points supported by the model comparison, not exact or unique ecological boundaries. We will also explain that the two-breakpoint model does not provide a clear improvement over the one-breakpoint model based on the fit comparison, while the GAM results show broadly consistent transition behavior. The results from the GAM and two-breakpoint models will be presented as supplementary evidence to support the robustness of the main threshold interpretation. At the same time, we will soften the wording in the Conclusion and avoid overstating the precision of the threshold values. We sincerely appreciate this comment, because it helps us present the nonlinear threshold results in a more transparent and cautious way.

Additional planned revisions

In addition to the revisions made in response to the two reviewers’ comments, we also plan to make several additional changes to improve the rigor and clarity of the manuscript. First, we will remove the analysis of canopy interception from Section 3.3. The reason is that the mechanisms linking VPD to transpiration and soil evaporation are relatively well established, whereas the influence of VPD on canopy interception is less direct and more difficult to interpret mechanistically. Canopy interception is mainly controlled by canopy structure, commonly represented by LAI, and by precipitation amount and intensity. Therefore, a statistical relationship between VPD and interception may not necessarily indicate a causal mechanism. To avoid overinterpretation, we will focus Section 3.3 on transpiration and soil evaporation, for which the physical and ecohydrological interpretations are more robust. Second, we will detrend the data used in Section 3.3 before conducting the correlation analysis, so that the results better reflect interannual covariation rather than shared long-term trends. We will also add significance markers in Figure 3 to indicate correlations significant at P < 0.05. Third, we will add the official citation for ERA5-Land to improve the completeness of the data-source description. Finally, we will add supplementary maps showing the spatial distributions of land-cover types and climate zones, which provide necessary background information for interpreting the spatial heterogeneity of the VPD–E relationship.

Manuscript Number: EGUSPHERE-2025-4820

Revision notes: We sincerely thank the Editor for recognizing the potential value of our manuscript and for the time and effort devoted to handling the review process, including coordinating the assessment and inviting expert reviewers. We also sincerely thank both reviewers for their careful reading of our manuscript and for providing constructive and insightful comments. These comments have helped us identify several aspects that require further clarification and improvement. In the following response, the reviewers’ comments are shown in black, and our point-by-point responses and planned revisions are shown in blue. At this final-response stage, we provide a detailed explanation of how we will revise the manuscript according to the Editor’s and reviewers’ suggestions. We will carefully improve the Methods, Results, Discussion, figures, captions, and uncertainty analysis to make the manuscript more rigorous, transparent, and internally consistent. We respectfully hope to be given the opportunity to revise the manuscript and incorporate the suggested modifications.

Responses to the reviewer’s comments

Response to Reviewer #2

We sincerely thank the reviewer for the careful evaluation of our manuscript and for providing detailed and constructive comments. We fully recognize that the current version still requires substantial clarification in several key aspects, including the distinction between established knowledge and our quantitative contribution, the assumptions and spatial scale of the moving-window approach, the consistency between temporal trends and sensitivity estimates, the readability of Figures 1 and 2, the reporting of nonlinear VPD thresholds, and the interpretation of the SEM pathways. In the responses below, we provide a point-by-point plan for how we will revise the manuscript in accordance with the reviewer’s suggestions. We will carefully check the data processing, clarify the temporal and spatial definitions of the analyses, correct inconsistent threshold values, strengthen the uncertainty discussion, and revise the figures, captions, Methods, Results, and Discussion accordingly. We sincerely appreciate the reviewer’s comments, which have helped us identify several aspects that require further clarification and improvement, and we respectfully hope to be given the opportunity to revise the manuscript and incorporate these suggested modifications in a more rigorous and transparent way.

Comment 1: The abstract should better distinguish established knowledge from the novel contribution of this study. Statements that VPD affects evapotranspiration, that soil moisture limitation constrains E in arid regions, and that transpiration is an important pathway are broadly consistent with existing understanding. The main potential contribution of this work is the global quantification of VPD–E sensitivity, regional heterogeneity, and apparent VPD thresholds. The abstract should emphasize this quantitative contribution rather than presenting well-known mechanisms as new findings.

Response: Thank you for this constructive comment. We agree that the abstract should more clearly distinguish established knowledge from the new contribution of this study. The effects of VPD on evapotranspiration, the role of soil moisture limitation in dry regions, and the importance of transpiration are already well recognized in previous studies. In the revision, we will avoid presenting these mechanisms as novel findings. Instead, we will use them only as background to introduce the scientific context.

We will revise the abstract to emphasize the quantitative contribution of this study. Specifically, we will highlight that our work provides a global estimate of VPD–E sensitivity, identifies the spatial heterogeneity of this sensitivity across climate and land-cover conditions, and quantifies apparent VPD thresholds in different aridity zones. We will also report the main numerical results more directly, including the global mean sensitivity, the proportion of land showing positive or negative VPD–E sensitivity, and the aridity-zone-dependent threshold values. These changes will make the abstract more focused on what this study adds beyond existing knowledge, rather than restating mechanisms that are already broadly understood.

Comment 2: The space-for-time moving-window approach is not sufficiently explained. This method can be useful for estimating the local apparent sensitivity of E to VPD, but it relies on strong assumptions: neighbouring pixels must share similar climate, soil, vegetation, topography, and management conditions, and spatial VPD gradients must be a reasonable proxy for temporal VPD changes. The manuscript should clarify how “similar background climate” was defined and whether differences in soil moisture, radiation, elevation, vegetation structure, and land management were controlled.

Response: We agree that the space-for-time moving-window approach needs a clearer explanation. This method is useful for estimating the local apparent sensitivity of E to VPD, but it does rely on important assumptions. In the revision, we will explicitly define this estimate as an apparent local sensitivity, rather than a strict causal response or a direct substitute for temporal VPD change. We will also clarify that “similar background climate” is mainly constrained by the local moving-window design. All calculations will be conducted on a common 0.1° grid, and each target pixel will be compared only with candidate pixels within a 5 × 5 grid-cell neighborhood, corresponding to a 0.5° × 0.5° local window. This local window is intended to reduce large-scale climatic differences among compared pixels.

We will further clarify how surface heterogeneity is controlled within this window. Candidate pixels will be retained only when they share the same dominant MODIS land-cover type as the target pixel. The difference in the fraction of the dominant land-cover type will be required to be less than 10%. Pixels with an elevation difference greater than 100 m from the target pixel will also be excluded. These screening rules are designed to reduce differences in vegetation type, land-cover composition, and topographic conditions. We will describe these criteria more explicitly in the Methods section, so that readers can understand how the moving-window samples are selected.

At the same time, we will state the limitations of this approach more clearly. Soil moisture, radiation, soil properties, vegetation structure, and land management cannot be fully controlled by the current moving-window screening. They may still differ among neighboring pixels, especially in irrigated croplands, heterogeneous landscapes, mountainous areas, and coastal regions. We will therefore avoid claiming that the method fully isolates the independent effect of VPD. Instead, we will explain that the method reduces some major sources of spatial heterogeneity but cannot remove all confounding factors. We will add this limitation to the Methods and Uncertainty sections, and we will interpret dE/dVPD as a spatially derived apparent sensitivity that should be considered together with the partial correlation, SEM, and independent dataset checks.

Comment 3: The spatial scale of the moving window needs to be explained. The authors mention using a 5 × 5 km window, but this seems quite small relative to the spatial resolution of the input datasets. The manuscript should clarify the final resolution of the analysis and the resampling procedure for each dataset. Sensitivity tests with different window sizes are recommended to demonstrate the robustness of the results.

Response: Thank you for this helpful comment. We agree that the original description of the moving-window scale was confusing. The analysis is not conducted using a 5 × 5 km window. In the revision, we will correct this wording and clearly state that all analyses are conducted on a common 0.1° grid. This grid is used for the ERA5-Land VPD data and the GLEAM E data in the main analysis. The moving window will therefore be defined as a 5 × 5 grid-cell neighborhood, corresponding to a 0.5° × 0.5° window centered on each target pixel, rather than a 5 × 5 km window.

We will also add a clearer description of the resampling and harmonization procedure. Continuous variables will be harmonized to the common 0.1° grid before the moving-window calculation. For the MODIS MCD12C1 land-cover product, we will use the dominant land-cover type and the fractional cover of that dominant type within each 0.1° grid cell. Candidate pixels will be retained only when they share the same dominant land-cover type as the target pixel, and when the difference in dominant land-cover fraction is less than 10%. Elevation will also be harmonized to the same grid, and candidate pixels with an elevation difference greater than 100 m from the target pixel will be excluded. We will describe these rules explicitly in Section 2.2.2 so that the spatial scale, input resolution, and sample selection procedure are clear.

Following the reviewer’s recommendation, we will add a window-size sensitivity test to evaluate the robustness of the moving-window result. In addition to the main 5 × 5 grid-cell window, we will test alternative window sizes, such as 3 × 3 and 7 × 7 grid-cell neighborhoods, while applying the same land-cover and elevation screening criteria. We will compare the resulting dE/dVPD patterns with the main result in terms of sign consistency and broad spatial distribution. These results will be reported in the Supplementary Information. This additional test will help show whether the main spatial pattern of VPD sensitivity depends strongly on the selected window size. We will also soften the interpretation of dE/dVPD and describe it as a local apparent sensitivity, rather than a fully isolated causal effect of VPD.

Comment 4: The abstract states that E increases with rising VPD over 60.7% of the land surface, whereas Section 3.1 states that E decreases in areas where VPD increases, suggesting that atmospheric drying has an inhibitory effect. Section 3.2 then reports that 60.71% of land areas show positive sensitivity. This description is confusing and should be clarified.

Response: We agree that the original wording may cause confusion because it mixed two different quantities: the temporal trend of E and the sensitivity of E to VPD. The statement in Section 3.1 refers to the spatial co-occurrence of long-term E trends and VPD trends, whereas the result in Section 3.2 refers to the estimated local apparent sensitivity, dE/dVPD. These two results are related but not identical. Therefore, the manuscript should not describe them using the same wording.

In the revision, we will clarify this distinction throughout the Abstract, Section 3.1, and Section 3.2. In the Abstract, we will avoid the ambiguous phrase “E increased with rising VPD.” Instead, we will state that 60.7% of the global land surface shows positive apparent sensitivity of E to VPD, while negative or weak sensitivity is mainly found in water-limited regions. In Section 3.1, we will describe only the long-term temporal changes of VPD and E. In Section 3.2, we will separately present the spatial dE/dVPD sensitivity results. This revision will make clear that the 60.7% value refers to sensitivity, not to the fraction of land where E and VPD both increased over time. We will also adjust the related wording to avoid implying that rising VPD universally enhances E, because the response depends strongly on soil moisture and aridity conditions.

Comment 5: Figure 1 is difficult to interpret. The colorbars in panels (a) and (b) should clearly indicate the variable, unit, and sign convention. The text states that significant increases in VPD cover 76.22% of the land surface, but the figure does not clearly distinguish significant from non-significant trends. Hatching, stippling, or a separate significance panel should be added. I also suggest revising the colorbars so that decreases and increases are clearly represented by blue and red, respectively, with zero set as the midpoint.

Response: The colorbars and significance information should allow readers to understand the variable, unit, sign convention, and statistical significance directly from the figure. In the revision, we will redraw Figure 1 to improve its readability. For the VPD trend map, the colorbar will be labelled as Sen’s slope of VPD, with the unit hPa yr^-1. For the E trend map, the colorbar will be labelled as Sen’s slope of E, with the unit mm yr^-1. We will also explicitly define the sign convention: positive values indicate increasing trends, while negative values indicate decreasing trends.

Following the reviewer’s suggestion, we will use a diverging color scale for the trend maps, with zero set as the midpoint. Decreasing trends will be shown in blue and increasing trends in red. This will make the spatial pattern of increases and decreases easier to interpret. We will also add stippling to indicate pixels with statistically significant trends at P < 0.05. This change will make it clear which regions contribute to the reported significant increase in VPD over 76.22% of the land surface, and which regions show non-significant trends.

We will further revise the figure caption to describe these changes explicitly. The revised caption will state that panel a shows Sen’s slope of VPD, with stippling indicating significant trends at P < 0.05, and that panel d shows Sen’s slope of terrestrial E with the same significance convention. We will also clarify that the insets show the percentage distribution of pixels among slope classes. For the time-series panels, we will clearly identify the monthly global mean, annual mean, 5-year moving average, anomalies relative to 1981–2020, and Theil-Sen or piecewise trends where applicable. These revisions will make Figure 1 more self-contained and will better connect the visual evidence with the quantitative statements in the text.

Comment 6: Figure 2b is particularly confusing. A large fraction of samples appears to fall between approximately -7 and -28 °C. While such temperatures may occur in high-latitude or high-elevation regions, their dominance raises concerns that cold, ice-covered, tundra, or very low-E regions may be overrepresented. The precipitation range also appears too low if it represents annual precipitation; humid and tropical regions commonly exceed 300 mm yr⁻¹. The authors should clarify whether the axes show annual means or monthly values, whether each point represents a grid cell or a grid cell–month sample, whether area weighting was applied, and whether permanent ice/snow, barren, or very low-LAI regions were excluded.

Response: We agree that Figure 2b is difficult to interpret in its current form and may raise concerns about sample representation and data consistency. In the revision, we will recheck the input data used for this panel and redraw the figure after confirming the temporal scale, spatial sampling unit, and masking procedure. In particular, we will verify whether the temperature and precipitation values represent annual means, monthly means, or grid cell–month samples. We will also check the precipitation calculation to ensure that the unit and aggregation period are consistent with the interpretation in the figure and text. We will revise the figure caption and Methods to state explicitly what each point represents. We will clarify whether the points are individual grid cells or grid cell–month samples, and whether area weighting is applied when summarizing the distribution. If the figure is based on monthly samples, we will label the axes accordingly and avoid interpreting the values as annual climatic conditions. If annual means are used, we will ensure that precipitation is aggregated as annual total precipitation rather than monthly precipitation, and we will revise the axis unit accordingly.

We also agree that cold, ice-covered, tundra, barren, or very low-vegetation regions may be overrepresented and may distort the relationship shown in Figure 2b. In the revision, we will apply additional screening before remaking this figure. Permanent snow and ice areas will be excluded. We will also examine whether barren land, tundra, and very low-LAI pixels should be removed or shown separately, because these regions often have very low E and may not represent the vegetation-mediated VPD-E response that is central to this study. We will clearly report the final masking criteria in the Methods and figure caption.

After these checks, we will remake Figure 2b and revise the related text accordingly. If the corrected figure still shows a large influence of cold or low-E regions, we will interpret this pattern explicitly rather than treating it as a general global response. If the original pattern is found to result from inconsistent temporal aggregation or insufficient masking, we will replace it with the corrected result. These revisions will make the figure more transparent and will ensure that the conclusions drawn from Figure 2b are supported by reliable and consistently processed data.

Comment 7: Some explanations in the Results section, particularly lines 284–287, are not fully supported. In addition, the threshold results are unclear. Line 333 mentions that the arid-zone threshold is about 1.31 kPa, lines 337–338 refer to “temperate zones” with values of 1.67–1.68 kPa, and lines 406–409 again mention “arid regions” with values of 1.67–1.68 kPa. These inconsistencies should be corrected.

Response: In particular, the text around Lines 284-287 mixed direct result description with mechanistic interpretation. In the revision, we will rewrite this part to ensure that the Results section only reports patterns that are directly supported by the data. Interpretations related to physiological regulation, soil moisture limitation, and land–atmosphere feedbacks will be moved to the Discussion, where they can be linked more explicitly to the partial correlation analysis, SEM path coefficients, and previous literature. We will also avoid causal wording unless it is directly supported by the corresponding analysis.

We also agree that the threshold values were inconsistent in the previous version. The manuscript originally mixed different classification schemes and therefore reported incompatible values, including an arid-zone threshold of about 1.31 kPa, “temperate-zone” thresholds of 1.67-1.68 kPa, and later references to arid regions using the same 1.67-1.68 kPa range. This inconsistency will be corrected throughout the manuscript. In the revision, we will use one consistent classification framework based on aridity zones. The threshold results will be reported as: arid zone, semi-arid zone, semi-humid zone and humid zone. The earlier inconsistent values will be removed from the Results, Discussion, Abstract, and Conclusion.

We will further clarify how these thresholds are derived and interpreted. We will state that the thresholds are estimated from the comparison of one-breakpoint piecewise regression, two-breakpoint piecewise regression, and GAM, rather than being assumed a priori. The reported values will be described as dominant transition points in the observed nonlinear VPD-E relationship, not as exact or universal ecological limits. We will also revise the wording in Section 4.1 and the Discussion to avoid switching between “climate zones,” “temperate zones,” and “arid regions” without a consistent definition. These changes will make the threshold analysis internally consistent and will prevent overinterpretation of the nonlinear response results.

Comment 8: Piecewise regression and GAM are major components of the study, yet their results are mainly presented and interpreted in the Discussion. The Results section should include a clear subsection on nonlinear VPD thresholds, including model fits, threshold estimates, uncertainty ranges, and differences between methods. The Discussion should then focus on the ecohydrological interpretation of these thresholds.

Response: We agree that the nonlinear threshold analysis is an important part of the study and should be presented more clearly in the Results section. In the previous version, the piecewise regression and GAM results were described mainly in the Discussion, which made it difficult for readers to distinguish the statistical results from their ecohydrological interpretation. In the revision, we will restructure this part of the manuscript.

Specifically, we will add a new subsection in the Results section focusing on nonlinear VPD thresholds. This subsection will report the fitted VPD–E relationships for different aridity zones, the threshold estimates, and the comparison among one-breakpoint piecewise regression, two-breakpoint piecewise regression, and GAM. We will also include model performance information, such as R², RMSE, AICc or ΔAICc where applicable, to show why the selected threshold model is used as the main result. The revised Results section will clearly report the threshold estimates for each aridity zone, including the arid, semi-arid, semi-humid, and humid zones. We will also add uncertainty information for the threshold estimates, such as confidence intervals or ranges derived from model fitting and sensitivity tests, depending on the final model output. We will then revise the Discussion so that it no longer serves as the main place for presenting the threshold results. Instead, the Discussion will focus on the ecohydrological meaning of these thresholds. For example, we will discuss why higher thresholds occur in arid and semi-arid zones, how water limitation may weaken the positive response of E to VPD under high atmospheric demand, and why humid regions show lower transition values. We will also emphasize that these thresholds should be interpreted as dominant transition points in the observed nonlinear VPD-E relationship, rather than exact or universal ecological limits. This revision will make the manuscript structure clearer and will separate statistical reporting from process-based interpretation.

Comment 9: The SEM analysis needs clearer reasoning and explanation. As I understand it, the model suggests that changes in VPD can explain 66% of the variation in LAI. I am skeptical that this reflects a direct effect of VPD on LAI; it may instead reflect spatial climate gradients or seasonal changes. High VPD does not necessarily promote leaf area development; it usually imposes atmospheric water stress. This pathway may simply indicate that warmer, more productive regions have both higher VPD and higher LAI. The authors should test the SEM using deseasonalized anomalies, climate-zone-specific models, and additional control variables such as radiation, soil moisture, vegetation type, and growing-season length. Also, what do the red lines in Figure 5 represent?

Response: We agree that a strong VPD-LAI pathway should not be interpreted simply as evidence that high VPD directly promotes leaf area development. High VPD usually represents atmospheric water stress and can suppress stomatal conductance and vegetation activity, especially under limited soil water availability. Therefore, the positive VPD-LAI pathway in the current SEM may partly reflect shared spatial climate gradients, seasonal covariation, or the fact that warmer and more productive regions tend to have both higher VPD and higher LAI. In the revision, we will avoid interpreting this pathway as a direct physiological effect of VPD on LAI.

To address this concern, we will revise both the SEM analysis and its explanation. First, we will test the SEM using deseasonalized anomalies to reduce the influence of seasonal cycles. Second, we will conduct climate-zone-specific SEM analyses to examine whether the VPD-LAI relationship is consistent across arid, semi-arid, semi-humid, and humid regions, rather than being driven by broad spatial climate gradients. Third, we will include additional control variables where data availability allows, including radiation, soil moisture, vegetation type, and growing-season length. These tests will help determine whether the VPD-LAI pathway remains robust after accounting for climatic seasonality, vegetation background, and water-energy constraints.

We will also strengthen the mechanistic discussion by reviewing relevant studies on the effects of VPD on vegetation growth, canopy development, stomatal regulation, and plant water stress. If the revised SEM shows that the VPD-LAI pathway weakens substantially after deseasonalization or within climate zones, we will explicitly state that the original pathway was likely influenced by covariation rather than a direct VPD effect. If the pathway remains significant in some regions, we will still interpret it cautiously as a statistical association under specific environmental conditions, not as general evidence that high VPD promotes vegetation growth. In addition, we will revise Figure 5 and its caption to make the path notation clear. We will explicitly explain what the red lines represent, including whether they indicate positive standardized path coefficients, negative path coefficients, or statistically significant pathways, depending on the final graphical convention. We will also add the standardized coefficients, significance levels, and model-fit statistics more clearly in the figure or caption. These revisions will make the SEM more transparent and will prevent readers from overinterpreting the VPD-LAI pathway as a direct causal mechanism.

Additional planned revisions

In addition to the revisions made in response to the two reviewers’ comments, we also plan to make several additional changes to improve the rigor and clarity of the manuscript. First, we will remove the analysis of canopy interception from Section 3.3. The reason is that the mechanisms linking VPD to transpiration and soil evaporation are relatively well established, whereas the influence of VPD on canopy interception is less direct and more difficult to interpret mechanistically. Canopy interception is mainly controlled by canopy structure, commonly represented by LAI, and by precipitation amount and intensity. Therefore, a statistical relationship between VPD and interception may not necessarily indicate a causal mechanism. To avoid overinterpretation, we will focus Section 3.3 on transpiration and soil evaporation, for which the physical and ecohydrological interpretations are more robust. Second, we will detrend the data used in Section 3.3 before conducting the correlation analysis, so that the results better reflect interannual covariation rather than shared long-term trends. We will also add significance markers in Figure 3 to indicate correlations significant at P < 0.05. Third, we will add the official citation for ERA5-Land to improve the completeness of the data-source description. Finally, we will add supplementary maps showing the spatial distributions of land-cover types and climate zones, which provide necessary background information for interpreting the spatial heterogeneity of the VPD–E relationship.