Sap flow as a function of variables within nested scales: ordinary least squares vs. spatial regression models

Understanding scale-dependent influential drivers of sap flow variability can help managers and policymakers to allocate resources within a particular scale to improve forest health and resiliency against water-stress stimuli such as drought and insects, e.g. bark beetle infestations. We defined a daily measure of sap flow as a function of variables within nested scales of landscape, stand, and tree, using ordinary least squares (OLS), spatial lag and error regression models. Model covariates were elevation, latitude (Y-coordinate), longitude (X-coordinate), neighborhood tree density, tree diameter at breast height, and bark temperature for 40-surveyed Norway spruce (Picea abies) in the Czech Republic, Central Europe. Trees were spatially distributed within 19-established subplots across five plots, with distances ranging 2–9 km, at which variations in soil water potential and temperature were limited. The daily measure of sap flow within the regional scale allowed us to avoid the temporal and spatial variability of climate effects on sap flow. A relatively flat terrain across subplots also allowed us to control the effects of slope, aspect, and topography-related solar incidence angle on sap flow. Sap flow was strongly spatially autocorrelated, so OLS models failed to take spatial autocorrelation into account unless to some extent, depending on the spatial distribution of samples, by including latitude and/or longitude in the models. Among spatial regression models, spatial error models performed better than lag models, allowing to capture the effects of unmeasured independent variables. Sap flow variability for the most part (∼70%) was explained by the landscape-level variable of elevation followed by the stand-level variable of tree density, and the remaining part by variables related to tree characteristics; a nested down-scaling function, defined and visualized for the first time. Therefore, thinning forest stands and future plantations with optimum distances, based on the elevation gradients, may be required to counterbalance the allocation of resources, e.g. water, nutrients, and light, among trees, leading to enhance forest health and resiliency against water-stress stimuli.


Introduction
Forest is a major terrestrial ecosystem that includes 30% of the land area on the earth (Schmitt et al 2009). Forest ecosystems play a critical role in carbon sequestration and rainfall infiltration that secures drinking water supplies on a global scale (Uy and Shaw 2012). Forest health and responses to water-stress stimuli can be assessed by monitoring sap flow to understand tree water-use and transpiration functions (Alarcón et al 2000, Özçelik et al 2022, West et al 2023.
Most studies to understand factors that influence sap flow variabilities either focused on variables acting at the level of an individual tree or a subset of environmental variables, with few on combined sets of these variables. Some of the tree characteristics found to be significant were angles of tree deviations from the axis, stem moisture, temperature inside the trunk, and the normalized difference vegetation index (NDVI; Safargalieva et al 2021), highly correlated with leaf biomass. Some of the environmental variables found to be significant include soil water potential (Ježík et al 2015, Cao et al 2020, Korakaki and Fotelli 2021, solar radiation and vapor-pressure deficit (VPD; Ježík et al 2015, Korakaki andFotelli 2021, Safargalieva et al 2021), photosynthetically active radiation (PAR; Cao et al 2020, Suárez et al 2021, air temperature (Cao et al 2020), surface terrains causing variations in slope (Kume et al 2016, Schoppach et al 2021 and aspect (Hassler et al 2018), and soil texture (Hassler et al 2018). Tree density was found to be another important driver of sap flow variability (Hassler et al 2018, Suárez et al 2021. Parametric models, such as regression models, were commonly used to define the relationships between sap flow and its driving factors even though the number of developed models were found to be few. For example, Hassler et al (2018) used multiple linear regression models to study the relative importance of tree-, stand-, and site-specific characteristics on the sap velocity of 61 beech and oak species across 24 sites in Luxembourg. A combination of the tree-and site-specific factors such as tree species and diameter, and site geology, e.g. soil texture and aspect, was found to be the main drivers of sap velocity across the landscape (Hassler et al 2018). In addition, Safargalieva et al (2021) used auto-regressive moving-average and linear regression models to predict sap flux density and to study the relationship between sap flux density and treeand environmental factors, respectively. A combination of air temperature and VPD, trunk deviations from the axis, temperature inside the trunk, stem humidity, and NDVI was found to be the best variables to explain the variability in sap flux density (Safargalieva et al 2021). Suárez et al (2021) also used a mixed-effects linear regression model to define the relationship between sap flow and several microclimatic attributes of different agroforestry systems in the western Colombian Amazon. The authors found PAR, and thus incident radiation, as the most important factor affecting sap flow among other investigated factors such as relative air humidity, air temperature, and VPD.
Given all these studies, information about the definition and importance of different scales and how factors, if spatially autocorrelated, at different levels of scales may influence sap flow has yet to be fully investigated. For example, as most environmental variables are spatially autocorrelated and may non-linearly influence sap flow, selecting a particular type of regression model (non-linear, OLS, or spatial regression models) is still a question of concern. When the variables of interest, e.g. sap flow, are spatially dependent, but the independency is assumed in a regression model, the type I error (false rejection of a true null hypothesis) will be inflated (Lennon 2000). We proposed to (1) define the scale-dependent factors, changing at the tree-, stand-, or landscape-level, that can explain the most sap flow variability in Norway spruce-dominated forests; (2) examine if latitude and longitude can be substitutes for space effects in the OLS models; (3) take the existing spatial autocorrelation into account by developing different types of spatial regression models; (4) investigate whether stand-and landscape-level driving factors of sap flow can compensate for the effects of factors acting at the level of an individual tree, with further management implications; and (5) providing foundations for future research, in any discipline, to choose an appropriate regression model.

Study area and field measurements
The study focused on five established plots within the Norway spruce-dominated forest at the university forest enterprise in Kostelec nadČernými Lesy, Czech Republic (Central Europe; figure 1).
Each plot consisted of four subplots, i.e. A, B, C, and D, at which 8-10 Norway spruce trees were selected with an average DBH of 40 cm (figure 1). Two plots (no. 3 and 4) were established within proximity at 430 m while distances among other plots and from these two plots were relatively far, ranged 2-9 km (figure 1(A)). Distances among subplots within a plot ranged 40-250 m. The six established plots, within a 100 m elevation range, were designed under EXTEMIT-K project. The project was proposed to conduct research experiments to deliver science-based approaches for dealing with the environmental challenges of protecting Norway spruce forests in the Czech Republic (Trubin et al 2023). Bark beetle infestations and periodic drought effects on trees and beetle infestation dynamics were the main challenges in this project (Trubin et al 2023).
The topography of the area is characterized by a flatter landscape with slope ranging from 0.8 • to 6 • . Sap flow (kg h −1 ) was continuously measured, at 10 min intervals, using EMS81 sap flow sensors (EMS Brno, Czech Republic (CZ)) installed on the North (N) side of each stem at ca. 2.5 m height. The sensors use coupled electrodes and thermocouples to measure the vertical water flux rate per cm of stem circumference using the trunk heat balance method (Čermák et al 2004).
The bark surface temperature ( • C) was monitored using two IRT sensors (Apogee, Logan, UT, USA) installed 3 m above ground on the North and South (S) sides of the tree trunks equipped with sap flow   1(B)). We additionally recorded the latitude, longitude, and elevations (m) at which trees were located. We compiled 12 h (6 am-6 pm) measurements of sap flow on sum and bark temperature (N and S) on average for 40 trees on a cloudless day, 19 April 2019, among several cloudy days before and after. The air temperature was 17 • C on average. A single-day measure of sap flow at the beginning of the growing season (May-October), when trees were not under water stress, allowed us to avoid drought effects and smoothing data variabilities if averaged for longer temporal scales. However, we additionally tested the correlation between the sap flow measures for this day and the 6 month growing season on average and found a very high association between them (r = 0.89). Soil water potential (Ψw) and temperature ( o C) were also measured for each subplot; however, we did not include them in the models as they were a single measurement for each subplot with very limited variabilities across subplots. None of the trees was found to have deviated from the axis. In addition to field-measured data, we counted the number of neighboring trees within a 10 m-radius buffer (314 m 2 ) from each tree using 20 cm UAV imagery, underlying tree positions, in ArcGIS Desktop (ver. 10.8.1, ESRI). This variable was considered as neighborhood tree density (hereafter tree density).
We further defined scales for each variable (table 1). For example, the scale on which we measured a variable based upon was considered a measurement scale. Significant changes in some of the variables measured at those scales happen within the extent (scale) to which those changes/variations can have Figure 2. Box plots of two trees, 5A1 and 5B7, were removed from the analyses. The sap flow sum (kg/12 h) for the day of analyses, 19 April 2019, was found to be outliers (black-filled circles) for these two trees.
considerable influences on most biological factors including sap flow. Thus, we defined that broader extent or scale as a biologically-defined scale (table 1).
The spatial and temporal scales of this study also allowed us to keep the effects of climate variability on hold in the proposed models while investigating the daily influential factors of sap flow at the regional scale (figure 1).

Data screening and consistency
We initially used datasets of all six-established plots (1A-D, 2A-D, 3ABD, 4A-D, 5AB, and 6A-D; figure 1(A)) and screened them for any missing data for bark temperature and sap flow for the day of analyses, 19 April 2019; trees with missing data were removed from the pool. We then used box plots to check for data consistency for the day of analyses, using 13 d of data; 6 d-before and -after 19 April 2019. If bark temperature or sap flow for the day of analyses were found to be outliers, the associated trees were removed from the dataset. The sap flow sum for two trees (5A1 and 5B7; figure 2) was found to be outliers; thus, those trees were removed from the dataset. After removing the two trees with sap flow outliers, we removed plot 5 from the analyses as the remaining three trees in plot 5 were found to be an unbalanced sample size compared to other plots with 7-9 trees.
One tree per subplot was randomly selected to display and compare the temporal patterns of sap flow sum for a day before and after the day of analyses, and a week later (figure 3). Twenty out of 40 trees (half of the sample size) were randomly selected for this purpose.
The temporal patterns of sap flow sum for the three consecutive days (18-20 April) were found to be very similar. The patterns showed more variations by moving a week from 18 April 2019 (figure 3). The observed temporal autocorrelation patterns of sap flow were likely to be driven by temporal autocorrelation of air temperature and subsequent water-pressure deficit. The air temperature on average was 15 • C, 17 • C, 18 • C, and 21 • C for 18-20 and 25 April, respectively.

Plot-based statistics
Plot-based statistics, such as mean and standard deviation, were provided for the response (sap flow sum) and explanatory variables. However, we were not able to provide this information for subplots within plots as the number of trees within subplots was found to be few. In addition, the non-parametric Kruskall-Wallis statistical test (Kruskal and Wallis 1952) was conducted to examine if the mean of variables significantly varies across plots at a 95% confidence interval (p-value = 0.05). April. 20 out of 40 trees (half of the sample size; one tree and four trees per subplot and plot, respectively) were randomly selected. Two trees (3A2 and 3A4) in subplot 3A were selected to compensate for a missing subplot (3C) from the original datasets. (1)), in the spatial lag model (equation (2)), a coefficient (ρ) for the spatial weight matrix of the dependent/response variable (Wy) is added to control the spatial dependency, as confounding effects, on the significance of all independent variables under considerations (β 1 (X 1 ) + β 2 (X 2 ) + . . . + β n (X n )).

Spatial regression models: lag vs. error model Different than OLS model (equation
In the spatial error model (equation (3)), the error (ε) is divided as independent identically distributed errors (u) and spatially autocorrelated error (spatially weighted error; Wε in equation (4)). λ is a coefficient for the spatially autocorrelated error term (equation (4)).

Regression diagnostics to choose an appropriate regression model
Following the prepared hierarchical steps (figure 4), we initially ran Pearson's correlation test among all variables in R Core Team (2021) to use one out of two highly correlated explanatory variables (r ⩾ 0.7) to avoid data redundancy in the proposed models. As elevation was highly correlated with longitude (r = −0.78; table 2), we avoided including both in any of the proposed models.
We developed a spatial weight matrix using the first-order Queen Contiguity weight algorithm (GeoDa; Anselin et al 2006) on a vector shape file of tree positions (ESRI 2021). The weight matrix was used in the OLS models to check the assumptions of normal distribution and spatial autocorrelation in model residuals. If the assumption of normality was met, we proceeded to the next step of checking the availability of spatial autocorrelation in model residuals. When there was no spatial autocorrelation in model residuals, we kept the OLS models as the best-fitted models; otherwise, we ran spatial lag and error models to take the spatial autocorrelation into account (figure 4). We kept the most parsimonious models, with fewer predictors while   resulting in a desired level of goodness of fit, e.g. R 2 ⩾ 50, and evaluated them based on Akaike Information Criterion (AIC; Akaike 1973). Different combinations of the explanatory variables were also proposed to address different objectives of this research. Models with an AIC-difference (∆AIC) less than and equal to two were considered equal-performing models. The Jarque-Bera, Moran's I, and Breusch-Pagan tests were used to assess normality, spatial autocorrelation, and heteroskedasticity of error, respectively. All the steps taken for regression diagnostics and modeling were done in the GeoDa open-source software (Anselin et al 2006). We additionally examined spatial autocorrelation (Univariate Moran's I) in sap flow measures using the developed spatial weight matrix.

Plot-based variability in sap flow and explanatory variables
The variations in sap flow and the proposed explanatory variables for sap flow were statistically significant across five plots; all p-value for the Kruskal-Wallis statistical test was significant at a 0.05 error rate (table 3). The highest sap flow on average was in plot 2 with the associated highest elevation and lowest tree density on average (table 3). The bark temperature was also lowest for plot 2 (table 3).

OLS vs. spatial regression models
Sap flow was strongly spatially autocorrelated (Moran's I = 0.50, pseudo p-value with 999 permutations = 0.001). All assumptions for model residuals including normality, heteroskedasticity, and multicollinearity were met. The assumption of spatial dependency of model residuals was violated in all OLS models except in some of them (no. 13-17) with having latitude and longitude included (significant probability for Moran's I; table 4). When latitude and longitude were not included, the spatial regression models were found to be appropriate models to take the spatial autocorrelation into account (models no. 1-8; table 4). In some of the OLS models (no. 9-12), with fewer variables included than other models, latitude and/or longitude could not provide substitutes for the space effects (significant probability for Moran's I; table 4). Among the spatial regression models, spatial error models were found to be better performing than spatial lag models (lower AIC; table 4). When variables such as DBH and bark temperature were added to the top-ranked spatial regression models (no. 1 vs. 2-4), the results were not significantly changed (∆AIC ⩽ 2; table 4). As longitude was highly correlated with the elevation, we did not include these two variables together in any of the proposed models to avoid the multicollinearity issue (table 4). A summary of how much each explanatory variable was competent to explain the sap flow variability is also provided at which elevation and DBH were found to be the highest-and lowest-ranked significant variable, respectively (table 5). Tree density was found to be the second most important variable (table 5). None of the variables in the models violated a measure of multicollinearity, the variance inflation factor (VIF < 3; table 5).
All regression coefficients for the top-ranked spatial regression model (no. 1; table 4) were also standardized using equation (5).

Standardized Regression Coefficient (SRC
where β i , σ (X i ) , and σ (Y) stand for the unstandardized regression coefficient, the standard error for an independent (explanatory) variable, and the standard error for the dependent (response) variable, respectively. The statistics, e.g. coefficient and significance (p-value) of parameters, for the top-ranked spatial regression model (no. 1) are provided for further applications, for example, to adjust field-measured sap flow based on the standardized coefficients of the associated variables and Y-intercept (table 6).

Sap flow as a function of variables within nested scales
Results from plot-based statistics and regression models (no. 1 and 17) indicate that sap flow can be defined as a function of variables, such as elevation, tree density, and individual tree traits such as bark temperature, acting within three scales that are nesting in each other. In models no. 1 and 17, elevation, tree density (significant at a 0.1 error rate for model 17), latitude, and bark temperature were significant factors explaining sap flow variability (table 4). The variable significance (%; table 5) for model no. 17 is also indicating the order of importance for landscape-, stand-, and tree-level factors.
Based on these results, we provided two different ways of visualizing the main drivers of sap flow variability within nested scales from the landscape-to stand-and tree-level, for future reference (figures 5 and 6). We also introduced a new legend to show how different variable gradients may act within nested scales (figure 5). Table 4. Proposed OLS and spatial, lag and error, regression models with associated Moran's I p-value and AIC to assess spatial dependency in model residuals and model performance, respectively. X, Y, NoTrees, and BarkTemp are designated names for longitude, latitude, tree density, and bark temperature on the North and South sides of trunks on average, respectively. Variables with a star ( * ) were not statistically significant (α = 0.05) in OLS, spatial lag, or error regression models with the associated digits of 1, 2, and 3, respectively. For example, in model no. 2, bark temperature was not significant in any of the models. Models in bold (no. 1 and 13) are top-ranked spatial regression and OLS models with the least AIC among other competing models in their categories. In models no. 13, 15, and 17, the variable in bold (NoTrees) was not significant at a 0.05 error rate, but it was significant at a 0.1 error rate (90% confidence interval). A group of models in rows were divided to provide the required information to address the research objectives (purpose of modeling). The spatial regression models were run after if the OLS model residuals were spatially dependent, assessed by the significant probability of Moran's I.
Purpose of modeling    Figure 5. A schematic representation of sap flow as a function of variables acting within three nested scales of landscape-, stand-, and tree-level. The study area follows the 100 m elevation gradients as landscape-level variations. The five established plots showed gradients of tree density, as stand-level variations, which were negatively correlated with the elevation. Individual tree traits such as bark temperature and DBH, as tree-level variations, were positively and negatively correlated with the tree density, respectively. A new graphical legend is introduced to convey the concepts of nested scales for variable gradients.

Sap flow as a function of variables within nested scales
We were able to explain 2/3 of sap flow variability at which elevation and neighborhood tree density, defined as landscape-and stand-level factors, respectively, were found to be the most influential variables on daily sap flow variability (table 4; spatial lag and error model no. 1). Factors measured at the tree-level such as Figure 6. A schematic representation of influential factors on sap flow variability at nested scales from the landscape-to standand tree-level. The size and direction of the blue arrows represent the magnitude and direction of factors acting or interacting at different scales. At the landscape level, when elevation increases, the decrease in air temperature and vapor pressure influences sap flow. Higher elevations are generally associated with the steeper landscape that causes variations in incidence angle. When the incidence angle decreases (either due to landscape terrains, as shown here, or sun positions at different times of the day, season, or year), the intensity of solar irradiance increases. Within the stand level, the density and spatial configurations of trees determine the resource (water, nutrients, and solar irradiance) allocations among trees and heat loads on trees (influencing bark temperature) and interspace forest grounds that all together influence sap flow. At the tree level, variations in DBH and aboveand below-ground tree structures are other potential driving factors of sap flow variabilities.
DBH did not add much to explaining sap flow variabilities, similar to Hassler et al (2018). However, bark temperature was found to be significant in model with having longitude or latitude included (models no. 13, 14, 16, and 17; table 4). The low variability of DBH and bark temperature across our plots may have caused these variables to be insignificant in other models even though both variables were found to be correlated with sap flow (table 2). The remaining part of sap flow variability may be due to some other tree-level attributes we did not fully account for. For example, variations in needle mass (Lagergren and Lindroth 2004) and underground structure of tree roots (Schoppach et al 2021) were found to be other potential factors to influence sap flow. Although the elevation range in our study area may not be considered very broad, the significance of the elevation in all the models (tables 4 and 5) may additionally emphasize elevation as a very important environmental factor to influence sap flow. Elevation gradients, as the major driving factor of the most environmental variables, may therefore determine how dense trees can naturally grow and survive across different ranges of elevation (table 5 and figures 5 and 6). For example, Mezei et al (2014) found elevation as one of the most important factors for tree mortality caused by the spruce bark beetle (Ips typographus (L.)), using 10 year data in Central Europe. Higher elevation was associated with higher tree mortality in all of the outbreak phases (Mezei et al 2014). Within an elevation range, tree density may further influence individual tree traits such as bark temperature due to crown closures, DBH, and its associated factors, e.g. tree height and crown area, due to limitations in accessing shared resources among trees (table 5, figures 5 and 6).

Can landscape-and stand-level drivers of sap flow represent or compensate for tree-level factors?
The elevation gradients are associated with the subsequent temperature and vapor-pressure gradients that influence sap flow (figures 5 and 6). Higher elevations are generally associated with a steeper landscape, causing higher incidence angles, which eventually decrease the intensity of solar irradiance (figure 6). These counter-playing environmental factors related to elevation gradients, e.g. temperature and VPD, and the intensity of solar irradiance, together influence the amount of sap flow (figures 5 and 6). The stand-level measure of neighborhood tree density, related to the compositions (proportions) of the tree cover and non-tree interspace gaps, was negatively correlated with sap flow (table 2; figures 5 and 6). The increase in tree density may sometime influence the spatial configurations (juxtapositions) of tree cover and non-tree gaps that, for the most part, define the spatial structure of forests (figures 5 and 6). The spatial configurations further define several important elements in the evergreen forests that including illuminated and shadowed trees and illuminated and shadowed gaps within and among tree canopies (Zabihi et al 2021). When the number of trees increases per area, the competition for using shared resources such as water, nutrients, and solar radiation increases, influencing photosynthesis and sap flow. On the other hand, the higher amount of interspace gaps among trees, for the most part, due to lower tree density, allows more exposure to solar radiation causing more heat loads on tree canopies and forest grounds that increase evapotranspiration. However, the low stand density was found to be a very important factor to improve branch diameter (Mäkinen and Hein 2006), and DBH and tree height (Katrevičs et al 2018). In our datasets, the correlation between DBH and, elevation, tree density (table 2) may indicate that elevation and tree density together are driving individual tree traits such as DBH and its associated factors. For example, even though DBH and bark temperature were positively and negatively correlated with sap flow, respectively (table 2), including them in the spatial regression models did not improve the model performances (models no. 1 vs. 2-8; table 4). These findings may indicate that the effects of individual tree traits on sap flow to some extent can be covered by their driving factors such as elevation and tree density. Therefore, measuring tree density ranges within elevation ranges may be used to evaluate forest health, related to individual tree traits, e.g. DBH, height, above-and below-ground biomass, to some extent, waiving to conduct labor-intensive and costly tree-level measurements of sap flow (figure 5). For example, Korolyova et al (2022) found areas with denser trees, likely to limit available resources for trees due to competition effects, negatively impacted Norway spruce survival in post-outbreaks of bark beetle infestations in Central Europe.

Can latitude and longitude be substitutes for space effects in the OLS models?
Including latitude in some of the OLS models resulted in spatially independent model residuals (models no. 13-17) that avoided needing to run spatial regression models. However, latitude and/or longitude were not able to remove spatial autocorrelation effects in models no. 9-12 (table 4). Given the regional scale of our study, the latitude significance in the OLS models (no. 13 and 17) is likely an underlying result of spatial autocorrelation effects rather than the actual patterns in the landscape. Usually, when latitude increases over very broad scales, the solar incidence angle increases, causing cooler temperatures and likely less flow of sap. The longitude significance in the OLS models (no. 14 and 16) is also likely due to the underlying result of spatial autocorrelation effects. Although latitude and longitude gradients accounted for spatial autocorrelation effects in the OLS models to some extent, the results were found to be misleading. For example, comparing models without and with latitude (no. 1 vs. 15, 2 vs. 13, and 4 vs. 17) revealed that latitude effects caused tree density to become an insignificant factor at a 0.05 error rate, likely due to multicollinearity effects. Therefore, if a study is spatially distributed over a very broad scale at which latitude and longitude are variables of interest, including them in an OLS model may account for the spatial autocorrelation effects while revealing patterns associated with latitude and longitude effects. The extent to which latitude and/or longitude may account for the space effects is highly driven by the spatial distribution of samples, e.g. the level of clustered patterns across latitude and/or longitude gradients. On the other hand, the confounding effects of space and environmental patterns associated with latitude and longitude gradients could be complex to tease apart, using an OLS model. Thus, using a spatial regression model is recommended for the situation where latitude and/or longitude are variables of interest while the datasets are spatially autocorrelated.

Can OLS regression models cause misleading results?
As expected, OLS models completely failed to take the spatial autocorrelation into account when latitude and longitude were not considered in the models (no. 1-8; significant p-value of Moran's I; table 4). In addition, some of the models (no. 9-12) having latitude and/or longitude included were not even able to take the spatial autocorrelation into account while other models (no. 13-17) were able to. Results from the top-ranked OLS model (no. 13) were different from the top-ranked spatial regression models (no. 1-3; ∆AIC ⩽ 2). In addition, in models no. 9 and 10, latitude became insignificant in spatial regression models whereas it was found to be significant in the OLS models even though the OLS models were not correct choices (significant p-value of Moran's I; table 4). Therefore, when the existing spatial autocorrelations are not taken into account, OLS models lead to inflated type I error (a false rejection of a true null hypothesis) and wrong inferences due to incorrect estimates of a parameter (coefficient).

What type of spatial regression model to use?
In most of the proposed models, spatial error models performed better than the spatial lag model with a lower AIC value, likely due to excluding more noises/errors associated with the additional estimates of space in the spatial lag models. However, if the inference about the space effects on the response variable is as important as other model predictors, the spatial lag model may be preferred over than spatial error model. The spatial error model assesses the spatial clustering on the error term, allowing the model to capture the effects of unmeasured independent variables.

Conclusions
The local-to regional-scale forest management strategy may be required to improve forest health and resiliency against water-stress stimuli as climate regimes, such as inter-annual and periodic variability in temperature and precipitations, may affect forests globally. For example, resource managers and policymakers may establish efficient silviculture practices and plantation strategies, such as forest thinning and further reforestation with optimum distances, to counterbalance the allocation of resources, e.g. water, nutrients, and light, among trees. As sap flow was spatially autocorrelated, and tree density and elevation were the major drivers of sap flow variability, the distance/spacing thresholds among trees may need to be adjusted based on the elevation gradients.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// datadryad.org/stash/share/HzNx0cUXsFtq1m65OzD4QH8ShnJk2wCTQ5mc-Ug33O4.