Modeling tree radial growth in a warming climate: where, when, and how much do potential evapotranspiration models matter?

Process-based models of tree-ring width are used both for reconstructing past climates and for projecting changes in growth due to climate change. Since soil moisture observations are unavailable at appropriate spatial and temporal scales, these models generally rely on simple water budgets driven in part by temperature-based potential evapotranspiration (PET) estimates, but the choice of PET model could have large effects on simulated soil moisture, moisture stress, and radial growth. Here, I use four different PET models to drive the VS-Lite model and evaluate the extent to which they differ in both their ability to replicate observed growth variability and their simulated responses to projected 21st century warming. Across more than 1200 tree-ring width chronologies in the conterminous United States, there were no significant differences among the four PET models in their ability to replicate observed radial growth, but the models differed in their responses to 21st century warming. The temperature-driven empirical PET models (Thornthwaite and Hargreaves) simulated much larger warming-induced increases in PET and decreases in soil moisture than the more physically realistic PET models (Priestley–Taylor and Penman–Monteith). In cooler and more mesic regions with relatively minimal moisture constraints to growth, the models simulated similarly small reductions in growth with increased warming. However, in dry regions, the Thornthwaite- and Hargreaves-driven VS-Lite models simulated an increase in moisture stress roughly double that of the Priestley–Taylor and Penman–Monteith models, which also translated to larger simulated declines in radial growth under warming. While the lack of difference in the models’ ability to replicate observed radial growth variability is an encouraging sign for some applications (e.g. attributing changes in growth to specific climatic drivers), the large differences in model responses to warming suggest that caution is needed when applying the temperature-driven PET models to climatic conditions with large trends in temperature.


Introduction
Tree-ring widths are used both to infer past climate variability and change (Fritts 1976, Jones et al 2009 and to detect climate influences on forest productivity (Williams et al 2013. Doing so requires a model linking tree radial growth to climate (Fritts 1976), with models ranging in complexity from strictly empirical statistical models to mechanistic models of tree-ring formation via known or hypothesized environmental limitations to cell division and expansion in the vascular cambium (e.g. Vaganov et al 2011, Anchukaitis et al 2020).
The process-based models generally simulate radial growth responses to solar radiation, temperature, and soil moisture, but direct soil moisture observations are not usually available at appropriate spatial and temporal scales for model simulations. They therefore rely on simple water budget estimates of rootzone soil moisture, partitioning inputs of precipitation (P) into losses of moisture by evapotranspiration and runoff and infiltration of moisture into one or more soil layers (Huang et al 1996, Tolwinski-Ward et al 2011.
To drive the atmospheric demand side of the water balance, water budget-based soil moisture models depend on reliable estimates of potential evapotranspiration (PET), but there are several common methods for doing so, each with varying degrees of complexity and physical realism resulting in different spatial patterns, mean conditions, seasonality, and long-term trends (Fisher et al 2011, Sheffield et al 2012, Trenberth et al 2013. However, the extent to which the choice of PET model influences tree growth simulations and responses to climatic forcing have not yet been established. The semi-empirical VS-Lite model (Tolwinski-Ward et al 2011), for example, uses PET from the Thornthwaite (Th) method (Thornthwaite 1948, Thornthwaite andMather 1955) to drive the atmospheric demand side of the water budget. An empirical model based solely on monthly mean air temperature, the Th model assumes covariation of temperature with the physical drivers of evapotranspiration (net radiation, wind speed, and vapor pressure deficit), which is reasonable for spatial patterns and short-term variation of PET but which is violated for long-term trend assessment (Sheffield et al 2012). The Th model therefore tends to overestimate PET responses to warming temperatures, leading to spurious trends in drought indices and water balance-based soil moisture estimates under warming (Sheffield et al 2012, Trenberth et al 2013, Berg and Sheffield 2018. While the influence of PET models on drought indices has been well established, it is not clear how these differences would translate to the growth simulations of tree-ring models since they also depend on both the magnitude and seasonality of PET change and are mediated by effects of solar radiation and temperature on growth. Establishing how much the PET model choice influences both accuracy and trends of ring-width simulations is therefore a necessary step for interpreting modeled growth responses to warming. Here, I examine how choice of PET model affects moisture stress and growth simulations in the VS-Lite model using more than 1200 tree-ring chronologies in the conterminous United States. I specifically compare Th-driven water balance estimates to those derived from three other PET models of varying complexity: the Hargreaves (Hg), Priestley-Taylor (PT), and Penman-Monteith (PM) models. I examine two aspects of VS-Lite performance with each PET model: (a) the ability to replicate interannual variability in observed ring widths, and (b) the sensitivity of VS-Lite moisture stress and growth simulations to warming. I hypothesize that growth simulations based on more physically realistic PET models (e.g. Penman-Monteith) will outperform simpler empirical models in their ability to capture the variability of observed ring widths. I also hypothesize that moisture stress responses to warming will differ among the PET models, with Th-driven radial growth decreasing most strongly with increased warming and radial growth simulations driven by the more physically realistic PET models (PM and PT) decreasing the least.

VS-Lite model
The VS-Lite model, a semi-empirical simplification of the process-based Vaganov-Shashkin model (Vaganov et al 2011, Anchukaitis et al 2020, simulates monthly growth responses to temperature (g T ), soil moisture (g M ), and solar radiation (g E ) (Tolwinski-Ward et al 2011). The main model inputs are site latitude, monthly mean temperature, and monthly total precipitation, from which soil moisture is estimated using a simple single-layer 'leaky' bucket water balance model (Huang et al 1996). Growth responses to temperature and soil moisture are then defined as piecewise linear functions ranging from 0 to 1 (figure 1): where T mean and M are monthly mean temperature and soil moisture (respectively), T 1 and M 1 are temperature and soil moisture thresholds below which growth cannot occur, and T 2 and M 2 are thresholds beyond which temperature and soil moisture are no longer limiting and do not stimulate additional growth. Consistent with the Principle of Limiting Factors and Liebig's Law of the Minimum (Fritts 1976), a monthly overall growth response (g) for a given site is calculated as the minimum of g T and g M , scaled by incoming solar radiation (g E , the ratio of monthly mean daylength to daylength at summer solstice): The growth index, g, is then annually integrated over a fixed growing season, by default starting in September of the prior year (reflecting dependence of growth on conditions late in the prior growing season) through December of the growth year, and standardized to a mean of zero and unit variance over the simulation period.

PET models 2.2.1. Thornthwaite (Th) model
Originally designed for climate classification (Thornthwaite 1948), the Th model is based on Figure 1. Hypothetical VS-Lite temperature (a) and soil moisture (b) response functions (g T and gM, respectively). T1 and M1 represent lower temperature and soil moisture thresholds (respectively) below which growth cannot occur, while T2 and M2 represent upper thresholds beyond which no additional growth is stimulated (i.e. where temperature or soil moisture conditions are optimal for growth).
empirical relationships between PET, daylength (L), and monthly mean temperatures (T mean ): where PET Th is Thornthwaite PET in mm month −1 , N is number of days in the month, I is a heat index based on a sum of scaled climatological monthly mean temperatures (T i ), and α is an empirical constant based on the heat index, I:

Hargreaves (Hg) model
The Hg model Samani 1985, Hargreaves andAllen 2003) estimates PET from potential solar radiation at the top of the atmosphere (R a ) (estimated daily and summed over each month) scaled by monthly mean (T mean ), mean minimum (T min ), and mean maximum (T max ) temperatures, all in degrees Celsius: The difference between T max and T min implicitly accounts for reduction of solar radiation due to clouds since the diurnal temperature range tends to decrease with increasing cloud cover (Hargreaves and Allen 2003). The model estimates potential latent heat flux in the same units as R a (MJ m −2 month −1 ), which I converted to mm month −1 of PET based on the latent heat of vaporization.

Priestley-Taylor (PT) model
The PT model (Priestley and Taylor 1972) represents the thermodynamic drivers of PET based on the slope of the saturation vapor pressure curve (∆) at mean monthly air temperature and the monthly net radiation, R n , available to drive evapotranspiration: where γ is the psychrometric constant (estimated from mean expected surface air pressure at a given elevation) and a is an empirical constant. I estimated R n as the sum of net solar radiation (R s ) and net longwave radiation (R l ) at the surface. Following Allen et al (1998), R s was estimated as: where α is surface albedo of a well-watered reference crop (set at 0.23 following Allen et al 1998), and k Rs is an empirical constant (0.19 for coastal locations and 0.16 for continental locations). Net longwave radiation was estimated from T max and T min using Stefan-Boltzman's Law, with surface and atmospheric emittances estimated from surface vapor pressure (using Teten's formula and monthly mean dewpoint temperature, T d,mean ) and the ratio of R s to potential clearsky solar radiation (Allen et al 1998, Campbell andNorman 1998). Like the Hg model, the PT model estimates latent heat flux, which I converted to PET in mm month −1 .

Penman-Monteith (PM) model
The PM model represents both thermodynamic and aerodynamic drivers of evapotranspiration, including net radiation, temperature, vapor pressure deficit (VPD), and aerodynamic and stomatal resistance to vapor transport (Penman 1948, Monteith 1965. Here, I used the simplified FAO-56 PM reference evapotranspiration model (Allen et al 1998): where e s and e a are monthly mean saturation and actual vapor pressure (in kPa), G is ground heat flux (in W m −2 ), and u 2 is monthly mean wind speed two meters above the surface (in m s −1 ). I estimated R n , γ, and ∆ using the same approaches described in section 2.2.3. Monthly mean e s was estimated from Teten's formula as the average of saturation vapor pressures at monthly mean maximum and minimum air temperatures (Campbell and Norman 1998). Actual vapor pressure (e a ) was calculated from monthly mean dewpoint temperature using Teten's formula. I estimated mean ground heat flux for each month based on the difference between mean air temperatures of the preceding and following months: Wind speed is not widely available from gridded meteorological products, including the climate data used here (section 2.4 below), but PET estimates using the PM model are not particularly sensitive to this term (Allen et al 1998). I therefore estimated PET PM assuming 1 m s -1 average wind speed.

Tree-ring data
To calibrate and evaluate the VS-Lite model, I obtained measured tree-ring widths from the International Tree Ring Data Bank (ITRDB) for all sites in the conterminous US with elevation estimates provided in the site metadata (required for the PT and PM models) and with complete measurements from 1895 to 1980, resulting in a total of 1279 tree-ring width chronologies covering more than 70 species (figure 2). Raw ring widths were detrended using a stiff 2/3 smoothing spline and averaged to site-level chronologies using Tukey's biweight robust mean in the dplR library (Bunn 2008) of the R statistical programming environment (R Core Team 2021).  (Daly et al 2008) to estimate PET with each of the four models for the period 1895-2020, along with monthly total precipitation (P) for water balance estimates of root zone soil moisture. The 4 km spatial resolution of the PRISM product makes it particularly well suited for application to forest growth models, where processes occur at much finer scales than typically represented in coarser resolution gridded meteorological products (Bontemps and Bouriaud 2014). Accuracy of the gridded PRISM product generally compares favorably to other gridded products ( For comparison to VS-Lite's leaky bucket soil moisture simulations, I also obtained daily (at 0:00 UTC), 9 km root zone (0-100 cm) soil moisture estimates from the Soil Moisture Active Passive (SMAP) Level 4 product for its operational period April 2015-December 2020, which I averaged to monthly scale (a total of 69 months). The SMAP Level 4 product assimilates satellite L-band passive microwave brightness temperature observations, which are sensitive to moisture content in the upper layers of the soil and vegetation, into a detailed land surface model forced with instrumental precipitation observations . The land surface model used by the SMAP Level 4 product is much more sophisticated than the leaky bucket model used by VS-Lite, and the assimilation of satellitebased radiometric brightness temperatures improves spatial representation of soil moisture compared to a model forced solely with station observations (Reichle et al 2017, Dong et al 2019, especially in regions with sparse meteorological data (Dong et al 2019). The combination of more sophisticated land surface model and satellite data assimilation makes the SMAP soil moisture a good benchmark for the leaky bucket soil moisture estimates driven by the different PET models.

Climate data and model projections
To test the sensitivity of VS-Lite to temperature change when driven by the different PET models, I obtained an ensemble of 1/8 • , statistically downscaled CMIP5 historical model runs  and RCP4.5 andRCP8.5 model projections (2006-2099) of monthly T min , T max , and P from the bias correction and spatial disaggregation (BCSD) dataset (table S1 (available online at stacks.iop.org/ERL/16/ 084017/mmedia)) (Maurer et al 2007). To avoid overrepresenting models with large ensembles, I only used a single ensemble member (r1i1p1) per model (Diffenbaugh et al 2018). I did not use the IPSL models because they included projections where monthly T min slightly exceeded T max , which produced complex numbers when used in equations (6) and (8) at many of the sites. The CMCC-CM model had a similar issue, but only at a small handful of sites in Washington, so I used the real component of the complex number at these few sites.
I extracted climate and soil moisture data from the PRISM, SMAP, and BCSD grid cells nearest to the geographic coordinates of each tree-ring site. From both PRISM and BCSD, I estimated monthly T mean as the average of T min and T max . Since the BCSD archive does not provide projections of humidity, VPD, or dewpoint temperature, I calculated sitespecific empirical adjustment factors from the PRISM data to estimate T d,mean from T min (with the constraint that T d,mean cannot exceed T min ), and then used this relationship to estimate projected change in T d,mean based on projected change in T min .

Model calibration and evaluation
For each tree-ring chronology (figure 2), I calibrated each of the four versions of the VS-Lite model (i.e. one version per PET model) using Bayesian parameter estimation (Tolwinski-Ward et al 2013) based on observed, detrended ring-width indices during the period 1895-1960. To test the ability of the different PET-driven VS-Lite models to represent observed growth variation, the calibrated models were then used to simulate growth over the period 1961-2016 and validated using the coefficient of determination (R 2 ) between the simulated growth and the overlapping tree-ring width observations at each site during this independent validation period. In addition to the site-level validation, I also calculated the median R 2 (with bootstrapped 95% confidence intervals) within each Level I EPA Ecoregion of the conterminous US (figure 2) (Omernik 1987).
To test the sensitivity of VS-Lite-simulated soil moisture, soil moisture stress, and ring width to warming-induced changes in PET, I ran each of the four calibrated VS-Lite models at each site using projected changes in T min , T max , and T d,mean under both the RCP4.5 and RCP8.5 scenarios over the period 1951-2099. To isolate just the PET effect, monthly precipitation and monthly growth responses to temperature (g T ) were held constant at their mean historical (1981-2010) monthly values, preserving their seasonality but not allowing long-term change. Thus, the only influence on simulated soil moisture, soil moisture stress, and growth is temperature-induced change in PET. I then calculated percent change (relative to a 1981-2010 baseline) in PET (∆PET), soil moisture (∆M), and soil moisture stress (∆g M ) for each ITRDB chronology using the CMIP5 model ensemble. I also calculated projected relative ring widths for each ITRDB chronology under all models and scenarios by summing the monthly growth index, g (equation (3)), over a fixed growing season each year (from the previous September through December in the year of growth) and then dividing the annual simulated radial growth index by its 1981-2010 mean growth index.

Model performance during the historical period
Across the US, the ability of VS-Lite to replicate observed variation in growth was not strongly influenced by the choice of PET model ( figure 3; table 1). Among the 1279 tree-ring chronologies, R 2 ranged from roughly 0-0.8 (figures 3(a)-(d)), with median values across the US ranging from 0.14 for the Th model to 0.16 for the PM model but with overlapping 95% confidence intervals indicating no significant difference in model performance across the US (table 1). This similarity in performance likely resulted from two sources: similar seasonal and interannual variabilities in the simulated soil moisture regardless of PET model (figures S1 and S2) and adjustments in the M1 and M2 parameters that accounted for the small differences in baseline soil moisture among the PET models (figures S3 and S4). While the Th model simulated higher baseline soil moisture than the other three PET models due to its lower baseline PET, the seasonality and interannual variability of simulated soil moisture were similar across all four PET models (figure S1), and they were nearly identically correlated with monthly SMAP soil moisture ( figure S2, table S2). That similarity resulted in similar M1 and M2 parameters across the four models (figures S3 and S4), though the Th-driven M2 parameter tended to be higher than the others in some regions to adjust for the higher baseline soil moisture simulated with the Th model ( figure S4). VS-Lite model performance also did not differ significantly among the four PET models at the ecoregion scale (table 1): in all nine ecoregions, the median R 2 across the four models did not differ significantly. Regardless of which PET model was used, VS-Lite performed best in dry ecoregions and worst in wet and/or cool ecoregions. All models performed best in Mediterranean California (median R 2 = 0.46; figure 3(k)) and the Temperate Sierras (median R 2 = 0.41; figure 3(m)), with no statistically distinguishable differences among the PET models (table 1). The models performed worst in the Marine West Coast Forests ( figure 3(g)), followed by Northern Forests (figure 3(e)), Northwestern Forested Mountains (figure 3(f)), and Eastern Temperate Forests ( figure 3(h)). However, model performance varied substantially across the Northwestern Forested Mountains, with very poor performance in the northern Rockies and Cascades but relatively good performance in the southern Rockies and Sierra Nevadas (figures 3(a)-(d)).

Model responses to projected warming
As hypothesized, the Th-modeled PET increased the most under projected warming across all ecoregions, followed by Hg-modeled PET, and then PM-and PT-modeled PET (figures 4 and S5). By 2099, the Th model projected roughly 25%-30% increases in PET relative to the 1981-2010 baseline across most ecoregions under RCP8.5, while in most cases the Hg model projected roughly 15%-25% increases and the PM and PT models projected roughly 10%-20% increases ( figure 4). Under RCP4.5, the simulated increases in PET were much lower (roughly 5%-15%)  but with similar differences among the four models ( figure S5). While the Th model projected consistently larger increases in PET than the other three models, they all started from different baseline PET conditions (also shown in Fisher et al 2011) and the seasonality of their projected increases differed, so their projected effects on soil moisture were not as consistent as the changes in PET would suggest (figures S6 and S7). Averaged across all sites and CMIP5 ensemble members, the Th and Hg models projected similarly large declines in soil moisture across all ecoregions (roughly 5%-10% lower than baseline historical conditions by 2100 in RCP8.5) while the PM and PT models generally projected smaller declines under RCP8.5 (figure S6). Magnitudes of simulated soil moisture loss were lower under RCP4.5, but also with larger decreases in the Th and Hg models than in the PT and PM models (figure S7).
These projected PET-driven decreases in soil moisture across all models led to increases in moisture stress (figures S8 and S9) and decreases in simulated radial growth (figures 5 and S10) over the 21st century in all regions, but especially in the driest regions. The projected declines in simulated radial growth resulted from both an extension of the period during which growth was primarily water-limited (rather than temperature-limited) and from increases in the intensity of moisture stress (i.e. a decrease in mean g M during the water-limited part of the year) (figures 6 and S11). However, projected soil moisture stress and radial growth among the different PET-driven models increasingly diverged as warming intensified during the 21st century, especially under RCP8.5 (figures 5 and S8). This divergence among the models was relatively minimal in the wettest and coolest ecoregions (Northern Forests, Marine West Coast Forests, and Eastern Temperate Forests). In the driest ecoregions, however, the Th and Hg models projected increases in soil moisture stress roughly double those of the PT and PM models by 2099 (figures S8(e)-(i)), with the differences among the models increasing towards the end of the century. This translated to large differences among the four PET models in projected growth trends (figure 5), with the Hg model typically projecting the greatest decreases in growth, followed by the Th model, and then the PT and PM models.

Discussion and conclusions
Models of tree ring formation, such as VS-Lite, are often used for attributing growth variation to climatic drivers (Lavergne et al 2015, Mina et al 2016, Chen et al 2017, Tumajer et al 2017, Björklund et al 2019, for projecting growth under climate change (Pompa-García et al 2017, Zeng et al 2019, and for reconstructing past climate , Tolwinski-Ward et al 2014. These models simulate temporal variability in ring width based on knowledge of environmental limitations to cell division and expansion in the vascular cambium, including limitations arising from non-optimal temperatures and inadequate soil moisture. For the latter, models must rely on water balance estimates of soil moisture based on moisture supply (precipitation), demand (evapotranspiration), and infiltration into and runoff from the soil profile. Tree ring models typically use temperature-based PET models to represent atmospheric demand for moisture, but the extent to which the PET model influences the accuracy and trends of ring width simulations, especially when there are strong secular trends in temperature, remains unclear. Here, I examined how four different PET models-varying in complexity and physical realism-affect both the ability of VS-Lite to replicate observed variation of growth and the simulated responses of soil moisture, moisture stress, and radial growth to warming.
Among the US tree-ring chronologies, there were no statistically significant differences in validationperiod R 2 among the models within any ecoregion, and thus no evidence of significant differences in the ability of VS-Lite to accurately simulate interannual variation of ring widths when driven by different PET models as long as the model is calibrated specifically for each PET formulation. These similarities in VS-Lite model skill across all PET models are an encouraging sign for many applications, particularly for data assimilation-based climate reconstructions, for which model inversion benefits from simple model structures with minimal, widely available input data (Tolwinski-Ward et al 2011. However, while much of the similarity in simulated ring widths can likely be traced to similar interannual variabilities in soil moisture regardless of PET method (figure S1), the similarities across different PET formulations could also suggest that VS-Lite's flexible parameterization and semi-empirical (rather than process-based) formulation may make it susceptible to overfitting. Regardless of PET model choice, VS-Lite performed best in dry ecoregions where moisture is strongly limiting to growth, consistent with the principle of limiting factors, one of the integrating theories of dendrochronology (Fritts 1976), and with higher observed sensitivities of tree growth to dry extremes than to wet extremes across much of the US, especially in dry regions . In cooler and more mesic ecoregions of the eastern and northern US, where climatic constraints to growth are less pronounced, the climate-driven VS-Lite model performed much worse. Poor performance in these regions could also partly result from lack of snow accumulation and melt dynamics in the leaky bucket hydrology model.
Despite the similar interannual variabilities in simulated ring widths regardless of PET method, the diverging warming responses of the different PET-driven VS-Lite models suggest caution when using tree-ring models in a nonstationary climate, as secular trends in temperature could lead to spurious trends in PET, moisture stress, and growth. In regions where solar radiation and temperature are the primary limiting factors for growth (e.g. in high latitude or high elevation forests), the influence of PET model choice is likely quite small even in extreme warming scenarios (e.g. figures 5(a), (c) and (d)), as the Liebig's Law formulation of VS-Lite (equation (3)) dictates that growth is limited only by the most limiting factor, which in cold or very wet regions would almost always be temperature or solar radiation, not soil moisture. However, in persistently moisture-limited regions (e.g. figures 5(e)-(i)), the choice of PET model could exert a large influence on growth simulations in a warming climate, with the Th and Hg model simulating much larger increases in soil moisture stress and decreases in growth than the PM or PT models. These differences among the models likely reflect how they interpret temperature change. Both the Th and Hg models include direct responses to temperature, but the PM and PT models respond more indirectly via effects of warming on net radiation (R n ) and the slope of the temperaturesaturation vapor pressure curve (∆), as well as VPD in the PM model (which may partly explain the small differences between the PT-and PM-modeled warming responses). However, the modeled responses of R n , ∆, and VPD to warming are smaller than the direct temperature effects included in the Th and Hg models.
For climate reconstructions that rely on inverting radial growth models, it is possible that the Th PET model could potentially result in significant over-or under-estimation of the target climate fields for periods when temperatures are systematically warmer or cooler than the calibration period, as moisture stress in the Th-based forward model simulations in dry regions may be oversensitive to temperature change. However, since these approaches depend on having a simple and easily invertible model of growth, this is an understandable and perhaps necessary tradeoff. However, using temperature-driven PET methods in radial growth models to project potential changes in forest growth due to climate change could be especially problematic. Both the Th and Hg models are very sensitive to temperature change, and this sensitivity may not reflect the actual physical mechanisms driving evapotranspiration, which is driven by both thermodynamic (net radiation) and aerodynamic (vapor concentration gradient and resistances to vapor transport) processes (Campbell andNorman 1998, Bonan 2016). In a rapidly warming climate, the empirical, temperature-driven PET models will thus likely overestimate change in PET, and thus overestimate loss of soil moisture and moisture stress, leading to potentially spurious changes in model-simulated radial growth. Most climate model ensembles (e.g. CMIP5 and CMIP6) include most or all of the necessary variables for more physical representation of PET. While the additional driver datasets required for the physical PET models may add some additional uncertainties (Berg and Sheffield 2018), forward modeling of growth responses to climate change should likely be performed with physically realistic PET models rather than the temperaturedriven empirical PET models, which are likely illsuited (and not originally intended) for this purpose.

Data availability statement
No new data were created or analyzed in this study. MATLAB code for all analyses is available at https://github.com/mpdannenberg/vs-lite-potentialevapotranspiration.