Determination of Fertility Rating (FR) in the 3-PG Model for Loblolly Pine Plantations in the Southeastern United States Based on Site Index

Abstract

Soil fertility is an important component of forest ecosystems, yet evaluating soil fertility remains one of the least understood aspects of forest science. We hypothesized that the fertility rating (FR) used in the model 3-PG could be predicted from site index (SI) for loblolly pine in the southeastern US and then developed a method to predict FR from SI to test this hypothesis. Our results indicate that FR values derived from SI when used in 3-PG explain 89% of the variation in loblolly pine yield. The USDA SSURGO dataset contains SI values for loblolly pine for the major soil series in most of the counties in the southeastern US. The potential of using SI from SSURGO data to predict regional productivity of loblolly pine was assessed by comparing SI values from SSURGO with field inventory data in the study sites. When the 3-PG model was used with FR values derived using SI values from SSURGO database to predict loblolly pine productivity across the broader regions, the model provided realistic outputs of loblolly pine productivity. The results of this study show that FR values can be estimated from SI and used in 3-PG to predict loblolly pine productivity in the southeastern US.

1. Introduction

The 3-PG model [1], Physiological Principles Predicting Growth, is a process model that predicts forest productivity based on radiation use efficiency, carbon balance, and partitioning. 3-PG and its variants have been calibrated and tested on many commerically important tree species around the globe. Work with eucalypt (Eucalyptus grandis W. Hill ex Maiden) [2,3], patula pine (Pinus patula Schlechdt. et Cham.) [4], loblolly pine (Pinus taeda L.) [5,6], radiata pine (Pinus radiata D. Don) [7], slash pine (Pinus elliottii Engelm. var. elliottii) [8], and ponderosa pine (Pinus ponderosa Dougl. ex. Laws.) [9] have found that 3-PG can accurately predict growth under a variety of management regimes and climatic and edaphic conditions.

3-PG calculates the amount of photosynthetically active radiation (PAR, ${\varphi }_p$) intercepted by a stand (APAR, ${\varphi }_{pa}$) which is then converted into gross primary productivity (GPP,P${}_G$) using canopy quantum efficiency (${\alpha }_c$), constrained by enviromental factors such as vapor pressure deficit (D), mean temperature (T), soil moisture (${\theta }_s$), frost days, and site nutrient status [2]. Canopy quantum efficiency, ${\alpha }_c$, is calculated as:

$${\alpha }_c=f_T{f}_N{f}_f\phi {\alpha }_{Cx}$$

where ${\alpha }_{Cx}$ is the theoretical maximum canopy quantum efficiency, ${\mathrm{f}}_{\mathrm{T}}$, ${\mathrm{f}}_{\mathrm{N}}$, ${\mathrm{f}}_{\mathrm{f}}$, and φ are temperature, nutrition, frost, and physiology related modifiers, respectively. The physiology related modifier (φ) is defined by the multiplication of an age-dependent modifier (${\mathrm{f}}_{\text{Age}}$) and the most restrictive of the vapor pressure deficit (${\mathrm{f}}_{\mathrm{D}}$) and soil water modifiers (${\mathrm{f}}_{\theta }$); such that $\phi ={\mathrm{f}}_{\text{Age}}$ min{${\mathrm{f}}_{\mathrm{D}}{\mathrm{f}}_{\theta }$}. Under non-limiting conditions, all these modifiers have values of 1 [1,2,10]. Net primary productivity (NPP, P${}_n$) is calculated as a fixed proportion of P${}_{\mathrm{G}}$ that accounts for autotrophic respiration. P${}_n$ is then allocated to aboveground and belowground biomass production.

Fertility rating (FR), is a site-specific variable in 3-PG that relates soil fertility to stand productivity. FR is used in 3-PG as an index of nutrient availability and the effects of FR on ${\mathrm{P}}_{\mathrm{G}}$ are determined by the impact of FR on leaf area index (LAI) which in turn affects APAR, canopy light use efficiency, and canopy conductance. The nutrient modifier (${\mathrm{f}}_{\mathrm{N}}$) in 3-PG is calculated as: ${\mathrm{f}}_{\mathrm{N}}=1-(1-{\mathrm{f}}_{\mathrm{N}0}){(1-\text{FR})}^{{\mathrm{n}}_{\text{fN}}}$, where ${\mathrm{f}}_{\mathrm{N}0}$ is the value of ${\mathrm{f}}_{\mathrm{N}}$ when FR = 0 and $n_{fN}$ determines the shape of the response. In 3-PG, $f_{N0}$ is set at 0.5 [6] and $n_{fN}$ is set to 1 so that the fertility modifier is a linear function of FR [11]. FR also has a direct effect on the partitioning coefficient, m, which is used to calculate partitioning to roots [1]. It is calculated as: $\mathrm{m}={\mathrm{m}}_0+(1-{\mathrm{m}}_0)\text{FR}$, where ${\mathrm{m}}_0$ is set to 0.1 [6]. Therefore when FR = 0, m = 0.1 and when FR = 1, m = 1.

Several studies have attempted to use the relationship between soil properties and tree growth to estimate FR. Sampson et al. (2008) [12] used clay and sand percentage obtained from State Soil Geograhic Survey (STATSGO) database to predict FR. Vega-Nieva et al. (2013) [13] used clay content and Ca, K, and Na in the soil profile to predict FR. Pérez-Cruzado et al. (2011) [14] used fertility limitation, water limitation, oxygen limitation, management limitation, and topographic limitation to calibrate FR. In all three cases 3-PG predictions of stand growth using the estimated FR values closely matched the observed data. However, they were not validated on independent sites. Stape et al. (2004) [3] used a different approach based on the observed growth response following fertilization to predict FR in clonal eucalyptus stands. They first calculated fertilizer response (FER) by subtracting the ratio of biomass growth to initial biomass in fertilized plots from control plots. Sites with the lowest observed FER were given an arbitrary FR value of 0.6 and the sites with FER = 0 were given FR values of 1. This method of FR estimation requires fertilization trials on all the soil types in a region which is both costly and time consuming [10].

Site index (SI) is one of the most common measures of site quality because the height growth of dominant trees is less correlated with stocking than diameter growth and highly correlated with productivity [15]. Dye et al. (2004) [16] used SI to estimate FR. However, they used FR as a tunable parameter to match the 3-PG predicted productivity with observed productivity. The USDA NRCS SSURGO database [17] contains information about SI of major soil series across the southeastern US. SI in the SSURGO dataset is defined as the average height that dominant and codominant trees of a given species attain in a specified number of years, and it is applied to well stocked and even-aged stands.

SI is a realized measure of site quality and is widely used in many traditional growth and yield models as a driver of productivity [15]. SI is a function of soil fertility, soil moisture, and climate. Although there is a strong relation between rainfall and forest productivity, previous research has shown that soil fertility has a stronger impact on productivity of loblolly pine in the southeastern United States [18,19,20,21,22,23,24]. Within the natural range of loblolly pine, rainfall exceeds evapotranspiration [25]. Actual evapotranspiration values for loblolly pine in the southeastern US range from 570 to 1050 mm year${}^{-1}$ while rainfall in this region ranges from 1050 to 1825 mm year${}^{-1}$. Furthermore, results from a number of irrigation and fertilization studies show that the response to fertilization is much greater than irrigation. Studies with loblolly pine in the Sandhills of North Carolina [26] and lower Coastal Plain and Piedmont of Georgia [23,27,28] have shown that fertilization has larger impact on growth than irrigation in loblolly pine stands. In a study of fertilization and irrigation in loblolly pine by Albaugh et al. (1998) [26] on an excessively drained, sandy soil in North Carolina, fertilization increased volume growth by 188% over four years comparison to control plots while irrigation only increased volume growth by 17%. In a similar feritlization and irrigation study in loblolly pine in the lower Coastal Plain of Georgia by Cobb et al. (2008) [28] revealed no significant effects of irrigation on height, basal area, and volume growth after six growing seasons, whereas fertilization significantly increased volume growth by 55% compared to control plots. In addition, recent work with rainfall exclusion by Will et al. (2015) [29] showed that an approximate 30% reduction in throughfall did not affect loblolly pine productivity in Florida and Virginia, whereas fertilization significantly increased productivity compared to nonfertilized treatment. Therefore, for loblolly pine in the southeastern US, we hypothesized that soil fertility is the main driver of productivity and that FR could be predicted from SI.

The objective of this study was to use SI for the a priori prediction of FR that could be used as a fixed rather than a tunable parameter in 3-PG${}_{\text{lob}}$ [6], a variant of 3-PG to simulate growth of loblolly pine across the southeastern US. This study also assessed the potential to use SI from the Soil Survey and Geographic Database (SSURGO) [17] to predict FR for individual soil series in a county.

2. Materials and Methods

2.1. Study Site Description

This study was conducted using data from the control plots of a loblolly pine fertilization trial installed at 21 sites [30] located in the southeastern United States (Figure 1). The fertilization trial was established in the early to late 1990s as an incomplete factorial design of nutrient dose (0 to 269 kg ha${}^{-1}$) and application frequency (0, 1, 2, 4, and 6 years) to evaluate the rates and frequencies of fertilization to optimize growth and fertilizer use efficiency. First generation open pollinated seedlings were planted at all sites at densities ranging from 1215 trees ha${}^{-1}$ to 2141 trees ha${}^{-1}$. Study sites represented six physiographic provinces. Four of the sites were located in the Lower Coastal Plain, one in the Upper Coastal Plain, six in the Piedmont, five in the Western Gulf Coastal Plain, four in the Eastern Gulf Coastal Plain, and one in the Valley and Ridge (Table 1). The majority of soils on the study sites have a udic moisture regime. Soil series present at each site were typical forest soils in each physiographic province.

Figure 1. Location of study sites used to develop and test Fertility Rating (FR) for loblolly pine in 3-PG.

Figure 1. Location of study sites used to develop and test Fertility Rating (FR) for loblolly pine in 3-PG.

During the study period, average monthly temperature in July ranged from 31.7 °C in Brunswick, Virginia to 34 °C in Angelina, Texas. The average monthly minimum temperature in January ranged from $-3$ °C in Brunswick, Virginia to 4.3 °C in Brantley, Georgia. Average annual precipitation ranged from 1165 mm in Brunswick, VA to 1491 mm in Montgomery, Mississippi. Highest and lowest monthly average rainfall were 167 mm in Marengo, Alabama and 62 mm in Marion, Georgia, respectively. Average annual precipitation on all the study sites exceeds the potential evaportranspiration rate reported in the southeastern US [25]. Solar radiation ranged from 7.1 MJ m${}^{-2}$ day${}^{-1}$ in December to 21.2 m${}^{-2}$ day${}^{-1}$ in June. Average annual solar radiation ranged from 14.62 MJ m${}^{-2}$ day${}^{-1}$ in Craven, North Carolina to 15.6 MJ m${}^{-2}$ day${}^{-1}$ in Talbot, Georgia.

This study used data from the control plots on each study site to determine the baseline FR and predict growth. The study sites in Brunswick, Virginia; Craven, North Carolina; and Berkley, South Carolina had four control plots and the rest of the study sites had two control plots each. Plot size ranged from 0.028 to 0.059 ha. On these control plots, all the living trees were measured annually for diameter at breast height and total tree height. Mortality and damages to trees were also recorded in each plot.

Table 1. Location (State, County), physiographic regions, latitude / longitude (decimal degrees), soil series, average annual precipitation, average maximum and minimum temperature (°C), planting density (stems ha${}^{-1}$), and plantation year for study sites used to determine relationship between site index and FR in 3-PG.

Study Site	State, County	Physiographic Province	Lat. / Long.	Soil Series	Precipitation (mm year ${}^{-1}$)	Temperature max/min	Density (stems ha${}^{-1}$)	Plantation Year
180101	South Carolina, Kershaw	Piedmont	34.45/-80.50	Lakeland	1208.5	32.6/-0.5	1445	1997
180301	Georgia, Oglethorpe	Piedmont	33.89/-82.91	Mecklenburg	1237.3	32.9/-0.6	1640	1993
180601	Virginia, Brunswick	Piedmont	36.68/-77.99	Cecil	1165.6	31.7/-3.0	1677	1993
180801	North Carolina, Craven	LCP	35.23/-76.97	Leaf	1235.9	32/0.9	1317	1992
181101	South Carolina, Berkeley	LCP	33.19/-80.19	Lynchburg	1289.5	33.4/1.3	1622	1994
181201	Alabama, Coosa	Piedmont	32.91/-86.38	Louisa	1412.6	32.6/0.3	1457	1996
181502	Georgia, Floyd	Valley & Ridge	34.15/-85.38	Townley	1371.8	32.1/-1.0	1655	1998
181503	Texas, Angelina	WGCP	31.13/-94.46	Kurth	1331.3	34.1/2.5	1267	2000
182201	Georgia, Wilkes	Piedmont	33.81/-82.96	Appling	1240.6	32.9/-0.3	1815	1997
183101	Louisiana, Sabine	WGCP	31.72/-93.56	Sacul	1370.6	33.8/1.2	2141	1993
183102	Louisiana, Vernon	WGCP	31.34/-93.18	Sacul	1489.8	33.8/2.1	1756	1994
183601	Mississippi, Kemper	EGCP	32.70/-88.58	Smithdale	1451.4	33.3/0.6	1632	1996
183901	Alabama, Marengo	EGCP	32.37/-87.84	Savannah	1431.5	33.3/1.2	1413	1998
184201	Georgia, Brantley	LCP	31.34/-81.82	Leon	1310.3	33.2/4.3	1781	1994
184202	Georgia, Brantley	LCP	31.34/-81.83	Leon	1308.2	33.2/4.3	1781	1995
184301	Georgia, Marion	UCP	32.17/-84.63	Troup	1275.9	33.1/2.0	1862	1996
184401	Arkansas, Bradley	WGCP	33.49/-92.13	Savannah	1402.1	33.5/0.8	1247	1996
184501	Alabama, Marengo	EGCP	32.25/-87.55	Brantley	1418.5	33.3/1.4	1415	1996
184801	Texas, Newton	WGCP	30.48/-93.78	Evadale	1496.8	33.7/3.7	1264	1999
185201	North Carolina, Montgomery	Piedmont	35.28/-79.94	Herndon	1212.7	32.5/-1.0	1215	1999
185301	Mississippi, Montgomery	EGCP	32.55/ -89.64	Shubuta	1491.2	33.2 / 0.7	1326	1997

Table 1. Location (State, County), physiographic regions, latitude / longitude (decimal degrees), soil series, average annual precipitation, average maximum and minimum temperature (°C), planting density (stems ha${}^{-1}$), and plantation year for study sites used to determine relationship between site index and FR in 3-PG.

Study Site	State, County	Physiographic Province	Lat. / Long.	Soil Series	Precipitation (mm year ${}^{-1}$)	Temperature max/min	Density (stems ha${}^{-1}$)	Plantation Year
180101	South Carolina, Kershaw	Piedmont	34.45/-80.50	Lakeland	1208.5	32.6/-0.5	1445	1997
180301	Georgia, Oglethorpe	Piedmont	33.89/-82.91	Mecklenburg	1237.3	32.9/-0.6	1640	1993
180601	Virginia, Brunswick	Piedmont	36.68/-77.99	Cecil	1165.6	31.7/-3.0	1677	1993
180801	North Carolina, Craven	LCP	35.23/-76.97	Leaf	1235.9	32/0.9	1317	1992
181101	South Carolina, Berkeley	LCP	33.19/-80.19	Lynchburg	1289.5	33.4/1.3	1622	1994
181201	Alabama, Coosa	Piedmont	32.91/-86.38	Louisa	1412.6	32.6/0.3	1457	1996
181502	Georgia, Floyd	Valley & Ridge	34.15/-85.38	Townley	1371.8	32.1/-1.0	1655	1998
181503	Texas, Angelina	WGCP	31.13/-94.46	Kurth	1331.3	34.1/2.5	1267	2000
182201	Georgia, Wilkes	Piedmont	33.81/-82.96	Appling	1240.6	32.9/-0.3	1815	1997
183101	Louisiana, Sabine	WGCP	31.72/-93.56	Sacul	1370.6	33.8/1.2	2141	1993
183102	Louisiana, Vernon	WGCP	31.34/-93.18	Sacul	1489.8	33.8/2.1	1756	1994
183601	Mississippi, Kemper	EGCP	32.70/-88.58	Smithdale	1451.4	33.3/0.6	1632	1996
183901	Alabama, Marengo	EGCP	32.37/-87.84	Savannah	1431.5	33.3/1.2	1413	1998
184201	Georgia, Brantley	LCP	31.34/-81.82	Leon	1310.3	33.2/4.3	1781	1994
184202	Georgia, Brantley	LCP	31.34/-81.83	Leon	1308.2	33.2/4.3	1781	1995
184301	Georgia, Marion	UCP	32.17/-84.63	Troup	1275.9	33.1/2.0	1862	1996
184401	Arkansas, Bradley	WGCP	33.49/-92.13	Savannah	1402.1	33.5/0.8	1247	1996
184501	Alabama, Marengo	EGCP	32.25/-87.55	Brantley	1418.5	33.3/1.4	1415	1996
184801	Texas, Newton	WGCP	30.48/-93.78	Evadale	1496.8	33.7/3.7	1264	1999
185201	North Carolina, Montgomery	Piedmont	35.28/-79.94	Herndon	1212.7	32.5/-1.0	1215	1999
185301	Mississippi, Montgomery	EGCP	32.55/ -89.64	Shubuta	1491.2	33.2 / 0.7	1326	1997

2.2. 3-PG Parameterization

This study used the parameter set developed for 3-PG by Bryars et al. (2013) [6] with several modifications (Table 2). The projected specific leaf area for mature stands (SLA1) was set to be 4 [31,32] and we set tSLA = 6 [33,34,35,36]. The light extinction coefficient for APAR was set to be 0.69 [37]. This study set the value of ${\alpha }_{Cx}$ at 0.053 which is the point between the values of 0.055 and 0.0485 used by Landsberg et al. (2001) [5] and Bryars et al. (2013) [6], respectively for loblolly pine.

2.3. FR Calculation from Site Index

The goal of this study was to determine if SI could be used to calculate a FR value that could be used as a fixed parameter in 3-PG to accurately predict stand growth. This was done in a 3-step process, because stand biomass production is not linearly related to SI [15] and plant productivity is not linearly related to soil fertility [38,39,40] and thus SI and FR should not be linearly related. First SI was calculated in each plot using the observed data on the tree height. Then the relationship between SI and volume in each plot was determined. Finally, the relationship between SI and stand volume was used to determine the relationship between SI and FR where FR varies from 0 to 1.

2.3.1. Calculation of Site Index

SI in each plot was calculated using a dynamic SI model for loblolly pine [41] with the following form:

$$Y=\frac{26.14+X_0}{1+\frac{1455}{X_0{t}^{1.102}}}$$

(1)

where $X_0$ = 0.5($Y_0-26.14+\sqrt{{(Y_0-26.14)}^2+4\times 1455\times Y_0{\mathrm{t}}_0^{1.107}}$), Y is the predicted height in meters at age t, and Y${}_0$ is the predictor height at t${}_0$ years. The dynamic SI model was derived with the generalized algebraic difference approach [42] using a large data set from permanent plots. This model is base-age invariant and estimates SI from any height age combination [41]. Average height of the tallest 80% of trees on a plot was used to derive dominant height [43]. Dominant height in each plot at age 11 or 12 was used to calculate SI.

2.3.2. Determine Relationship between Stand Volume and Site Index

Most of the plots in the study areas were thinned at age 12, so volume at age 11 was used to derive the relationship between SI and yield. Volume for each tree was calculated using the equation: Volume = 0.1365 + 0.0024 × DBH${}^2$ H [44], where volume is in cubic feet, DBH is in inches, and H, total tree height, is in feet. Values were calculated in English units and converted to metric units. Stand volume was determined by summing the volume of individual trees in the plot and using the plot area to adjust to a per hectare value. The relationship between SI and stand volume in the plots was fit using three different equations: linear, exponential, and sigmoidal with the following forms.

• Linear ($\text{Volume}={\beta }_0+{\beta }_1\text{SI})$
• Exponential (Volume = e${}^{{\alpha }_0+{\alpha }_1\text{SI}}$)
• Sigmoid (Volume = ${\displaystyle \frac{{\gamma }_1}{1+e^{{\gamma }_2+{\gamma }_3\text{SI}}}}$)

where ${\beta }_0$ and ${\beta }_1$ represent the intercept and slope parameter for the linear model; ${\alpha }_0$ and ${\alpha }_1$ represent parameters for the exponential function; and ${\gamma }_1,{\gamma }_2$, and ${\gamma }_3$ represent parameters for the sigmoid function. Nonlinear least square fits of sigmoidal and exponential equations were accomplished by using ‘nls’ procedure in “R version 3.1.1”. To ensure the solution was global and not local, different initial values for the parameters were provided for the fits. R${}^2$ and the leave-one-out cross validation method were used to determine the equation with the best fit using the following selection criteria: Root Means Square Error (RMSE), Mean Absolute Error (MAE), and Predicted Residual Sums of Square (PRESS) statistics.

RMSE was calculated as the standard deviation of the difference between y${}_i$-${\widehat{\mathrm{y}}}_{\mathrm{i},-\mathrm{i}}$, where y${}_i$ is observed volume per hectare and ${\widehat{\mathrm{y}}}_{\mathrm{i},-\mathrm{i}}$ is the predicted value for the volume using the model with all observation in the fitting data except the observed volume on that plot. MAE was calculated as:

$$MAE=\sum _{i=1}^n{\displaystyle \frac{\mid y_i-{\widehat{y}}_{i,-i}\mid }{n}}$$

where n is the number of observations and y${}_i$ and ${\widehat{y}}_{i,-i}$ are as described above.

Similarly, the PRESS statistic was calculated as:

$$PRESS=\sum _{i=1}^n{(y_i-{\widehat{y}}_{i,-i})}^2$$

2.3.3. Calculate Relationship between FR and Site Index

After developing a relationship between stand volume and SI, an expression for FR based on SI was derived. SI of loblolly pine in the southeastern United States can approximately range from 10.7 m to 30.5 m [45,46,47]. Therefore, based on this range of site indices reported for loblolly pine in the southeastern United States, we assumed that SI of 10.7 m corresponded to FR = 0 and SI of 30.5 m corresponded to FR = 1. Then the volume at SI values between the maximum and minimum site indices was used to derive the relationship between SI and FR. We fit a sigmoidal equation: FR = $\frac{{\beta }_0}{1+{\mathrm{e}}^{({\beta }_1+{\beta }_2\text{SI})}}$, where ${\beta }_0$, ${\beta }_1$, and ${\beta }_2$ are coefficients derived from data. We used nonlinear least square method to fit the sigmoidal equation using ‘nls’ procedure in “R version 3.1.1”. To ensure the solution is global and not local, different initial values for the parameters were provided for the fits.

2.4. Validation of the Relationship between FR and Site Index

The FR values derived using the above procedure were used as input in 3-${\text{PG}}_{\text{lob}}$ [6] to predict growth in the control plots at the 21 study sites. Weather data including mean monthly temperature, rainfall, frost days, and vapor pressure deficit were obtained from Daymet Surface Weather and Climatic Summaries (http://daymet.ornl.gov). Stand initialization data were obtained from data collected in each plot. Soil texture and moisture information were obtained from Soil Survey Geographic Database (SSURGO). Species specific parameters in Table 2 were used. These data and the FR value calculated from measured SI were input into 3-${\text{PG}}_{\text{lob}}$ and the output was generated for total aboveground biomass, stem number, and LAI annually for each plot. The predicted values from 3-${\text{PG}}_{\text{lob}}$ were compared with measured values in each plot. A linear model was fitted between observed aboveground biomass and model simulated aboveground biomass and the null hypothesis of slope is equal to 1 was tested. Aboveground biomass for each tree was calculated using the equation: Aboveground biomass = 0.026 DBH${}^{2.015}$ HT${}^{0.864}$, where aboveground biomass is in kg, DBH is in cm, and HT, total tree height, is in m [48]. Stand aboveground biomass was determined by summing the aboveground biomass of individual trees in the plot and using the plot area to adjust to a per hectare value. Similarly, a linear model was fitted between observed stem density and fitted stem density and the null hypothesis of slope is equal to 1 was tested using the predicted value as a regressor [49]. Additionally, for simulated and observed values of aboveground biomass production, model efficiency (EF) [50] was calculated.

Table 2. 3-PG parameters and their values for loblolly pine used in this study.

	Parameters	Meaning	Unit	Value
1	pFS2	Ratio of foliage:stem partitioning at stem diameter = 2 cm		0.40
2	pFS20	Ratio of foliage:stem partitioning at stem diameter = 20 cm		0.25
3	StemConst	Constant in stem mass diameter relationship		0.10
4	StemPower	Power in stem massvdiameter relationship		2.50
5	pRx	Maximum fraction of NPP to roots		0.40
6	pRn	Minimum fraction of NPP to roots		0.20
7	SLA0	Projected specific leaf area at the beginning of plantation	m${}^2$ Kg${}^{-1}$	6.40
8	SLA1	Projected specific leaf area for mature stand	m${}^2$ Kg${}^{-1}$	4.00
9	tSLA	Age at which SLA is mean of SLA0 and SLA1	year	6.00
10	k	Extinction coefficient for APAR by canopy		0.69
11	fullCanAge	Age at full canopy cover	year	4.00
12	MaxIntcptn	Maximum proportion of rainfall intercepted by canopy		0.20
13	LAImaxIntcptn	LAI for maximum rainfall interception		5.00
14	${\alpha }_{Cx}$	Maximum canopy quantum efficiency	molC molPAR${}^{-1}$	0.053
15	MaxCond	Maximum canopy conductance	m s${}^{-1}$	0.006
16	LAIgcx	Canopy LAI for maximum canopy conductance		3.00
17	CoeffCond	Defines stomatal response to VPD	mbar${}^{-1}$	0.02
18	BLcond	Canopy boundary layer conductance	m s${}^{-1}$	0.10
19	wSx1000	Maximum stem mass per tree at 1000 trees ha${}^{-1}$	kg tree${}^{-1}$	235.00
20	thinPower	Power in self thinning law		1.60
21	mF	Fraction of mean foliage biomass per tree on dying trees		0.00
22	mR	Fraction of mean root biomass per tree on dying trees		0.20
23	mS	Fraction of mean stem biomass per tree on dying trees		0.40
24	fracBB0	Branch and bark fraction at stand age 0		0.40
25	fracBB1	Branch and bark fraction for mature stand		0.10
26	tBB	Age at which brak fraction is mean of fracBB0 and fracBB1		15.00
27	gammaFx	Maximum litterfall rate	month${}^{-1}$	0.042
28	gammaF0	Litterfall rate at age 0	month${}^{-1}$	0.001
29	tgammaF	Age at which litterfall rate is mean of gammaFx and gammaF0		18.00
30	Rttover	Average monthly root turnover rate		0.0168
31	m0	Value of m when FR is zero		0.10
32	fN0	Value of fN when FR is zero		0.50
33	Tmin	Minimum temperature for growth	${}^{°}$C	4
34	Topt	Optimum temperature for growth	${}^{°}$C	25
35	Tmax	Maximum temperature for growth	${}^{°}$C	38
36	kF	Number of days production lost for each frost day		1
37	MaxAge	Maximum stand age used to compute relative age		40
38	nAge	Power of relative age in age modifier		3
39	rAge	relative age to make age modifier 0.5		0.20
40	y	NPP to GPP ratio		0.47

Table 2. 3-PG parameters and their values for loblolly pine used in this study.

	Parameters	Meaning	Unit	Value
1	pFS2	Ratio of foliage:stem partitioning at stem diameter = 2 cm		0.40
2	pFS20	Ratio of foliage:stem partitioning at stem diameter = 20 cm		0.25
3	StemConst	Constant in stem mass diameter relationship		0.10
4	StemPower	Power in stem massvdiameter relationship		2.50
5	pRx	Maximum fraction of NPP to roots		0.40
6	pRn	Minimum fraction of NPP to roots		0.20
7	SLA0	Projected specific leaf area at the beginning of plantation	m${}^2$ Kg${}^{-1}$	6.40
8	SLA1	Projected specific leaf area for mature stand	m${}^2$ Kg${}^{-1}$	4.00
9	tSLA	Age at which SLA is mean of SLA0 and SLA1	year	6.00
10	k	Extinction coefficient for APAR by canopy		0.69
11	fullCanAge	Age at full canopy cover	year	4.00
12	MaxIntcptn	Maximum proportion of rainfall intercepted by canopy		0.20
13	LAImaxIntcptn	LAI for maximum rainfall interception		5.00
14	${\alpha }_{Cx}$	Maximum canopy quantum efficiency	molC molPAR${}^{-1}$	0.053
15	MaxCond	Maximum canopy conductance	m s${}^{-1}$	0.006
16	LAIgcx	Canopy LAI for maximum canopy conductance		3.00
17	CoeffCond	Defines stomatal response to VPD	mbar${}^{-1}$	0.02
18	BLcond	Canopy boundary layer conductance	m s${}^{-1}$	0.10
19	wSx1000	Maximum stem mass per tree at 1000 trees ha${}^{-1}$	kg tree${}^{-1}$	235.00
20	thinPower	Power in self thinning law		1.60
21	mF	Fraction of mean foliage biomass per tree on dying trees		0.00
22	mR	Fraction of mean root biomass per tree on dying trees		0.20
23	mS	Fraction of mean stem biomass per tree on dying trees		0.40
24	fracBB0	Branch and bark fraction at stand age 0		0.40
25	fracBB1	Branch and bark fraction for mature stand		0.10
26	tBB	Age at which brak fraction is mean of fracBB0 and fracBB1		15.00
27	gammaFx	Maximum litterfall rate	month${}^{-1}$	0.042
28	gammaF0	Litterfall rate at age 0	month${}^{-1}$	0.001
29	tgammaF	Age at which litterfall rate is mean of gammaFx and gammaF0		18.00
30	Rttover	Average monthly root turnover rate		0.0168
31	m0	Value of m when FR is zero		0.10
32	fN0	Value of fN when FR is zero		0.50
33	Tmin	Minimum temperature for growth	${}^{°}$C	4
34	Topt	Optimum temperature for growth	${}^{°}$C	25
35	Tmax	Maximum temperature for growth	${}^{°}$C	38
36	kF	Number of days production lost for each frost day		1
37	MaxAge	Maximum stand age used to compute relative age		40
38	nAge	Power of relative age in age modifier		3
39	rAge	relative age to make age modifier 0.5		0.20
40	y	NPP to GPP ratio		0.47

After the model performance was tested using data from the control plots at the 21 sites, the performance of 3-${\text{PG}}_{\text{lob}}$ was evaluated against independent data not used in the FR model development. We used control plots at six sites of a loblolly pine mid-rotation fertilization trial located in the southern United States. This fertilization trial was established in the mid 1980s as a factorial design of N (0,112, 224, and 336 kg ha${}^{-1}$) and P (0, 28, and 56 kg ha${}^{-1}$) fertilization (see [43,51]). These independent study sites included a wide range of stand ages, soil types, and climate across the southeastern US (Table 3). We used control plots on each study site to independently validate FR derived a priori from SI in the 3-PG model.

Table 3. Location (County, State), latitude / longitude, planting density, range of age, FR, average maximum temperature, average annual precipitation, and soil series of the sites used for independent validation of FR derived from site index.

Location (County, State)	Lat. / Long.	Planting Density (Trees ha${}^{-1}$)	Age (Years)	FR	Avg max T (${}^{°}$C)	Precipitation (mm year${}^{-1}$)	Soil Series
Lancaster, South Carolina	34.55/-80.63	1097	12-16	0.34	36.13	1157.60	Appling
Covington, Alabama	31.20/-86.25	1147	11-19	0.38	35.51	1547.86	Florala
Kemper, Mississippi	32.40/-81.44	865	11-19	0.51	36.76	1429.97	Wilcox
Effingham, Georgia	36.21/-76.94	1509	10-20	0.35	36.28	1208.37	Leefield
Bertie, North Carolina	32.75/-88.45	1442	10-20	0.34	33.88	1280.63	Norfolk
Howard, Arkansas	34.03/-94.02	1000	10-16	0.33	38.20	1375.83	Sacul

Table 3. Location (County, State), latitude / longitude, planting density, range of age, FR, average maximum temperature, average annual precipitation, and soil series of the sites used for independent validation of FR derived from site index.

Location (County, State)	Lat. / Long.	Planting Density (Trees ha${}^{-1}$)	Age (Years)	FR	Avg max T (${}^{°}$C)	Precipitation (mm year${}^{-1}$)	Soil Series
Lancaster, South Carolina	34.55/-80.63	1097	12-16	0.34	36.13	1157.60	Appling
Covington, Alabama	31.20/-86.25	1147	11-19	0.38	35.51	1547.86	Florala
Kemper, Mississippi	32.40/-81.44	865	11-19	0.51	36.76	1429.97	Wilcox
Effingham, Georgia	36.21/-76.94	1509	10-20	0.35	36.28	1208.37	Leefield
Bertie, North Carolina	32.75/-88.45	1442	10-20	0.34	33.88	1280.63	Norfolk
Howard, Arkansas	34.03/-94.02	1000	10-16	0.33	38.20	1375.83	Sacul

2.5. Application of FR Model in SSURGO Database

The wider application of our approach to estimate loblolly pine productivity requires readily available SI values. The SSURGO dataset has detailed data on all soils mapped in the United States including SI values for loblolly pine on the majority of soil series in the southeastern US. We tested the potential to predict loblolly pine productivity across broader regions in the southeastern US using SI derived for individual soil series for three counties in the southeastern US using SSURGO data. SSURGO has SI values of loblolly pine at base age 50 years which were converted to SI base age 25 using a base age invariant SI model [41]. To illustrate the utility of this approach, Kemper County, Mississippi was selected and loblolly pine productivity was estimated spatially in that county with the 3-PG model using the FR values derived from the SSURGO dataset.

3. Results

3.1. Relationship between Site Index and Volume

SI ranged form 16 to 30 m (base age 25) in the control plots. The exponential, sigmoidal, and linear relationship between volume and SI had R${}^2$ values of 0.93 (Figure 2). All parameters used to derive the relationship between site index and volume were highly significant (p < 0.0001). When leave-one-out cross validation was carried out among the three model forms, MAE, RMSE, and PRESS statistics were lowest for the linear model, followed by the sigmoidal model (Table 4). However, the sigmoidal model performed better when extrapolated outside the regressor variable hull. Both the exponential and linear models performed poorly when extrapolated outside the regressor variable hull. For example, the linear model predicted negative yield when site index was below 12.9 m. Therefore, the sigmoidal model was selected as the best model to predict volume (m${}^3$ ha${}^{-1}$) from site index at base age 25 (m):

$$\text{Volume}={\displaystyle \frac{379.57}{1+{\mathrm{e}}^{(4.556-0.185\text{SI})}}}$$

(2)

Figure 2. Scatter plots between stand volume (m${}^3$ ha${}^{-1}$) and site index (m). Left, middle, and right diagram show the exponential, sigmoidal, and linear relationships, respectively.

Figure 2. Scatter plots between stand volume (m${}^3$ ha${}^{-1}$) and site index (m). Left, middle, and right diagram show the exponential, sigmoidal, and linear relationships, respectively.

Table 4. Comparison of model fit statistics for exponential, sigmoid, and linear relationships between stand volume outside bark (m${}^3$ ha${}^{-1}$) and site index (m) at base age 25 based on leave-one-out cross validation, Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Predicted Residual Sums of Squares (PRESS).

Model Name	Model Form	RMSE	MAE	PRESS
Exponential	Volume = e${}^{(2.855+0.094\text{SI})}$	20.60	15.51	13157.4
Sigmoid	Volume = ${\displaystyle \frac{379.57}{1+{\mathrm{e}}^{(4.556-0.185\text{SI})}}}$	19.70	16.36	12029.6
Linear	Volume = -207.693 + 16.115 SI	18.51	15.33	10622.7

Table 4. Comparison of model fit statistics for exponential, sigmoid, and linear relationships between stand volume outside bark (m${}^3$ ha${}^{-1}$) and site index (m) at base age 25 based on leave-one-out cross validation, Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Predicted Residual Sums of Squares (PRESS).

Model Name	Model Form	RMSE	MAE	PRESS
Exponential	Volume = e${}^{(2.855+0.094\text{SI})}$	20.60	15.51	13157.4
Sigmoid	Volume = ${\displaystyle \frac{379.57}{1+{\mathrm{e}}^{(4.556-0.185\text{SI})}}}$	19.70	16.36	12029.6
Linear	Volume = -207.693 + 16.115 SI	18.51	15.33	10622.7

3.2. Relationship between Site Index and FR

Using the sigmoidal relationship to predict volume from SI, values for volume at SI values of 10.7 m and 30.5 m were calculated as 26.7 and 281.8 m${}^3$ ha${}^{-1}$ at age 11. Based on the assumption that these represent the minimum and the maximum SI for loblolly pine in the southeastern United States, the value of FR was set to 0 at 26.7 m${}^3$ ha${}^{-1}$ and the FR value was set to 1 at 281.8 m${}^3$ ha${}^{-1}$. Based on the observation that the relationship between plant productivity and soil nutrient concentration between the nutrient deficient stage and nutrient sufficient stage is linear [18,38,39,40], the difference between 281.8 and 26.7 was distributed evenly between the FR values of 0 and 1 (Table 5) and the relationship between FR and SI was determined using a sigmoidal model (Figure 3):

$$\text{FR}={\displaystyle \frac{1.190}{1+e^{-(-5.899+0.245\text{SI})}}}$$

(3)

Table 5. Value of volume (m${}^3$ ha${}^{-1}$) that corresponded to the value of FR based on the assumed linear incremental relationship and the respective value of site index (m) based on the sigmoidal relationship between stand volume and site index.

Volume (m${}^3$ ha${}^{-1}$)	Site Index (m)	FR
26.7	10.7	0
52.0	14.7	0.1
77.6	17.3	0.2
103.0	19.4	0.3
128.5	21.1	0.4
154.1	22.7	0.5
179.7	24.2	0.6
205.6	25.6	0.7
230.7	27.1	0.8
256.2	28.7	0.9
281.8	30.5	1

Table 5. Value of volume (m${}^3$ ha${}^{-1}$) that corresponded to the value of FR based on the assumed linear incremental relationship and the respective value of site index (m) based on the sigmoidal relationship between stand volume and site index.

Volume (m${}^3$ ha${}^{-1}$)	Site Index (m)	FR
26.7	10.7	0
52.0	14.7	0.1
77.6	17.3	0.2
103.0	19.4	0.3
128.5	21.1	0.4
154.1	22.7	0.5
179.7	24.2	0.6
205.6	25.6	0.7
230.7	27.1	0.8
256.2	28.7	0.9
281.8	30.5	1

Figure 3. Sigmoidal relationship between Fertility Rating (FR) and site index based on the relationship between stand volume and site index illustrated in Table 3.

Figure 3. Sigmoidal relationship between Fertility Rating (FR) and site index based on the relationship between stand volume and site index illustrated in Table 3.

All parameters used to derive the relationship between site index and FR were highly significant (p < 0.0001).

3.3. Model Evaluation

Using the FR values predicted from SI, 3-PG accurately predicted yield and mortality of loblolly pine in the control plots from the 21 sites in this study. Overall, using FR determined a priori from site index using Equation (3), the aboveground biomass derived from DBH and Height (observed aboveground biomass) and predicted aboveground biomass from 3-PG had an R${}^2$ of 0.89 and the measured stand tree density and predicted stand tree density also had an R${}^2$ value of 0.89 (Figure 4). The model efficiency for aboveground biomass was 0.88. The slope of the relationship between observed aboveground biomass and predicted aboveground biomass was not significantly different from 1 (p < 0.001). Similarly, the slope between measured values of stand density and modeled values of stand density was not different from 1 (p < 0.001). Although 3-PG predited yield well at most sites, there was more variation when observed vs predicted growth was examined on individual sites. Figure 5 illustrates results from several study sites where predictions matched well with the observed values and the sites where predictions did not match well with the observed values. Among the sites where prediction did not match well with the observed values were the sites dominated by Spodosols in the Lower Coastal Plain and the sites where some stochastic events caused high level of mortality.

Figure 4. Relationship between observed aboveground biomass (Mg ha${}^{-1}$) and modeled aboveground biomass (left) and measured stem density (stem ha${}^{-1}$) and modeled stem density across the 48 control plots at the 21 study sites located across the southeastern United States. Solid lines represent linear fit between observed and predicted values and dotted lines represent the 1:1 line between observed and predicted values. Group of points around (1400,1000) in left diagram were due to lightning caused mortality.

Figure 4. Relationship between observed aboveground biomass (Mg ha${}^{-1}$) and modeled aboveground biomass (left) and measured stem density (stem ha${}^{-1}$) and modeled stem density across the 48 control plots at the 21 study sites located across the southeastern United States. Solid lines represent linear fit between observed and predicted values and dotted lines represent the 1:1 line between observed and predicted values. Group of points around (1400,1000) in left diagram were due to lightning caused mortality.

Figure 5. The relationship between observed (dotted line, filled rectangles and circles) and 3-PG predicted (solid line, unfilled rectangles and circles) values for aboveground biomass (rectangle) and stand density (circle) in some of the study sites. The species specific parameters in 3-PG to predict biomass and stocking were used from [6]. Planting density (PD, stems ha⁻¹), FR, and location (county, state) for each installation are given inside each plot.

Figure 5. The relationship between observed (dotted line, filled rectangles and circles) and 3-PG predicted (solid line, unfilled rectangles and circles) values for aboveground biomass (rectangle) and stand density (circle) in some of the study sites. The species specific parameters in 3-PG to predict biomass and stocking were used from [6]. Planting density (PD, stems ha⁻¹), FR, and location (county, state) for each installation are given inside each plot.

3-PG also performed well when evaluated against the independent data from the mid-rotation stands not used to develop the relationship between site index and FR. The predicted values from 3-PG explained 73% of the variation in observed aboveground biomass and 86% of the variation in measured stand density in the independent data (Figure 6). However, there was a positive bias in the predicted aboveground biomass in the data set with a significant difference from the 1:1 line. The slope of the relationship between observed stand density and predicted stand density was not significantly different from 1 (p < 0.001). 3-PG predicted stem density explained 88% of the variance in observed stem density. 3-PG predicted mortality accurately on the majority of the sites. The largest discrepancies in mortality predictions were observed in Kershaw, South Carolina; Brantley, Georgia; Vernon, Louisiana; and Sabine, Louisiana. The site in Kershaw, South Carolina suffered higher rates of mortality due to lightning, the site in Brantley Georgia suffered higher mortality due to unknown reasons, and the site in Vernon, Louisiana had no mortality.

The predicted LAI value from 3-PG increased as FR increased, reaching a peak LAI of 5.5, 4, and 2 for FR values of 0.96, 0.64, and 0.21, respectively. LAI measurements were available on 4 of the 21 study sites; 3-PG predicted LAI reasonably well on these sites. About 53% of the variation in observed LAI was described by the predicted LAI (Figure 7). 3-PG also predicted a reasonable pattern for LAI development throught time. Predicted LAI increased as FR increased and the age of maximum LAI was a function of FR. Higher FR values tend to shorten the time to reach maximum LAI (Figure 7).

Figure 6. Relationship between observed aboveground biomass (Mg ha${}^{-1}$) and modeled aboveground biomass (left) and measured stem density (stem ha${}^{-1}$) and modeled stem density in the control plots of mid-rotation sites which are independent from the study sites listed in Table 1. Solid lines represent linear fit between observed and predicted values and dotted lines represent the 1:1 line between observed and predicted values.

Figure 6. Relationship between observed aboveground biomass (Mg ha${}^{-1}$) and modeled aboveground biomass (left) and measured stem density (stem ha${}^{-1}$) and modeled stem density in the control plots of mid-rotation sites which are independent from the study sites listed in Table 1. Solid lines represent linear fit between observed and predicted values and dotted lines represent the 1:1 line between observed and predicted values.

Figure 7. Relationship between measured LAI and 3-PG${}_{\text{lob}}$ predicted LAI (left) and LAI simulations from 3-PG${}_{\text{lob}}$ across fertility gradients (right).

Figure 7. Relationship between measured LAI and 3-PG${}_{\text{lob}}$ predicted LAI (left) and LAI simulations from 3-PG${}_{\text{lob}}$ across fertility gradients (right).

3.4. County Level Productivity Estimation

SI values in SSURGO and observed SI at the study sites were highly correlated. An R${}^2$ value of 0.75 was observed when linear regression was carried out between SI values in this study and SSURGO site index (Figure 8). Table 6 shows the major soil series mapped in Kemper County, Mississippi; Brunswick County, Virginia; and Brantley County, Georgia, the SSURGO estimate of site index, SSURGO site index adjusted for base age 25 using base age invariant SI model developed by Diéguez-Aranda et al. (2005) [41], and the predicted FR values from adjusted SI using Equation (3). We used Kemper County, Mississippi as a case study to evaluate the potential to use SSURGO based data on site index to predict FR for individual soil series, and then used 3-PG to spatially predict productivity of loblolly pine across a county. Figure 9 shows the spatial pattern of aboveground biomass for 12-year-old loblolly pine in Kemper County, Mississippi predicted by 3-PG using SSURGO and SI information.

Figure 8. The relationship between SSURGO site index at base age 25 derived from base age 50 using [41] and the site index in the study sites at base age 25.

Figure 8. The relationship between SSURGO site index at base age 25 derived from base age 50 using [41] and the site index in the study sites at base age 25.

Figure 9. Map of 3-PG predicted aboveground biomass (Mg ha${}^{-1}$) for 12-year-old loblolly pine in Kemper County, Mississippi. The top left diagram shows soil series mapped in Kemper County Mississippi based on SSURGO database. The top right diagram shows site index at base 25 calculated from SSURGO site index at base age 50 using the base age invariant method developed by [41]. The bottom left diagram shows FR values based on site index, and the bottom right diagram shows 3-PG predicted aboveground biomass based on FR values calculated from SSURGO site index.

Figure 9. Map of 3-PG predicted aboveground biomass (Mg ha${}^{-1}$) for 12-year-old loblolly pine in Kemper County, Mississippi. The top left diagram shows soil series mapped in Kemper County Mississippi based on SSURGO database. The top right diagram shows site index at base 25 calculated from SSURGO site index at base age 50 using the base age invariant method developed by [41]. The bottom left diagram shows FR values based on site index, and the bottom right diagram shows 3-PG predicted aboveground biomass based on FR values calculated from SSURGO site index.

Table 6. Soil series, site index at base age 50, site index adjusted for base age 25 using [41], and FR values derived from site index using Equation (3) in Kemper County, Mississippi; Brunswick County, Virginia; and Brantley County, Georgia based on SSURGO database.

		Kemper				Brunswick				Brantley
	Series	SI	Adjusted SI	FR	Series	SI	Adjusted SI	FR	Series	SI	Adjusted SI	FR
	Bibb	30.5	20.7	0.36	Appling	25.9	17.0	0.18	Albany	29.0	19.4	0.29
	Daleville	29.0	19.4	0.29	Ashlar	22.9	14.6	0.11	Bladen	28.7	19.2	0.28
	Freest	27.4	18.2	0.23	Badin	24.4	15.8	0.14	Bonifay	25.9	16.7	0.18
	Jena	30.5	20.7	0.36	Cecil	25.3	16.5	0.16	Centenary	25.9	16.7	0.18
	Kinston	30.5	20.7	0.36	Chewacla	25.6	16.7	0.17	Eulonia	27.4	18.2	0.23
	Kipling	27.4	18.2	0.23	Emporia	25.9	17.0	0.18	Florala	27.4	18.2	0.23
	Kirkville	29.0	19.4	0.29	Enon	20.4	12.8	0.07	Foxworth	24.4	15.8	0.14
	Mantachie	29.9	20.2	0.33	Fluvanna	23.2	14.8	0.11	Fuquay	25.9	17.0	0.18
	Mayhew	27.4	18.2	0.23	Georgeville	24.7	16.0	0.15	Hurricane	27.4	18.2	0.23
	Mooreville	29.0	19.4	0.29	Goldston	22.6	14.4	0.10	Kinston	30.5	20.7	0.36
	Oktibbeha	23.2	14.8	0.11	Helena	25.6	16.7	0.17	Lakeland	22.9	14.6	0.11
	Ora	25.3	16.5	0.16	Herndon	24.4	15.8	0.14	Leefield	25.6	16.7	0.17
	Prentiss	26.8	17.7	0.21	Iredell	20.4	12.8	0.07	Leon	22.9	14.6	0.11
	Quitman	28.0	18.7	0.25	Lignum	23.2	14.8	0.11	Lynn Haven	24.4	15.8	0.14
	Ruston	25.6	16.7	0.17	Madison	22.0	13.9	0.09	Mascotte	21.3	13.4	0.08
	Savannah	24.7	16.0	0.15	Mattaponi	24.4	15.8	0.14	Meggett	30.5	20.7	0.36
	Smithdale	26.2	17.2	0.19	Pacolet	23.8	15.3	0.12	Meldrim	25.9	17.0	0.18
	Sweatman	25.3	16.5	0.16	Rion	24.4	15.8	0.14	Ogeechee	27.4	18.2	0.23
	Wilcox	24.7	16.0	0.15	Riverview	30.5	20.7	0.36	Olustee	24.4	15.8	0.14

Table 6. Soil series, site index at base age 50, site index adjusted for base age 25 using [41], and FR values derived from site index using Equation (3) in Kemper County, Mississippi; Brunswick County, Virginia; and Brantley County, Georgia based on SSURGO database.

		Kemper				Brunswick				Brantley
	Series	SI	Adjusted SI	FR	Series	SI	Adjusted SI	FR	Series	SI	Adjusted SI	FR
	Bibb	30.5	20.7	0.36	Appling	25.9	17.0	0.18	Albany	29.0	19.4	0.29
	Daleville	29.0	19.4	0.29	Ashlar	22.9	14.6	0.11	Bladen	28.7	19.2	0.28
	Freest	27.4	18.2	0.23	Badin	24.4	15.8	0.14	Bonifay	25.9	16.7	0.18
	Jena	30.5	20.7	0.36	Cecil	25.3	16.5	0.16	Centenary	25.9	16.7	0.18
	Kinston	30.5	20.7	0.36	Chewacla	25.6	16.7	0.17	Eulonia	27.4	18.2	0.23
	Kipling	27.4	18.2	0.23	Emporia	25.9	17.0	0.18	Florala	27.4	18.2	0.23
	Kirkville	29.0	19.4	0.29	Enon	20.4	12.8	0.07	Foxworth	24.4	15.8	0.14
	Mantachie	29.9	20.2	0.33	Fluvanna	23.2	14.8	0.11	Fuquay	25.9	17.0	0.18
	Mayhew	27.4	18.2	0.23	Georgeville	24.7	16.0	0.15	Hurricane	27.4	18.2	0.23
	Mooreville	29.0	19.4	0.29	Goldston	22.6	14.4	0.10	Kinston	30.5	20.7	0.36
	Oktibbeha	23.2	14.8	0.11	Helena	25.6	16.7	0.17	Lakeland	22.9	14.6	0.11
	Ora	25.3	16.5	0.16	Herndon	24.4	15.8	0.14	Leefield	25.6	16.7	0.17
	Prentiss	26.8	17.7	0.21	Iredell	20.4	12.8	0.07	Leon	22.9	14.6	0.11
	Quitman	28.0	18.7	0.25	Lignum	23.2	14.8	0.11	Lynn Haven	24.4	15.8	0.14
	Ruston	25.6	16.7	0.17	Madison	22.0	13.9	0.09	Mascotte	21.3	13.4	0.08
	Savannah	24.7	16.0	0.15	Mattaponi	24.4	15.8	0.14	Meggett	30.5	20.7	0.36
	Smithdale	26.2	17.2	0.19	Pacolet	23.8	15.3	0.12	Meldrim	25.9	17.0	0.18
	Sweatman	25.3	16.5	0.16	Rion	24.4	15.8	0.14	Ogeechee	27.4	18.2	0.23
	Wilcox	24.7	16.0	0.15	Riverview	30.5	20.7	0.36	Olustee	24.4	15.8	0.14

4. Discussion

Our hypothesis that FR can be estimated from SI and used as a fixed parameter in 3-PG was supported by the results of these simulations. When the independent estimate of FR from SI was input into 3-PG${}_{\text{lob}}$, the simulated value of yield matched the observed yield well, with an R${}^2$ value of 0.89. The slope of the relationship between observed aboveground biomass and predicted aboveground biomass was not significantly different from 1. Previous work on 3-PG to predict loblolly pine growth and yield by Landsberg et al. (2001) [5] at a single location in North Carolina and by Bryars et al. (2013) [6] on multiple locations across Georgia also showed that 3-PG could predict loblolly pine growth and yield accurately. The present study was carried out on a wider range of sites across the South. This study used 48 control plots located at 21 study sites across 9 states in the southeastern United States and found that aboveground biomass and mortality were predicted reasonably well in the majority of sites in all the physiographic regions. The model also performed well when used against the independent data from mid-rotation stands not used in the model development. The predicted values from the model explained 73% of the variance in the observed aboveground biomass and 86% of the variance in measured stem density, which indicates that the model can be applied across a wide range of ages and stand conditions.

The model to predict FR from site index provided a priori estimates of FR in 3-PG. This is an improvement over previous work with 3-PG that has used FR as a tunable parameter to match 3-PG simulations to observed data [10]. Tuning FR makes detecting the proper value of FR a subjective decision. The validity of 3-PG growth predictions would be improved if a quantitative method could be developed to determine FR based on measured input parameters [6,13]. Determining soil fertility by measuring soil profiles will not always be practical, affordable or accurate for management applications [52]. It is also difficult to quantify soil fertility from soil properties in forested ecosystems because of the complex relationship between soil properties and stand productivity, deep rooting systems of trees that can explore soil deep into the profile, and variability in soil properties and characteristics within relatively small geographic areas [10,16]. 3-PG${}_{\text{lob}}$ is highly sensitive to FR. For example, FR value of 0, 0.5, and 1 for 12-years-old loblolly pine stand at 1677 stems ha${}^{-1}$ in Brunswick, Virginia corresponded to 16.4 Mg ha${}^{-1}$, 90.0 Mg ha${}^{-1}$, and 167.2 Mg ha${}^{-1}$ of aboveground biomass, respectively. FR was also accurately estimated from the SI values in the SSURGO dataset and used in 3-PG to model productivity of loblolly pine across the landscape. The 3-PG model produced realistic values of aboveground productivity when FR values were derived using SI values from the SSURGO dataset. For example in Kemper County, Mississippi 3-PG predicted aboveground productivity ranging from 40-85 Mg ha${}^{-1}$ at 12 years. These values corresponded well with the previously reported range in the southeastern US [57,66,67].

The maximum site index for loblolly pine at base age 25 is approximately 30.5 m [47] and some of the infertile sites in the southeastern United States can have site index values as low as 10.7 m [45,46]. Conceptually, the relationship between soil nutrient concentration and plant growth can be classified into four steps: nutrient deficiency, sufficiency, luxury consumption, and toxicity in plants [38,39,40] and the relationship between soil nutrient supply and plant growth is curvilinear [39,40]. Plant productivity is low when soil nutrient concentrations are low but increases considerably at higher concentrations until nutrients no longer limits growth. After reaching the sufficiency stage, additional nutrient applications result in luxury consumption. The relationship between plant productivity and soil nutrient concentration between the deficiency stage and sufficiency stages is linear [18,38,39,40]. We used a linear function to define the relationship between FR and productivity, and a three parameter sigmoidal function to address three stages of relationship between soil nutrient concentration and plant productivity, namely, deficiency, sufficiency, and luxury consumption.

Mortality has been ignored or considered zero in many previous studies with 3-PG, such as those in clonal eucalyptus stands [2,3,53]. 3-PG predicted mortality poorly in a study with Pinus radiata in southeastern New Zealand [54]. Bryars et al. (2013) [6] also found poor predictions of mortality for loblolly pine in Georgia. In contrast, when used to simulate stocking in native eucalyptus forests in southeastern Australia, 3-PG performed well and explained about 89% of the variability in observed data [55]. Similar results were observed in our study, where the relationship between observed and predicted mortality had an R${}^2$ of 0.88. In the 21 study sites in this study, there were only 4 sites where large discrepancies were observed. Two of these sites suffered higher mortality due to causes not related to stand density. At 10 of 21 study sites 3-PG predicted stocking within ±30 trees ha${}^{-1}$.

Results from this study also demonstrate that 3-PG can be used to accurately reproduce observed LAI in loblolly pine stands. Our model predictions of LAI as a function of FR were reasonable based on observed values [31,56]. Good correspondence was obtained between simulated LAI across fertility gradients and reported LAIs. For example, LAI values of 5.5 on sites with higher FR values matched well with LAI values observed by Akers et al. (2013) [31] and Zhao et al. (2012) [57]. LAI simulated on medium and low fertility ratings corresponded well with LAI reported by Peduzzi et al. (2012) [58] for loblolly pine. Both the maximum LAI and the time to reach it were observed to be functions of FR. Our results support the assertion from Vose et al. (1994) [56] that slow growing stands reach maximum LAI later than fast growing stands. An average value for growth efficiency on control plots for loblolly pine in the southeastern United States is approximately 7.2 m${}^3$ ha${}^{-1}$ yr${}^{-1}$ LAI${}^{-1}$ [26,59]. In this study, the sites with low fertility, such as the control plots in Oglethorpe, Georgia, had current annual increment (CAI) of 17 m${}^3$ ha${}^{-1}$ yr${}^{-1}$ at age 11. The predicted LAI for this site was 2.40, which translates into a growth efficiency of 7.1 m${}^3$ ha${}^{-1}$ yr${}^{-1}$ LAI${}^{-1}$. The sites with high fertility such as the control plots in Marengo, Alabama had a CAI of 34 m${}^3$ ha${}^{-1}$ yr${}^{-1}$ at age 9. The predicted LAI for this site was 5.4, which translates into a growth efficiency of 6.3 m${}^3$ ha${}^{-1}$ yr${}^{-1}$ LAI${}^{-1}$.

Results from this study show that while predictions of loblolly pine aboveground biomass by 3-PG across the broader region were reasonably accurate, predictions at individual specific sites may be poorer. This result has been observed in other models of pine growth such as empirical growth and yield models [15]. In this study 3-PG performed well at most of the sites with finer textured Alfisols and Ultisols. 3-PG predictions of aboveground biomass were less accurate in sandy Spodosols and Entisols. On these soils 3-PG predictions were generally more accurate at younger ages and then declined with time. Temporal variation in nutrient availability may cause the poorer prediction on the sandy soils.

Nutrient availability in pine plantations varies through time, especially for N, which is often the nutrient most limiting productivity [18] and therefore should be highly related to FR. Nitrogen availability is generally high early in the rotation due to the rapid decomposition of organic matter and mineralization of N contained in slash and logging debris [18]. The Assart effect [60,61] causes a relatively short-lived pulse of N during the early phase of stand development as the forest floor and logging debris decompose. The Assart effect is more pronounced in sandy Spodosols [18]. After the Assart effect disappears, N availability declines as the forest floor accumulates and acts as a sink for N [62,63,64]. The impact of the Assart effect and the accumulation of N in the forest floor on N availability is greater on sandy soils with low organic matter that are inherently less fertile than fine-textured soils. Tree growth is initially rapid but then slows as nutrient availability declines and nutrient deficiency develops [18,62]. Albaugh et al. (2006) [65] showed that height growth decreased through time at a sandy site in North Carolina that was similar to the soil at the site in Kershaw, South Carolina where the observed height growth decreased significantly in the later stages of stand development in the control plots. This suggests that FR should vary through time to more accurately predict soil nutrient availability on the sandy soils. Management practices such as N fertilization temporarily increase nutrient availability in soils, which increases growth on nutrient deficient soils [18]. Therefore, including a temporal component to FR may also enable 3-PG to predict the response to silvicultural practices such as N fertilization.

5. Conclusions

In most previous work with 3-PG, FR has been used as a tunable parameter that is adjusted so that predicted values match observed data. An unbiased a priori estimate of FR would greatly enhance the utility of the 3-PG model [6]. The results presented in this study indicated that FR can be estimated a priori from SI and used to accurately predict the growth and stocking of loblolly pine across the southeastern US. This study suggests that FR can also be estimated from the SI values in the SSURGO dataset and used in 3-PG to model productivity of loblolly pine across the landscape. Future refinements that enable FR to vary through time would likely improve 3-PG predictions in loblolly systems on sandy soils of the southeastern United States.