A New Approach for Modeling Volume Response from Mid-Rotation Fertilization of Pinus taeda L. Plantations

Mid-rotation fertilization presents an opportunity to increase the economic return of plantation forests in the southeastern United States (SEUS). For this reason, the Forest Productivity Cooperative established a series of mid-rotation fertilization trials in Pinus taeda L. plantations across the SEUS between 1984 and 1987. These trials identified site-specific responses to nitrogen (N) and phosphorus (P) fertilizers, resulting in increased stand production for 6–10 years after fertilization. There are successful volume response models that allow users to quantify the gain in stand productivity resulting from fertilization. However, all the current models depend on empirical relationships that are not bounded by biological response, meaning that greater fertilizer additions continue to create more volume gains, regardless of physiological limits. To address this shortcoming, we developed a bounded response model that evaluates relative volume response gain to fertilizer addition. Site index and relative spacing are included as model parameters to help provide realistic estimates. The model is useful for evaluating productivity gain in Pinus taeda stands that are fertilized with N and P in mid-rotation.


Introduction
Silvicultural practices are constantly changing in plantation forests [1]. Hence, forest growth and yield systems that are able to predict the effect of different silvicultural practices on stand growth are essential tools needed to guide operational decision-making [2]. Systems with enough flexibility to estimate forest productivity, under a variety of silvicultural practices and management regimes, are available in the literature for common production species [3][4][5].
Soil nutrient availability typically limits forest production [6,7], emphasizing the importance of ameliorating soil nutrient deficiencies in order to enhance forest productivity [8]. To overcome this limitation, the application of N and P for mid-rotation stands has become a regular practice to increase Forests 2020, 11 forest productivity [9]. Hence, volume response models for fertilization are essential to describe forest stand dynamics and to guide decision making in mid-rotation stands [10].
There is a wide range of literature with respect to response models that successfully estimate productivity gains in plantation forests through mid-rotation fertilization [11][12][13]. For instance, [14] proposed the evaluation of mid-rotation fertilization gains using two equations added on top of stand dominant height and basal area projections. The authors indicated that these functions included all fertilization effects concerning stand volume gain in Pinus taeda stands and assumed there was no interaction between dominant height and basal area responses. Ref [15] expanded this framework to include separate effects for different N fertilization levels and P, but they still used a set of decoupled response equations. That is, the effect on growth and stand development for a combination of treatments was the sum of the responses for each separate treatment.
Ref [10] included the concepts of Type 1 and Type 2 responses to silvicultural treatment in growth and yield models for Pinus radiata in New Zealand, in which the author advocated that the yield form of the new response model was better compared to the approach of [14]. Ref [16] developed a stand-level volume response model as a function of the leaf area index (LAI), foliar N, growth efficiency, quadratic mean diameter, and age for Pinus taeda plantations. The author highlighted the accuracy of the predictions and that mid-rotation fertilization responses were well accounted by cumulative volume gains after 4 years across a wide environmental gradient in the SEUS. Ref [17] developed volume response models for Pinus silvestris stands in Finland using mixed modeling and reported the sites with the expected higher gain in stand volume.
However, all modeling efforts lack the fundamental biological response where stand growth gains decrease as nutrient additions increase, which is a property that is well studied in agricultural crops [18,19]. It is reasonable to assume that other factors limit stand productivity after nutritional constrains have been ameliorated [20]. Hence, a bounded response model is required to limit productivity gain as a maximum fertilization level is reached. At the same time, complexities in the nonlinearities present in a growth and yield model make an additive response an unreasonable assumption, particularly with interacting nutritional factors.
This study aimed to develop a post mid-rotation fertilization volume response model for Pinus taeda stands across the SEUS that is biologically sound and approaches an asymptotic response to fertilization when accounting for differences in initial stand conditions.

Study Area and Database Description
The dataset originated from the Regionwide 13 (RW13) trial of the Forest Productivity Cooperative (FPC) Figure 1.
The sites averaged 13 years of age at study establishment. Study sites were established in site prepared, previous loblolly pine plantations (no old field sites) between 1984 and 1987. Sites could have been thinned or unthinned (if thinned at least one year prior to study establishment), had < 1 m 2 of hardwood basal area as competing vegetation, and no mid-rotation fertilization history (lower coastal plain sites could have received P at establishment). No silvicultural treatments were permitted during the study.
At each site, seven fertilization treatments were applied in a completely randomized block design: control (no fertilization), and a factorial combination of three levels of N (112, 224, 336 kg ha −1 ) and two levels of P (0 and 28 kg ha −1 ), as shown in Table 1. At each site, three replications of each treatment were installed. The sites averaged 13 years of age at study establishment. Study sites were established in site prepared, previous loblolly pine plantations (no old field sites) between 1984 and 1987. Sites could have been thinned or unthinned (if thinned at least one year prior to study establishment), had < 1 m 2 of hardwood basal area as competing vegetation, and no mid-rotation fertilization history (lower coastal plain sites could have received P at establishment). No silvicultural treatments were permitted during the study.
At each site, seven fertilization treatments were applied in a completely randomized block design: control (no fertilization), and a factorial combination of three levels of N (112, 224, 336 kg ha -1 ) and two levels of P (0 and 28 kg ha −1 ), as shown in Table 1. At each site, three replications of each treatment were installed. Treatment N (kg ha -1 ) P (kg ha −1 ) 1  0  0  2  112  0  3  112  28  4  224  0  5  224  28  6  336  0  7  336  28 Soils ranged in drainage class from well drained to poorly drained and in family particle class size from very fine to sandy. Subsoil texture included clays and sandy clay loams to silt loams. Table  2 summarizes selected soil and foliar characteristics prior to treatment. Foliage and soil samples were collected in each plot prior to treatments, composited, and analyzed using a CHN analyzer (CE Instruments NC2100 elemental analyzer). Soil phosphorus was determined by Mehlich 3 extraction followed by analysis using inductively coupled plasma-atomic emission spectrophotometry on a Vista-MPX (Varian, Palo Alto, CA). Soil pH was determined with a 1:1 water solution.  Soils ranged in drainage class from well drained to poorly drained and in family particle class size from very fine to sandy. Subsoil texture included clays and sandy clay loams to silt loams. Table 2 summarizes selected soil and foliar characteristics prior to treatment. Foliage and soil samples were collected in each plot prior to treatments, composited, and analyzed using a CHN analyzer (CE Instruments NC2100 elemental analyzer). Soil phosphorus was determined by Mehlich 3 extraction followed by analysis using inductively coupled plasma-atomic emission spectrophotometry on a Vista-MPX (Varian, Palo Alto, CA, USA). Soil pH was determined with a 1:1 water solution. All the sites were measured at establishment, as well as 2, 4, 6, and 10 years after fertilization (years since treatment or YST). Site index (base age 25 years) prior to fertilization was estimated for all the sites using a site index equation developed by [15]. Individual tree and stand measurements Sites included in the study represented a wide range in initial stand characteristics, as shown in Table 3. Overall, the volume response for N fertilization treatments (compared to the control treatment) varied by P application over time, as shown in Figure 2. All the sites were measured at establishment, as well as 2, 4, 6, and 10 years after fertilization (years since treatment or YST). Site index (base age 25 years) prior to fertilization was estimated for all the sites using a site index equation developed by [15]. Individual tree and stand measurements included diameter at 1.37 m aboveground (DBH, cm), total height (ht in m), dominant height (H in m), and the number of live trees (trees per hectare or TPH).
Sites included in the study represented a wide range in initial stand characteristics, as shown in Table 3. Overall, the volume response for N fertilization treatments (compared to the control treatment) varied by P application over time, as shown in Figure 2. Table 3. Summary statistics of the RW13 study at establishment for stand age, number of live trees (TPH), dominant height (H), site index (S), basal area (B), and stand volume (V).   compared to the control treatment over time for treatments where P was not applied (a) and where 28 kg ha −1 of elemental P was applied (b). These data represent average volume response in SEUS.

Fitting Procedure
Parameters for the response model were found iteratively using an optimization algorithm (Nelder-Mead) that maximized the likelihood of the data (Maximum Likelihood framework): where λ is the log-likelihood; f(·) is a model function with expected parameter θ; y i is the vector of independent observations for i = 1...n; x i is the vector of dependent observations for i = 1...n; σ 2 is the expected parameter for the variance; and π is a constant. Least squares is a special case of the Maximum Likelihood (ML), in which the errors are assumed normally distributed, independent, and with constant variance. One of the main advantages of using ML as the fitting procedure is related to its flexibility in allowing for more complex model specifications [4]. The possibility of explicitly modeling the variance (σ 2 ) is needed to further explore the inference around fertilization response and to directly solve for uneven variance along the response period. We used optimization methods to iteratively find parameters to maximize the likelihood (minimize the log-likelihood) using the "optimx" package in R [21].

Data Standardization
Conventionally, treatment response has been calculated as the difference in the target variable for the treated plot minus the control [14]. For volume, the response would be the following: where v xt is the volume response in m 3 ha −1 of treatment x at t years since the treatment was established; V xt is the stand level volume in m 3 ha −1 of treatment x at t years since the treatment was established; and V 1t is the stand level volume in m 3 ha −1 of the control treatment at t years since the fertilizer was applied in the other treatments. The approach works well for experiments followed from establishment, where a carry-over effect in plant reserves might be negligible. However, for treatments applied in established stands at mid-rotation, it is unlikely to have the same starting volume between the control and the fertilization treatments. Stands with different initial sizes are likely to have different amounts of leaf area, and their responses will be correlated to the initial size, potentially obscuring the response of interest. Hence, prior to any evaluation, the dataset was standardized as a relative response with respect to each initial value: where rV is the relative volume of treatment x at t years since the treatment was established; V x1 is the stand-level volume of treatment x at the time when the treatment was established (t = 0); and other variables are previously defined. After calculating the relative volume increase for each treatment (rV), the relative volume response is calculated as the difference with respect to the control treatment: where rv xt is the relative volume response of treatment x at t years since the treatment was established; rV xt is the stand-level relative volume of treatment x at t years since the treatment was established; and rV 1t is the stand-level relative volume of control treatment at t years since the fertilizer was applied in the other treatments.

Model Development
In crop modeling, the Mitscherlich equation is typically applied for quantifying yield response for different fertilization levels: where rv x is the relative volume response of treatment x; N is the level of nitrogen in kg ha −1 applied to treatment x; a is the maximum yield at high levels of N; and b is the nutrient response coefficient. The use of Equation (5) in crop modeling implies that there is a maximum yield (volume response) after a certain level of N applied, as shown in Figure 3.
where is the relative volume response of treatment x; N is the level of nitrogen in kg ha −1 applied to treatment x; a is the maximum yield at high levels of N; and b is the nutrient response coefficient.
The use of Equation (5) in crop modeling implies that there is a maximum yield (volume response) after a certain level of N applied, as shown in Figure 3. However, in forest growth modeling, we typically model over time (t): where is the relative volume response of treatment x at t years since the treatment was established; a is the maximum yield at high t; and b is the time (years, t) response coefficient.
It is assumed through Equation (6) that increases over time until reaching a maximum response (asymptote), regardless of the different levels of N and P applied in the mid-rotation stands, as shown in Figure 4. However, in forest growth modeling, we typically model rv x over time (t): where rv xt is the relative volume response of treatment x at t years since the treatment was established; a is the maximum yield at high t; and b is the time (years, t) response coefficient. It is assumed through Equation (6) that rv xt increases over time until reaching a maximum response (asymptote), regardless of the different levels of N and P applied in the mid-rotation stands, as shown in Figure 4.  Let's consider the decomposition of the a parameter of Equation (6) as a function of the different levels of N and the presence or absence of P:  Let's consider the decomposition of the a parameter of Equation (6) as a function of the different levels of N and the presence or absence of P: where P is applied (P = 1) or not (P = 0); all other variables are previously defined. The asymptote is now constrained by different levels of N and P, which means that the model asymptote changes for different levels of N and their interaction with the application of P is shown in Figure 5.  Let's consider the decomposition of the a parameter of Equation (6) as a function of the different levels of N and the presence or absence of P: where P is applied (P = 1) or not (P = 0); all other variables are previously defined. The asymptote is now constrained by different levels of N and P, which means that the model asymptote changes for different levels of N and their interaction with the application of P is shown in Figure 5.  Hence, the proposed response model is now expressed as: where rv xt is the relative volume response of treatment x at t years since the treatment was established; N is level of nitrogen in kg ha −1 applied in treatment x; P is a binary variable indicating the application of phosphorus in treatment x; and the other variables are previously defined. Equation (8) still displays some shortcomings. The relative volume response is typically higher at sites with lower quality; however, site quality is not represented in the equation. In our dataset, large differences in the initial stand density exist ranging from 551 to 2059 trees ha −1 , as shown in Table 3, which are not represented in the equation.
To overcome the lack of site quality information in the model, the site index prior to any fertilization treatment could be incorporated in the model: where S is the estimated site index in m of each site estimated from the site index equation developed by [15] at the base age of 25 years; other variables are previously defined. The crown ratio is a variable well correlated with canopy size (commonly assessed as the leaf area index) and gives both a measure of the photosynthetic capacity of a tree or stand density. This variable is conventionally used to provide a measure of tree growing space as well as its competition for nutrients, water, and light. It is worth mentioning that while the leaf area index and crown ratio measurements are not easily available in most forest inventories, [22] found that the crown ratio is accurately predicted by relative spacing (Hart-Becking index). Hence, we may assume that the relative spacing should be monotonically related to the leaf area index, which is a variable that is previously reported to be directly related with the response level [16]. We included relative spacing in the model slope: where RS is the relative spacing-RS = 10000 TPH H and H is the dominant height in m; all the other variables were previously defined.
Typically, the linearity of the model errors, independence, and normality of error distribution are assumed [4], although this is rarely the case. We modeled variance using the equation: where βs are the coefficients to be estimated; all the other variables were previously defined. This function enables proper inferences and model extrapolation. This variance function allows the regression assumptions to be met for these data [4]. Variance modeling was completed using a generalized nonlinear regression approach.

Validation
Ref [23] suggested that the simulation techniques displayed the same power for modeling accuracy evaluation as data splitting. In this study, we used 1000 random samples with repetition generated from the non-parametric bootstrapping technique to evaluate modeling accuracy. This technique was employed by [5], who indicated that it worked well, particularly with extreme cases.
Statistical indices bias (T), mean absolute error (MAE), root mean square error (RMSE), model efficiency (EF), and graphical analyses using the R package ggplot2 [24] were used to evaluate the model.
where Pr is the predicted values; O is the observed values; O is the mean of the observed values; and n is the number of observations.

Overall Relative Volume Response to Different Fertilization Levels
The  (17) where all variables are previously defined. All the coefficients are significant at alpha = 0.05 and display proper signs. The fitted Equation (16) implies that the stand-level relative volume response increases with higher levels of N and the relative volume response is higher when N and P are applied together. Lower site qualities, represented here by the site index, also display higher relative volume responses for any fertilization practice compared to sites with higher quality. The relative response increases over time for any treatment, and the relative response increases for stands with higher RS, which means that more growing space (lower competition for nutrients x water x light) is available for the trees.
The statistical indices T = 0.006, MAE = 0.07, RMSE = 0.11, and EF = 0.60 indicated an acceptable level of accuracy for the estimates provided by Equation (16). The observed and estimated stand-level relative volume response were well correlated, as shown in Figure 6, but there is still some random variation that the model was not able to explain.

Overall Relative Volume Response to Different Fertilization Levels
The fitted relative volume response equation and the variance function are where all variables are previously defined. All the coefficients are significant at alpha = 0.05 and display proper signs. The fitted Equation (16) implies that the stand-level relative volume response increases with higher levels of N and the relative volume response is higher when N and P are applied together. Lower site qualities, represented here by the site index, also display higher relative volume responses for any fertilization practice compared to sites with higher quality. The relative response increases over time for any treatment, and the relative response increases for stands with higher RS, which means that more growing space (lower competition for nutrients x water x light) is available for the trees.
The statistical indices T = 0.006, MAE = 0.07, RMSE = 0.11, and EF = 0.60 indicated an acceptable level of accuracy for the estimates provided by Equation (16). The observed and estimated standlevel relative volume response were well correlated, as shown in Figure 6, but there is still some random variation that the model was not able to explain.  Finally, a variance function (Equation (17)) indicated that the error assumptions were not violated, as shown in Figure 7. Forests 2020, 11, x; doi: FOR PEER REVIEW www.mdpi.com/journal/forests Finally, a variance function (Equation (17)) indicated that the error assumptions were not violated, as shown in Figure 7.

Sensitivity Analysis for a Variety of Relative Spacing and Site Index Scenarios over a Variety of N and P Levels
Since the model fit well and the underlying assumptions were not violated, Equation (16) can be used with some confidence for making inferences. Thus, Equation (16)

Sensitivity Analysis for a Variety of Relative Spacing and Site Index Scenarios over a Variety of N and P Levels
Since the model fit well and the underlying assumptions were not violated, Equation (16) can be used with some confidence for making inferences. Thus, Equation (16) was tested for a variety of site index (S) values ( Figure 8) and relative spacing (RS) (Figure 9) scenarios for different levels of N (112, 224, 336, 560, and 896 kg ha −1 ) and P (0 and 28 kg ha −1 ). Finally, a variance function (Equation (17)) indicated that the error assumptions were not violated, as shown in Figure 7.

Sensitivity Analysis for a Variety of Relative Spacing and Site Index Scenarios over a Variety of N and P Levels
Since the model fit well and the underlying assumptions were not violated, Equation (16) can be used with some confidence for making inferences. Thus, Equation (16)      Ref [22] stated a relative spacing bounded between 0.2 and 0.3 as desired in loblolly pine stands. It is worth mentioning at a high level of stocking RS = 0.1 in Figure 9, stands do appear to slightly react to the fertilization treatments, which is due to the limited tree growing space, and this finding agrees with the model hypothesis.
Equation (16) is bounded, which means that the productivity gain of a given site increases with fertilization practices, but there is a limit associated with the fertilization level. These N applications were tested to provide an overview of how the model is bounded, although it is worth mentioning that some of the N applications in Figures 8 and 9 far exceed what was applied in the experiment. Finally, the sensitivity analysis supported the biological soundness of the model approach.  Ref [22] stated a relative spacing bounded between 0.2 and 0.3 as desired in loblolly pine stands. It is worth mentioning at a high level of stocking RS = 0.1 in Figure 9, stands do appear to slightly react to the fertilization treatments, which is due to the limited tree growing space, and this finding agrees with the model hypothesis.
Equation (16) is bounded, which means that the productivity gain of a given site increases with fertilization practices, but there is a limit associated with the fertilization level. These N applications were tested to provide an overview of how the model is bounded, although it is worth mentioning that some of the N applications in Figures 8 and 9 far exceed what was applied in the experiment. Finally, the sensitivity analysis supported the biological soundness of the model approach.

Discussion
We introduced a new relative volume response model in this paper to address mid-rotation fertilization gains in Pinus taeda stands across the SEUS. Relative volume response increases with increasing fertilization levels, but there is a maximum relative volume response after a certain fertilization level. The model produced a biologically sound response curve.
The modeling approach developed in this study provided consistent estimates. Ref [25] discussed how important it is to have models with accurate information, because this information is critical for operational decision-making. Additionally, [2] and [14] demonstrated how useful response model estimates are as an effective tool to address silvicultural questions for a range of site conditions and fertilization levels, especially for mid-rotation Pinus taeda L. stands in southeastern USA (SEUS). These stands are known to be limited by nitrogen (N) and phosphorus (P) availability [26,27]. Hence, the application of N and P at this stage of stand development has become a regular practice to increase forest productivity [9].
Budgetary constraints in the forestry sector indicate how important it is to select the best fertilization regime for a given stand [7,19]. Dramatic growth gains can be obtained with the use of a proper fertilization regime. Consequently, it is essential to have the best analytical tool to allow the decision maker to choose a fertilization regime that maximizes forest growth and profitability [28]. The response model developed in this study will be useful for evaluating the potential productivity gain for Pinus taeda stands that are fertilized with N and P at mid-rotation, because its predictions are accurate and biologically sound. The model structure is also recommended to be applied to

Discussion
We introduced a new relative volume response model in this paper to address mid-rotation fertilization gains in Pinus taeda stands across the SEUS. Relative volume response increases with increasing fertilization levels, but there is a maximum relative volume response after a certain fertilization level. The model produced a biologically sound response curve.
The modeling approach developed in this study provided consistent estimates. Ref [25] discussed how important it is to have models with accurate information, because this information is critical for operational decision-making. Additionally, [2,14] demonstrated how useful response model estimates are as an effective tool to address silvicultural questions for a range of site conditions and fertilization levels, especially for mid-rotation Pinus taeda L. stands in southeastern USA (SEUS). These stands are known to be limited by nitrogen (N) and phosphorus (P) availability [26,27]. Hence, the application of N and P at this stage of stand development has become a regular practice to increase forest productivity [9].
Budgetary constraints in the forestry sector indicate how important it is to select the best fertilization regime for a given stand [7,19]. Dramatic growth gains can be obtained with the use of a proper fertilization regime. Consequently, it is essential to have the best analytical tool to allow the decision maker to choose a fertilization regime that maximizes forest growth and profitability [28]. The response model developed in this study will be useful for evaluating the potential productivity gain for Pinus taeda stands that are fertilized with N and P at mid-rotation, because its predictions are accurate and biologically sound. The model structure is also recommended to be applied to commercial plantations with other tree species, since its properties ensure biologically sound estimates.
Our modeling approach examined the mean effect of fertilization response in stands with a large variation in initial stand density and different site qualities. Some random variation could not be explained by our modeling approach in this study, as previously highlighted. We recommend that future studies include LAI in volume response modeling. LAI combined with RS and site index provide extra information to characterize site conditions and especially stand health. For example, in stands with the same RS and site index, but with different relative volume response for the same level of applied N and P, the difference could be explained by including LAI in the analysis, as demonstrated by [4]. LAI is an important variable that relates directly to the ability of a given stand to photosynthesize [1].
Sites with same site index, LAI and RS respond differently based on soils [19] and limitations other than N and P may exist, and these factors should also be examined in future studies.

Conclusions
We successfully developed a relative volume response model to predict productivity gains due to mid-rotation fertilization in Pinus taeda stands across the SEUS. The model will help select the appropriate site-specific fertilizer level to apply in a given Pinus taeda stand. The model produces biologically sound output for a variety of relative spacing and site index scenarios. This information will help forest managers maximize forest growth and profitability.