Modelling growth of Pinus taeda and Eucalyptus grandis as a function of light sums modified by air temperature, vapour pressure deficit, and water balance

Background: Hybrid mensurational and physiological models seek to combine precision, process explanation, simplicity in parameter definition, and ability to estimate wood products. The aim of this study was to assess the suitability and the advantages of the hybrid mensurational-physiological approach where time has been substituted for light sums in growth equations, to replace traditional time-based models in forecasting systems for Eucalyptus grandis W.Hill and Pinus taeda L. Methods: Using 974 permanent sample plots from plantations in Uruguay, we adjusted growth equations to project dominant height, net basal area, maximum diameter breast height, and standard deviation of diameters as a function of accumulated light restricted by modifiers that account for principal physiological limitations on photosynthesis. We analysed: i) the inclusion of terrain aspect and slope information for computing radiation; ii) the use of modifiers for temperature, vapour pressure deficit and water balance; iii) bias and precision of hybrid models with respect to time-based equations. Results: Growth equations showed a good fit for both species when modelled as a function of light sums modified by vapor pressure deficit, air temperature, and water balance. Accounting for slope orientation when computing light sums did not increase precision. Compared to time-based formulations, hybrid models presented a root mean squared error reduction of 10.7% and 4.5% on average for Eucalyptus grandis and for Pinus taeda, respectively, and the relationship between growth and resource availability was consistent with eco-physiological principles for both species. Conclusions: The hybrid methodology can be applied as a basis of forecasting systems for the species studies with significant advantages over time-based models, such as: (i) an increase in precision; (ii) an increase in spatial and time resolution; and (iii) the possibility of simulating the effect of changes in air temperature and water availability on tree growth. New Zealand Journal of Forestry Science Rachid-Casnati et al. New Zealand Journal of Forestry Science (2020) 50:3 https://doi.org/10.33494/nzjfs502020x17x


Introduction
Hybrid growth and yield models have been developed to include the best features of mensurational and physiological approaches while attempting to avoid the shortcomings of complicated physiological models (e.g. limitations from a forest management perspective).
Keywords: forest modelling, hybrid mensurational-physiological models, modified light sums, stand dynamics, Pinus taeda, Eucalyptus grandis more than two models (Peng et al. 2002;Robinson & Ek 2003). A hybrid mensurational-physiological approach using equations based on potentially useable light sums (PULSE) was proposed by Mason et al. (2007). It combines light use efficiency with mensurational approaches to describe growth using difference sigmoidal models. The PULSE approach consists of the substitution of time for modified light sums in sigmoidal growth equations (Figure 1), to model directly the descriptive variable of interest (e. g. dominant height or basal area). Modified light sums are computed as the monthly light sum multiplied by adapting the modifiers used by the 3-PG model (Landsberg & Waring 1997) as follows: (1) where R t (MJ) is the total light sum from month 1 to T, R is radiation in month T, and f θ , f D , and f T are the soil water balance, vapour pressure deficit (VPD), and temperature modifiers calculated for month T.
The hybrid approach seeks to improve the explanation of the growth process by introducing information about several key factors regulating the proportion of radiation that could be potentially used by stands to grow. This is a major advantage compared to traditional mensurational models, which fail to represent growth under changing environmental conditions. Compared to physiological models, an advantage of this methodology when applied in commercial forests (where leaf area is rarely measured in forest inventories), is that the use of leaf area index (LAI) to estimate intercepted radiation can be avoided by employing total radiation (instead of absorbed PAR).
At an experimental scale, the hybrid modelling approach was first applied to estimate ground-line diameter growth in a Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco.) competition control experiment (Mason et al. 2007). Under several competition situations, parameters of the equation using potentially useable light sum equations (PULSE) proved to be stable in all the cases suggesting that environmental changes were "absorbed" by the modifiers. Salekin (2019) applied this approach to model growth and mortality of Eucalyptus globoidea Blakely and Eucalyptus bosistoana F.Muell.. Light-based equations were more precise than timebased equations and were also able to accommodate micro-site factors influencing growth and survival of both species. In a similar approach, Montes (2012) modelled increments in height and basal area, and mortality as a function of absorbed PAR (APAR) sums restricted by the same modifiers (temperature, soil water, and vapour pressure deficit) for Pinus taeda L., using a state space approach developed by García (1984). Their study showed good agreement between predicted and observed data at an experimental scale.
Studies at a regional scale are scarce. Mason et al. (2011) documented a gain in precision for basal area modelled applying PULSE with respect to time-based models for Pinus radiata D.Don plantations in New Zealand. However, for dominant height improvements were modest. In Sweden, the approach was applied to predict site index for Pinus sylvestris L. (Mason et al. 2018). The authors assessed several variables and found that modified light sums was the one that best explained site index.
Results of applying PULSE seem promising with respect to precision and outputs for predicting height and especially basal area, however there is still much to explore with respect to the potential performance of this approach, especially seeking its full application as a basis of forecasting systems. Analysis of the methodology in different species at regional scale are rare, and its application for modelling other variables used in projection systems, such as diameter variation, are unknown. We hypothesise that models for projecting maximum diameter at breast height and standard deviation of diameters could be improved when using potentially useable radiation sums. Previous studies showed that the inclusion of variables related to site quality or site characteristics, improved the precision of models for pine and eucalypts (Methol 2008(Methol , 2001Rachid-Casnati et al. 2019).
Eucalypts and pine plantations occupy 9,102,300 ha and 4,432,000 ha respectively in Brazil, Argentina, Chile and Uruguay (Di Marco 2014;IBGE 2017;Gysling et al. 2018;MGAP-DGF 2018). In Uruguay, Eucalyptus grandis W.Hill is the most planted species (over 230, 000 ha) managed in long rotations for saw logs and plywood or short rotations for pulp. Pinus taeda is the most planted pine species (nearly 180,000 ha) (MGAP-DGF 2018) for producing saw logs and plywood. Growth and yield systems predicting the effect of factors limiting productivity are of increasing importance for intensively-managed plantations of these species, and their development has been an increasingly important research goal (Soares & Leite 2000;Maestri 2003;Temps 2005;Carrero 2012;Munhoz 2016;Rachid-Casnati et al. 2019;Schons et al. 2019;Scolforo et al. 2019). However, the performance of PULSE applied to those species has not been analysed. FIGURE 1. Sigmoidal relationship between yield (y) and modified light sums at age t (R t ), where modified light sum is monthly light sum (R) in the top of the canopy multiplied by air temperature modifier (f T ), and the minimum between VPD modifier (f D ) and water balance modifier (f θ ) accumulated from plantation to age t.
The aim of this study was to assess the suitability of potentially useable light sums modified by water balance, vapour pressure deficit, and air temperature, to model four state variables (dominant height, basal area, maximum diameter, and standard deviation of diameters) that account for stand growth and structure in two contrasting and widely planted species: Pinus taeda and Eucalyptus grandis. For this purpose, we established three objectives: (1) To adjust sigmoidal equations for modelling dominant height (h dom ) , net basal area (G), maximum diameter (d max ), and standard deviation of diameters (SD d ) using different light sum formulations, to select the best alternative.
(2) To test the feasibility of using hybrid models to simulate the effect of changes in temperature and rainfall regimes on stand growth.
(3) To assess the advantages of hybrid mensurationalphysiological equations with respect to time-based counterparts for each state variable by comparing precision and bias from both sets of equations.

Study location and sources of data
The study was conducted in the northern half of Uruguay, which is the region of the country with the highest concentration of commercial forestry ( Figure 2). Data describing tree growth were obtained from 974 geo-referenced permanent sample plots (PSP) ( Climate data were obtained from 28 meteorological stations and 57 rain gauges (Figure 2). Climate and soil data were provided by the Uruguayan Institute of Meteorology (INUMET), Group of Agroclimate and Information Systems (GRAS) of INIA Uruguay, and the Ministry of Cattle, Agriculture and Fisheries (MGAP). Monthly average values for daily maximum and minimum temperature and sunshine duration from each meteorological station as well as accumulated precipitation from the rain gauges were interpolated into grids of 500 m x 500 m for each month corresponding to the time ranges of plots used in the study. For rainfall, the interpolation method applied was Inverse Distance Weighting, as it is a simple method which give similar results compared to Thin Plates Splines and Co-kriging (Hartkamp et al. 1999). However, these authors found that Thin Plates Splines were more effective for interpolating temperature, and were therefore used for this variable. Distribution of average monthly rainfall and soil water holding capacity are depicted in Figure 3. Data on aspect, slope, and elevation of plots were obtained from a digital elevation model with a resolution of 30 m x 30 m (MGAP-DNRNR 2016).

Calculation of modified light sums
Global horizontal radiation was computed using the Angstrom equation and applying parameters adjusted for Uruguay by Abal et al. (2010). Radiation was decomposed into direct and diffuse radiation considering tilted surfaces as proposed by Revfeim (1978) and Erbs et al. (1982), respectively. A detailed description of the procedure is provided in Additional File 1.
Radiation was calculated in two ways: i) considering the tilt of a stand's terrain (using aspect and slope of the terrain), and ii) assuming that stands are on a flat surface, in order to assess the utility of this site information for improving predictions. Monthly restrictions of light were applied using modifiers for water balance (f θ ), and VPD (f D ) in the same way as used by the 3-PG model (Landsberg & Waring 1997), whereas the modifier for temperature (f T ) was adapted for each permanent sample plot available in the time range of the study (1979 and 2012). Seven alternatives combining VPD, temperature and water balance modifiers, were compared. Formulations included light sums modified by: water balance, R θ ; VPD, R D ; temperature, R T ; the minimum between water balance and VPD (referred here as waterMod), R θD ; water balance and temperature, R θT ; VPD and temperature, R DT ; and all three modifiers, R θDT .
The modifier that considers vapour pressure deficit (VPD) describes a relationship where the modifier declines exponentially when VPD increases. It is TABLE 1. Summary of the dataset used in the analysis. computed following Landsberg and Waring (1997): where k g is a coefficient based on the relationship between stomatal conductance and VPD. VPD is calculated as follows: (3) where DT max and DT min are the saturated vapour pressure when temperature = T max or T min , respectively. These variables are calculated using minimum or maximum temperature per month (T i ) using the following equation: (4) The soil water-dependent modifier was calculated as (Landsberg & Waring 1997 where c θ is the coefficient and n θ is the power of the water modifier. Both take different values for different soil types. r θ is moisture ratio calculated as: being θ T the soil water balance, and is the potentially available soil water. Information provided by the map of PASW developed for Uruguay by Molfino and Califra (2001) (see Figure 3) was used in order to account for the maximum available water in the root zone.
Monthly soil water balance was estimated through the equation: where θ T-1 is root zone water balance in the FIGURE 3: Spatial variation of potentially available soil water (Molfino & Califra 2001), average for the studied period of accumulated annual rainfall, and plot location.
previous month; R is rainfall; I is canopy interception; E is evapotranspiration and D is drainage. When θ T-1 + R − I − E > PASW, there is excess water assumed to be lost as drainage. The 3-PG model estimates evapotranspiration using the Penman-Monteith "big leaf" model: where λ is the latent heat of water vaporisation (J kg ᵒC -1 ); S is the slope of saturation vapour pressure curve for water (kPaᵒC -1 ); R n is net radiation absorbed by canopy (J m -2 month -1 ); ρ a is air density (kg m -3 ); VPD is vapour pressure deficit (mbar); γ is the psychrometric parameter (kPa ᵒC -1 ); g b is boundary layer conductance (m s -1 ); and g c is canopy conductance (m s -1 ).
Boundary layer conductance depends on wind speed as well as the size and shape of leaves, and density of foliage (Landsberg & Sands 2011), however fixed values are commonly used for practical purposes. Based on the work of Mielke et al. (1999) a fixed value of 0.2 m s -1 was assumed. Mielke et al. (1999) observed wind speeds around 2 m s -1 leading to canopy conductance values of 0.2 m s -1 for E. grandis. According to the National Institute of Meteorology of Uruguay (INUMET), the historical average wind speed in the study area was 3.5 m s -1 . Martin et al. (1999) observed that boundary layer conductance did not increase markedly when wind speed ranged from 1 to 2 m s -1 . Therefore, it was assumed that wind speed was spatially and temporally uniform, and boundary layer conductance values assumed in this study did not seem to lead to significant error.
Canopy conductance was calculated as: (9) where g Sx is maximum stomatal conductance, assumed as 0.02 ms -1 (Almeida et al. 2004;Sands 2004), L is leaf area index, L gC is leaf area index at maximum conductance, and other terms are as specified earlier.
The value for L gC was assumed to be 3.33 (Sands 2004). Because L measurements were not available, it was assumed L/L gC = 1 to calculate canopy conductance.
Net radiation was estimated using a linear relationship with radiation (total shortwave radiation) as it follows: where H s is global horizontal radiation (W m 2 ), q a and q b are the intercept and the slope parameters. The values applied for P. taeda were those used in 3-PG by Landsberg and Sands (2011).  registered much lower intercept values for E. grandis growing in Brazil and since the sensitivity to this parameter is high, their value was applied for this species (see Table 2).
The temperature-dependent growth modifier was based in the assumption that production increased with increasing temperature and started declining after an optimum was reached (Landsberg & Waring 1997): where T a is mean temperature for each month, T min , T max , and T opt are minimum, maximum, and optimum temperature required for growth, respectively. In this study the mean daytime temperature was employed instead of mean temperature as Mason et al. (2011) found that this modification gave better precision than daily mean temperature.
Accumulated radiation for each month from time of planting was multiplied by the different combination of modifiers (temperature, VPD, and water balance): (12) where R θDT is the total radiation sum (MJ), and f θ , f D , and f T are the soil water balance, vapour pressure deficit, and temperature modifier previously calculated for month t. Light sums and modifiers calculations were programmed in the R software (R Development Core Team 2014).

Development of hybrid growth equations
Two sigmoidal differential expressions corresponding to polynomial forms of Schumacher (Eq. 13) and von Bertalanffy-Richards (Eq. 14), selected based on previous time-base studies (Rachid-Casnati et al. 2019), were fitted for each state variable of interest.
where Y 1 and Y 2 are yield at time 1 and time 2, respectively; R t1 and R t2 are modified radiation sums at time 1 and time 2, respectively, a is the asymptote parameter, and b, and c are shape parameters. Seventy percent of the plots were used for fitting while thirty percent were utilised for validation using all possible age intervals per plot. The inclusion of new variables was tested using a correlation-free dataset, where only one interval per plot was included. In this way, serial correlation that could lead to loss of precision of the estimates and invalidation of tests of significance was avoided (Rawlings et al. 1998).
For each state variable (dominant height, h dom ; basal area, G; maximum diameter, d max ; and standard deviation of diameters, SD d ), seven restricted light sum alternatives (previously described) were computed. Each one of these alternatives was developed depending on whether the net radiation calculation considered slope and aspect, or not. Candidate equations were compared through their root mean square error (RMSE), as a measure of precision to select the best expression for modelling each state variable. The best ranked model for each variable was selected. Normality was analysed graphically through histograms and Q-Q plots, while plots of residuals against the fitted and the independent variables were also assessed. For the validation stage, plots of residuals versus predicted values were examined in order to detect bias. Confidence intervals for the slopes of predicted versus independent values with a level of significance of 5% were computed to compare the actual slope with the ideal slope of 1 (Goulding 1979).

Testing the effect of changes in site conditions on growth
The behaviour of hybrid mensurational-physiological models was assessed by simulating growth under different combinations of rainfall and temperature regimes. To compare scenarios based on realistic but contrasting conditions, data from three areas of the country (Figure 4) with a range of minimum and maximum temperature and rainfall, were used. In this way, growth under historical conditions (central location) was compared to growth under increased temperature from a western plot, and increased rainfall from a northern plot. Monthly climate information for a 17-year period was applied, and the series was repeated  Table 3.
The R software (R Development Core Team 2014) was used for all statistical calculations. Procedures involving georeferenced information were developed using QGIS (QGIS Development Team 2015) and ArcGIS (ESRI 2013).

Comparing precision and bias of the hybrid approach with time-based equations
A comparison of the hybrid approach with timebased formulations in terms of precision and bias was undertaken given that an identical dataset was used by Rachid-Casnati et al. (2019) to model the same state variables through time-based equations. In this other study, two formulations: (1) simple (based solely on time); and (2) augmented equations were developed. The latter used physiographic and soil-based predictors namely: water holding capacity, aspect, slope, and elevation. Comparisons were undertaken using the validation dataset. The comparison of precision was made through the root mean square error (RMSE) and the comparison of bias through mean absolute bias (MAB).

Modified light sums
For each plot, the monthly average light sums was calculated using the three modifiers (i. e. R θDT ) and categorised into three classes per species (Figure 5). Regional differences in the amount of light that could be utilised were found for both species with higher values at the north of the country compared to the west and centre. This corresponds to larger accumulated rainfall per year, and generally higher water holding capacity of the soils in this area compared to the other studied sites (see Figure 3).

Growth and yield equations
Previous studies showed that potentially available soil water influences the asymptote of dominant height and basal area growth equations for E. grandis and P. taeda, and maximum diameter growth equations for E. grandis (Rachid-Casnati et al. 2019). Therefore, this possibility was also assessed for light-based equations. This was done by adding a water-related index computed as a 30-year average of monthly light sums constrained by waterMod, R θD . In this way, R θD has a specific value for each plot location. R θD was found to be significant in the models for all the state variables, except for SD d for both species, when analysed using the non-correlated hypothesis testing dataset. For E. grandis, the inclusion of the water-related index in the asymptote gave logical model behaviour and led to a substantial improvement of the models. However, for P. taeda no clear effect was observed, and the water-related index was not included. Therefore, all the state variables for this species were modelled using a fixed asymptote. For net basal area and maximum diameter, the polymorphic form of the Schumacher equation with two parameters showed the best fit, whereas for dominant height and standard deviation of diameters, the von Bertalanffi-Richards equation was selected (Table 4). Residuals plots (see Additional file 2) did not show any strong pattern of bias.

Differences between light sums formulations
PULSE fitted using all three modifiers (R θDT ) showed little improvement over formulations that only applied the water-related factor (R θD ). However, the full formulation was chosen since it comprised a more complete description of growth conditions. Table 5 shows the top two formulations with smallest RMSE, and the formulation with largest RMSE for each state variable. Equations using modified light sums applying only the temperature modifier were the least precise in general for both species with significant differences in RMSE. Very small differences were found in the errors associated with the different alternatives for computing radiation based on slope and aspect. Parameters estimated for all the fitted equations are presented in Additional file 3.

Effect of changes in site conditions
When testing the models combining two different levels of temperature and rainfall for both species, we found that growth was higher for the environmental conditions combining higher rainfall and lower temperature (1480 mm and 17.5 o C), and lower for the combination of lower levels of rainfall and higher temperature (1230 mm and 19.1 o C). In turn, the highest growth rates in P. taeda were obtained at lower temperatures, regardless of FIGURE 5: Permanent sample plots categorised according to monthly averages of accumulated potentially useable light sums per species.
--- R θDT1 and R θDT2 : radiations modified by water balance, VPD, and temperature modifiers and aggregated up to t 1 and t 2 respectively R θD : 30-year average of monthly radiation for the site modified by water balance and VPD.
-RMSE: root mean square error for each equation selected  (Figure 6), while in E. grandis they were obtained at higher rainfall, regardless of temperature (Figure 7). Moreover, in E. grandis differences between conditions were small for the standard deviation of diameters and marked for maximum diameter. For P. taeda, differences between regimes were smaller compared to the eucalypt species, and minimal for basal area when stand age approached 25 years. For E. grandis, larger contrasts in growth under different water regimes were defined by the water factor in the asymptote for all the state variables but SD d . Given that the pine species did not include such descriptor in the asymptote, response of the species to water availability was less significant.

Comparing precision and bias of hybrid approach with time-based equations
Comparison of root mean square error (RMSE) and mean absolute bias (MAB) of the three approaches is depicted in Tables 6 and 7, respectively. Differences in precision and bias were observed among the modelling approaches taken and between the two species. For Pinus taeda, base mensurational equations in general showed the largest errors, except for dominant height. For this variable, the hybrid approach was less precise with a small difference compared to the augmented time-based approach (3.4%). For all the other state variables, the hybrid methodology was the most precise; especially for modelling net basal area and maximum diameter (precision gains were 7 and 9 %, respectively over the base mensurational approach). The hybrid formulation provided the smallest calculated bias except for dominant height. The improvements achieved by modelling all state variables using hybrid models were considerable for Eucalyptus grandis; the greatest improvement over the base mensurational approach was more 14% for basal area. According to the MAB rank, this approach was the least biased as well. For both species and considering all the state variables analysed, the hybrid approach was predominantly first in the rankings for precision and bias with the exception of dominant height for Pinus taeda (rank 3).

Discussion
In our study, changes in dominant height, basal area, maximum diameter, and standard deviation of diameters, were modelled as functions of accumulated light that was potentially available after considering cardinal temperature for a species´ growth, available soil water, and vapor pressure deficit for each month during the lifetime of a stand. It was assumed that this relationship followed a sigmoidal curve, where functions comprise a growth component and a decline component (Zeide 1993), as typically considered in mensurational models. Simultaneously, the modelling method ensured that mathematical properties such as path invariance and consistency (Clutter 1963) are complied with, guaranteeing robustness to the forecasting system.
As opposed to time-based models, the explanatory variable (potentially available light sums) considered the structural characteristics of sites (soils characteristics, aspect, slope, etc), and climatic information. Site information changes spatially, whereas climatic information changes spatially and temporally. Consequently, the independent variables were not directly measured, as is the case with age, and were therefore susceptible to error. However, summarising several growth factors that interact and vary with location and plantation period is an advantage over time-based formulations that consider only one growth factor (Woollons et al. 1997a;Maestri 2003;Pinjuv et al. 2006;Scolforoet al. 2019a) or that consider surrogate variables (Rachid-Casnati et al. 2019).

Models formulations and use of growth modifiers
The simultaneous inclusion of the three modifiers provided the best formulation for projecting growth, which reflected the importance of VPD, air temperature and water balance on growth and yield. The modifiers that contributed the most to gains in precision were the water-related modifiers. This was observed by Campoe et al. (2016) for both species. For P. taeda, water and temperature are key drivers of productivity (Teskey et al. 1987), although across its natural range soil nutrient availability is considered to be one of the most influential factors (Hebert & Jack 1998;Allen & Albaugh 1999;Jokela et al. 2004). The sensitivity of P. taeda to VPD, soil water content, and photosynthetically active radiation was verified by Ford et al. (2004). The sensitivity to the first two factors was indicated by Manogaran (1973) and confirmed in Uruguay by Gándara et al. (2014). For E. grandis hybrids, soil water content and atmospheric humidity was shown to be the most important growth modifiers in Brazil (Almeida et al. 2007;Stape et al. 2008Stape et al. , 2010. Ryan et al. (2010) demonstrated that wood primary productivity increased due to increases in light interception, as well as photosynthetic efficiency, and partitioning to wood, when eucalypts were irrigated. An extensive review documenting the eco-physiology of productivity in eucalypts was presented by Whitehead and Beadle (2004) who emphasised the enormous increase in light use efficiency that the species can achieve in well-watered soils. However, recent studies conducted by Scolforo et al. (2019b) stressed that the impact of water deficit can vary among eucalyptus clones, uncovering the increasing difficulties that tree breeding imposes on the generalisation of growth models. The results of this study with respect to the relationship between tree growth and temperature, VPD, and soil moisture are coherent with physiological principles reported in the literature for both species. Basing growth formulation on biological hypotheses improves generalisation capacity of models and increases model utility (Maestri et al. 2013).
Considering slope and aspect for computing radiation Coops et al. (2000) showed that differences in incoming radiation among surfaces tilted to a variety of slopes and different orientations are worth considering when estimating incoming radiation for its application in forest eco-physiology. The study also showed that largest slopes gave larger estimation errors when  Growth curves for the modelled variables as a function of time calculated using light sums obtained for the range of growth conditions analysed for Eucalyptus grandis.
reflected radiation was not considered. In this study, the methodology used was the one proposed by Tian et al. (2001), who also registered differences in estimated global radiation considering slopes and aspects at a daily time step. Computing radiation considering a tilted surface did not improve the quality of the model fit in our study. However, applying slope and aspect information in time-base equations improved precision for E. grandis when modelling dominant height, standard deviation of diameters, and maximum diameter, and for P. taeda when modelling basal area in a previous study (Rachid-Casnati et al. 2019). The fact that the maximum slope value found in this study was 10 degrees and most of the plots were located on sites with gentle slopes could have contributed to these results. Therefore, it is possible that the methodology loses sensitivity in flatter terrain. The digital terrain model from which slope and aspect information were extracted had a medium resolution (30 m x 30 m). Consequently, considering that slopes are gentle, resolution of the digital terrain model combined with radiation calculated using sunshine duration, could have contributed to blur possible effects of slope and aspect on tree growth. This should be tested in the future by applying a finer-scale digital terrain model.

Time-based models vs hybrid models
Precision gain compared to time-based models was substantial for most state variables, however there were differences in gains among variables. Precision gain was larger in basal area compared with dominant height for both species, consistent with observations in previous studies where a similar methodology was applied (Mason et al. 2011;Montes 2012), and in studies where traditional models were augmented with site variables (Woollons et al. 1997b;Pinjuv et al. 2006;Rachid-Casnati et al. 2019). The hybrid approach was more precise with respect to stand structure changes, since maximum diameter and standard error of the diameter were also more precisely modelled with the hybrid formulations. Relationships between site quality and maximum diameter have been verified in Uruguay for E. grandis (Methol 2001) and P. taeda, with a positive effect on growth. A similar relationship between site index and standard deviation of diameter was found for E. grandis (Methol 2003) and E. globulus (Methol 2006). For both species, an increase in site index led to an increase in diameters variation in Uruguay. Sanquetta et al. (2014) found that parameters that characterised the shape of a Weibull distribution for modelling diameter variation in Acacia in Brazil were more precisely estimated when rainfall was included. With increasing rainfall, populations had larger mean diameters and greater diameter variability. Those results agree with our findings. However, the water-related index was not significant for modifying the asymptote in the function to predict the standard deviation of diameters. For eucalypt species this was the only modelled variable that did not include the water-related index, since it was significant when modelling dominant height, basal area, and maximum diameter. We hypothesise that competition is a significant factor influencing standard deviation of diameters, and information about potential share of resources among individuals may be as relevant as the general resource availability of sites in explaining tree size variability. According to Burkhart and Tomé (2012), growth rates of individual trees are influenced by micro-environment, genetics, stand density, and average growth potential influenced   (1) by local neighbours. In the hybrid models, information about stand density could be incorporated in the future to improve hybrid formulations for modelling standard deviation of diameters. This information has been reported as significant when included as a predictor of the asymptote with a negative relationship in E. grandis in Uruguay (Methol 2003), and Pseudotsuga menziesii and Pinus radiata in New Zealand (Methol 2001). In our study, this was not assessed. There are no (known) reports of the use of light sums for modelling growth of maximum diameter or standard deviation of diameters. Both variables are important for estimating diameter distributions if using the inverse Weibull distribution (Kuru et al. 1992), which is currently applied in Uruguay. Using resource-driven diameter distributions would show whether changes in growth factors affect stand structure and how those changes influence volumes of assorted products. It must be considered that management practices, such as thinning, influence stands structure greatly, and algorithms that modify the frequency of diameter classes should be accordingly introduced in forecasting systems. Hybrid formulations were more precise compared to time-based models, however there were differences between species. For E. grandis, average precision gain almost doubled that of the P. taeda for most of state variables, except for maximum diameters, which showed similar increases in precision between species. A potential hypothesis to explain these results is the lack of information of mature stands (more than 16 years) for accurately modelling the asymptote for the pine species. The water-related index did not show a clear effect when included in the asymptote for P. taeda, as opposed to its inclusion in the asymptote for E. grandis. For the latter, it was statistically significant for three of the four state variables modelled with a clear positive effect and it was responsible for a large part of the improvement in precision. For P. taeda, past research stressed the importance of soil fertility on growth (Albaugh et al. 2004(Albaugh et al. , 2005. Munhoz (2016) concluded that soil variables had greater influence on the use-efficiency of natural resources than climate variables for P. taeda plantations growing in subtropical areas of Brazil. Therefore, future studies should tackle this question in Uruguay to determine whether a fertility modifier should be included in hybrid models. A longer age-span of P. taeda data should help us to understand whether those differences are related to the species' requirements or to the data used in this study, or both.

Model applications
When assessing the behaviour of the model for contrasting combinations of rainfall and temperature, the approach showed sensitivity to a change of 2 ᵒC in average temperature and 250 mm of yearly rainfall, with a positive effect for rainfall increase and temperature decrease on growth for all response variables. These results agree with those reported by Binkley et al. (2017), who also obtained a positive relationship between rainfall and growth, and a negative relationship between temperature and growth when considering an average response of several Eucalyptus clones planted across Brazil. Our analysis also showed some differences between species regarding the effect of climate factors on growth for the ranges of the values studied. The inclusion of a 30-year monthly average of light sums restricted by water balance or VPD on the asymptote, contributed to increasing growth differences at the end of the projection periods between contrasting regimes for E. grandis.
The analysis undertaken served as an example of the possibilities of the hybrid equations demonstrating their capacity to accommodate the effects of soil moisture, VPD, and air temperature on growth. This offers innovative opportunities for forest management such as: i) assessing the effect of changes in factors such as climate change, or irrigation on growth; ii) assessing the consequences of applying tillage techniques that modify the water holding capacity of soils (Mason 2013); iii) analysing possible effects of changing growth factors on the structure of stands; iv) undertaking within-year analysis to assess the effect of seasonality on growth which could affect short-rotation plantations (Campoe et al. 2016); v) classify sites according to limiting factors; and vi) analyse the influence of microsite on growth differentiation, to plan differential management (González-Barrios et al. 2015;Salekin 2019).
We also identify possible shortcomings of this methodology, viz: i) the light-based approach would rely on historical climate data when applied by managers for projecting future scenarios; ii) making the models available to forest managers would necessitate a more complex information system to provide meteorological and soil-based information, as well as computing restricted light sums based on date of planting and extent of growth projection; and iii) precision of formulations rely on the quality of the information needed to compute modified light sums.

Conclusions
The hybrid mensurational-physiological approach presented here was clearly more precise than timebased equations. It also showed consistency with ecophysiological knowledge of both species, regarding the relationship of growth with air temperature, VPD, and soil water availability.
Modelling growth of maximum diameter and standard deviation of diameters as a function of modified light sums, enables the system to assess changes in population structure with changes in resource availability.
Slope and aspect information were not useful in this framework, however, assessing the methodology using a finer-scale digital terrain model is recommended before discarding its usefulness.
Future research should seek to model growth using a full life-span plots for P. taeda to understand the importance of resource availability near the asymptote age. The inclusion of a modifier for fertility could be assessed to understand its impact on model improvement.
Equations offered in this work can be included in a growth and yield modelling system that would be able to provide accurate plot production projections under changing growth conditions with higher temporal and spatial resolution. This greatly improves utility of models traditionally used by foresters, where usual input variables are only dbh, dominant height, and stand density at the start of the projection period. In the proposed approach, the only extra information that needs to be provided is the geographical locale of the stand.

Funding
Funding for this study was provided by the Ministry of Foreign Affairs and Trade, through NZAid Scholarship, and INIA Uruguay.