Analyzing Trade-offs, Synergies, and Drivers among Timber Production, Carbon Sequestration, and Water Yield in Pinus Elliotii Forests in Southeastern Usa

Managing Pinus elliotii forests for timber production and/or carbon sequestration is a common management objective, but can negatively affect water yield due to high losses from evapotranspiration. Thus, understanding the trade-offs and potential synergies among multiple ecosystem goods services, as well as the drivers influencing these interactions, can help identify effective forest management practices. We used available data from 377 permanent plots from the USDA Forest Service Forest Inventory and Analysis Program for 2002–2011, and a forest water yield model to quantify provision levels and spatial distribution and patterns of carbon sequestration, timber volume and water yield for Pinus elliotii ecosystems in North Florida, USA. A ranking-classification framework and statistical analyses were used to better understand the interactions among ecosystem services and the effect of biophysical drivers on ecosystem service bundles. Results indicate that increased biomass reduced water yield but this trade-off varied across space. Specific synergies, or acceptable provision levels, among carbon sequestration, timber volume and water yield were identified and mapped. 1410 Additionally, stand age, silvicultural treatment, and site quality significantly affected the provision level of, and interactions among, the three ecosystem goods and services. The framework developed in this study can be used to assess, map, and manage subtropical forests for optimal provision of ecosystem services.


Introduction
Ecosystem services are the benefits humans receive from natural processes and structures of an ecosystem [1,2].However, because of human alteration, many ecosystem services are increasingly threatened [3].As such, there has been a growing interest in the role of conservation, planning, and management of forests to sustain and optimize ecosystem service provision.Key to this approach is the need to quantify bundles of ecosystem services, their interactions, and their biophysical drivers in a spatially explicit manner [4][5][6].
Quantifying and mapping ecosystem services and goods are gaining interest as means to inform conservation and management approaches, including programs to financially compensate owners of biological resources (e.g., forest property owners) for conserving or managing their land to provide sustained levels of ecosystem services (i.e., payment for ecosystem services) [7].For example, simulation models coupled with spatial analyses have been used to quantify the capacity and economic value of a forest ecosystem to regulate water flows [7].Similar studies have also used this approach to analyze and identify specific areas where conservation interventions were most feasible and cost-effective in providing bundles of ecosystem services [5,8].Additionally, mapping of ecosystem services has been used in environmental assessments to identify locations of specific ecosystem structures and their functional contribution to key goods and services [9,10] and for identifying hotspots for such ecosystem provisioning [6].
Forests provide multiple ecosystem services (e.g., carbon sequestration, water yield, and timber), although landowners have historically been compensated solely for timber products.The timber industry is economically important to the US southeast [11,12].Indeed, forest industry in the State of Florida, United States (USA) contributes to its economy by providing more than 64,000 jobs and accounts for US $13 billion in the state's gross domestic product [13].While timber production is the primary forest industry's objective, many studies recognize the role of pine plantations in contributing to carbon sequestration [14].Carbon sequestration is also a key ecosystem service provided by forests and landowners can be incentivized for efforts to conserve, manage, and improve their forest lands [15] in the form of carbon offsets [11].
Water yield (i.e., water made available to aquifers and water bodies) of a system results from the difference between precipitation and evapotranspiration (ET), and is also a key ecosystem service provided by forests [16].Several studies discuss how this ecosystem service can be negatively affected by increased forest stocking and plantation management [17][18][19], which can vary depending on specific forest structures [20].McLaughlin et al. [21] suggested that management schemes that control key structural attributes (e.g., canopy and understory leaf area index, basal area) can help increase water yield from both plantation and natural forests.
Forest structure attributes like leaf area affect forest functions such as water use, tree growth, and primary productivity [22,23].Leaf Area Index (LAI), defined as the ratio of leaf area to ground area for broadleaf plant canopies or projected needle forests [23,24], is a key metric used to estimate forest ecosystem functions and is strongly related to other structural attributes, such as basal area, which also influences carbon storage processes and timber yields [21,25].Thus, basal area and LAI together are metrics for estimating the effects of biophysical and anthropogenic influences, or drivers, on water yield and aboveground biomass (e.g., carbon, timber) in forests.As such, they can be used to map, quantify, and understand interactions among these three ecosystem services [26].
Understanding ecosystem service interactions represents a challenge in land management decision-making [27,28].Specifically, ecosystem services interact when the level of provision of one service directly affects the level of another service [27].The interaction can either be a trade-off or synergy.Trade-off occurs when production of one ecosystem service compromises the provision of other services.Conversely, a synergy occurs when the productions of multiple ecosystem services either increase or are provided at similar levels.Changes in ecosystem services provisioning and their interactions result from influences of anthropogenic (e.g., management activities) and biophysical (e.g., forest structures and disturbances) drivers [27,[29][30][31].As such, analyses of the interactions among key ecosystem services can provide useful information on how drivers or managing for one service can affect the provision of other services [16].
Studies have previously defined ecosystem services trade-offs and identified spatial patterns at the landscape scale among bundles of ecosystem services for several ecosystem types [6,27,32].Such studies, however, are often based on land cover data, ecosystem service proxies, and graphical analyses.This use of land use/land cover (LULC) data is most commonly cross-sectional and dependent on satellite imagery acquired at one point in time.Thus, the LULC data does not reflect within class variability as they fail to depict small grain temporal and spatial variability due to factors such as age or disturbance.This is specifically true when broad LULC classes and/or low frequency updates are used, which is often the case in ecosystem service quantification and trade-off studies [32].As such, LULC data have been reported to neglect temporal and spatial changes in ecosystem functions [4] and are a poor fit for actual field data and identifying local or regional trends [33].A recent study has utilized geostatistical analysis to identify clusters (e.g., hot spots) of ecosystem services [25].In this study, once ecosystem service clusters were identified, they were analyzed against management drivers, socio-ecological or regulatory dynamics, and/or spatial concordance.Therefore, georeferenced, measured forest structure data for a specific ecosystem type can in-turn better identify specific areas of trade-offs and synergies and provide for a better understanding of the biotic and abiotic drivers affecting these changes [8,[34][35][36].
To better understand trade-offs and synergies among a bundle of ecosystem services, it is imperative to also analyze the anthropogenic and biophysical drivers behind their ecosystem interactions [26,27], but studies that account for these interactions and influences are still limited [4].The use of Geographic Information System (GIS) [29,37] coupled with available plot-level measurements (e.g., forest inventory data), however, can facilitate the development of models that use local-scale ecosystem structural attributes and ecosystem functional relationships (e.g., allometric biomass-carbon equations and LAI-ET relationships) to better quantify ecosystem services provision in specific areas [26,38].
The aim of this study was to analyze carbon sequestration, timber production, and water yield interactions in slash pine (Pinus elliottii) ecosystems in north Florida, USA.The analysis was based on plot-level forest inventory data, ecological function models, and a ranking-classification framework.Additionally, we statistically analyzed common biophysical and forest management drivers affecting these interactions.The specific objectives of this study were to: (1) quantify plot-level timber, carbon, and water yield provision levels; (2) develop a ranking-classification framework to analyze interactions (i.e., trade-offs and synergies) among these ecosystem services and their spatial distribution; and (3) statistically analyze the effect of forest management and biophysical drivers on individual ecosystem services and their interactions.

Carbon Sequestration Estimation
Carbon sequestration estimates were based on aboveground carbon (C) storage changes according to the USDA Forest Service Forest Inventory and Analysis (FIA) program.For specific details on the allometric biomass equation and C storage methods used in the data, see Hoover et al. [39] and Woudenberg et al. [40].Data were obtained for FIA cycles 8 and 9; cycle 8 inventory refers to 2002 to 2007 data, and cycle 9 refers to the 2009 to 2011 data.Often, C markets focus on offsets or temporal changes in tree C stocks [39]; therefore, our analyzed ecosystem service was aboveground tree gross C sequestration, or changes in plot-level C stocks during measurement cycles 8 and 9.We used Pinus elliottii stands in the northeastern and northwestern Florida FIA survey units and plot identification numbers to match plots sampled in cycle 8 to ones re-measured in cycle 9 to estimate plot-level changes during this analysis period.The result was 377 study plots, which are well-distributed over the entire study area (48,967 km 2 ) (Figure 1).
We used FIA plot-level aboveground tree C sequestration and timber volume data obtained from the Florida Geographic Data Library [41].Additionally, since the plot dataset did not include tree-level C estimates, we obtained tree-level data from the FIA website [42].The tree-level dataset provided information on each individual tree sampled including the status condition (i.e., whether the tree is live, felled, or dead), estimated C values measured from the dry biomass portions of the tree for aboveground components (i.e., stump, bole, saplings) and estimated timber volumes according to the market standard [41].The plot-level data also provided ecosystem service driver metric data such as ownership class, stand age, site productivity class, disturbance, silvicultural treatment regimes, and measured basal area, which were used for subsequent analyses.The acronyms used in the following sections are described in Appendix Table A1.
Carbon stores provided by the FIA database at the tree-level were converted to megagram of C per hectare (Mg Cha −1 ) using the conversion factor TPA_UNADJ (i.e., trees per acre unadjusted) according to the three sampling size units used in an FIA plot [40], which are: macroplots (i.e., 0.10 ha plots that originate at the center of subplots), subplots (i.e., plots of approximately 0.02 ha on which trees 7.50 cm and greater in diameter at breast height are measured), and microplots (i.e., 0.001 ha plots, nested within the subplots, on which trees smaller than 7.49 cm in diameter at breast height are measured).These specific factors were 0.999 in macroplots, 6.01 in subplots, and 74.96 in microplots [41].The C storage values for each plot (y) were therefore the sum of the aboveground C in all trees and was determined using Equation (1): where is the dry biomass in the merchantable bole; , the dry biomass in the tree stump; , the dry biomass in the top of the tree; , dry biomass of saplings; and , the dry biomass of woodland tree species.To determine C sequestration during the analysis period, we used the following Equation ( 2): where CSQTNET was the annual carbon sequestered during the time period, in (Mgha −1 yr −1 ), which can be positive or negative depending on management between cycles.Both carbon storage (CSTG1) in the first cycle and carbon storage (CSTG2) in the second cycle were Mgha −1 .Finally, REMPER was the remeasurement period, the number of years between measurements for remeasured plots during cycles 8 and 9.

Timber Volume Estimation
We used the same tree-level plot data from FIA cycle 9 for our timber volume estimates.Specifically, we used the VOLCFSND variable from the FIA database, which is defined as sound volume in the merchantable stem of trees and applies to trees with Diameters at Breast Height (i.e., DBH; diameter at 1.4 m above the surface) greater than 7.5 cm [40].We converted tree-level values of timber volume to m 3 ha −1 with the TPA_UNADJ conversion factor.We summed the values of all trees per plot, and these data were then matched to the plot table using the PLT_CN identifier in the FIA dataset.

Water Yield Estimation
As previously mentioned, forest basal area and LAI have been documented as key forest structures influencing water yield, with LAI as the best predictor of forest water use [21].Leaf Area Index is however not directly measured in the FIA database.Therefore, we used FIA's plot-level basal area data to predict plot-level LAI using McLaughlin et al.'s relationship between LAI and basal area for Pinus elliottii ecosystems in Florida [21]: where LAI is the leaf area index in m 2 m −2 , and BA is reported basal area (from the FIA dataset) in m 2 ha −1 .Modeled LAI was then used to predict the ratio of ET to precipitation (ET/PPT) and ultimately water yield using the models developed by McLaughlin et al. [21]: where WY is annual water yield (m 3 ha −1 yr −1 ); MAP is mean annual precipitation (m) during our analysis period (2002-2011); and, 10,000 is conversion factor yielding WY in m 3 ha −1 yr −1 .Annual precipitation data were obtained from a grid downloaded from the PRISM Climate Group's Oregon State University's website [43].As expected, coupling the relationships in Equations ( 3) and (4) will propagate uncertainty in each of these relationships when predicting water yield.We note that this uncertainty can be large [21], while also noting the applicability of this model for our study since it was developed for Florida P. elliottii stands and requires only basal area data.

Interactions in Supply of Ecosystem Services and Goods
The Land-Use Conflict Identification Strategy (LUCIS) model developed by Carr and Zwick [44] is regularly used to classify lands based on suitability analysis and to determine preferences or appropriateness for agriculture, conservation, or urban land uses.Since provision levels for several ecosystem service goals (e.g., timber volume, C sequestration, or water yield) are linked to ecosystem structural attributes, this ranking method can be used to prioritize areas where a particular goal (i.e., ecosystem service provision level) is being achieved at the highest level when compared to others.Therefore, we used this ranking method as the conceptual basis for a ranking classification framework for ecosystem service trade-offs.In contrast to the LUCIS model and Ecosystem Service mapping approaches, which use raster analyses and proxies [32], our approach was based on vector (point feature) analysis and measured plot data since our variables represent discreet phenomena and measured structures and modeled processes [45].
Since the data for the three ecosystem service variables present different units and ranges in values, we normalized the values using a scale from 0 to 1 by dividing all the values by the maximum values for each ecosystem service based on Carr and Zwick's [44] approach (Figure 2).Specifically, using ArcGIS (Version 10, ESRI, Redlands, CA, USA), we classified the normalized values using the Natural Break (i.e., Jenks) classification method rather than the Manual, Equal Interval, or Standard Deviation classification methods because the data values of the three services were not normally distributed.The Jenks-natural-breaks method classifies data based on Jenk's optimization procedure [46] where the algorithm arranges the values into different classes, relying on an iterative process where different breaks in the dataset are used to minimize the variance within classes and maximize the variance between classes [47].This method produced three provision level classes for each of the three ecosystem services, which were coded 3, 2, and 1, representing high, medium, and low provision levels, respectively (Figure 3).Since we were working with net estimates, negative values (e.g., negative C-sequestration values from land forest conversion activities) were classified as 1 (low).Moreover, we determined interaction codes (IC; i.e., synergy or trade-off) by combining the level of the individual ecosystem service into an "ecosystem service bundle" using the following formula: where IC is the three digit code defining the type of interaction (a priori defined); CSL is the C sequestration provision level; TVL, the timber volume provision level; and WYL, the water yield provision level (Figure 3).The order, or sequence, of the three digits in the IC of our framework is always C sequestration as the first digit, timber the second digit, and water yield as the third.The output codes for the interactions were IC numbers between 111 and 333 (Figure 3) and represented an ecosystem service bundle for each plot and were further classified into plots where one of the ecosystem service is dominant over the others (i.e., trade-off) and plots where at least two of the ecosystem services are dominant (i.e., synergy; Figure 3).Using ArcGIS we then created maps to display the spatial distribution of each ecosystem service, its provision level and the service bundle at the FIA plot-level to better identify the interactions.Table 1 presents the calculated provision levels of the three ecosystem services used in the classification framework.In our framework we defined a synergy as when two of the ecosystem services being analyzed in the IC are produced at the same or higher (high-high, medium-medium) level than the third one [27,32].More explicitly, synergies occur when two of the services are generated at: (1) the same level 3 (high and high) while the third is at level 2 (medium); (2) when two of the services are generated at the same level 2 (medium and medium) while the other at the level 1 (low); or (3) when the three services are yielded at the same level, either 1 (low), 2 (medium), or 3 (high).We defined a trade-off as when one of the ecosystem services was dominant compared to the other ones [27,32], specifically when: (1) one ecosystem service is generated at the level 3 and the other two at level 1 or 2; or (2) when one of services is produced at the level 2 and the other two at level 1.

Table 1.
Percentage of plots in each category of interactions among carbon sequestration, timber volume and water yield.The order, or sequence, of the three digits in the Code field is always carbon sequestration as the first digit, timber the second digit, and water yield as the third.

Analysis of Drivers
As mentioned in Section 2.1.1,we used FIA plot-level ecosystem service driver data, as defined in this study, from the Florida Geographic Data Library [41], such as ownership class, stand age, site productivity class, disturbance, and silvicultural treatments.We tested the effect of these different driver metrics on the provision of each plot-level ecosystem service, as well as the resulting interactions in the plot-level ecosystem service bundle based on our classification framework codes.We tested forest stand age based on the reported influence of this variable on C and timber [21,25].Stand age was determined from the second inventory period and ranged from 2 to 138 years with 75% of plots being less than 49 years old.
We also analyzed land tenure, site productivity class, disturbance regime and presence of silvicultural treatments as defined by FIA criteria [42].Ownership (i.e., land tenure) was reported as public (e.g., state, federal, national park service, national forest system, etc.) or private forests, but was reclassified into two dummy variables: 0 for public ownership and 1 for private ownership.Specific private land tenure classes (e.g., corporate and non-industrial private forests) were not considered.In addition, we assumed that the plots remained under the same land tenure status over the analysis period.The analyzed plots consisted of approximately 72% under private ownership and 28% as public lands.For the site quality (i.e., productivity class) driver, stands were ranked using FIA criteria as 1, 2, 3, 4, 5, and 6, which indicate stand productivities of: >15.8, between 15.7 and 11.6, 8.4 and 11.5, 6.0 and 8.3, 3.5 and 5.9, and 1.4 and 3.4 (m 3 /ha/yr), respectively [42].Disturbance drivers were based on Cycle 9 data and included presence of plot-level diseases, prescribed or natural fire at the crown or ground level, windstorm, livestock or animal grazing, or any other anthropogenic influence.Plots with disturbance drivers were coded as 1 and undisturbed plots were coded as 0. Approximately 89% of the analyzed plots showed no indication of disturbance.Finally, we tested the influence of silvicultural treatments as drivers of ecosystem service provision [14,20].The FIA plots documented treatments using several stand improvement activities, such as cutting, fertilization, herbicides, or other activities with the objective of enhancing the commercial value of a stand [40].About 78% of the plots were not treated, with only 22% receiving silvicultural treatments at least once during the analysis period.Among the silvicultural treatment used, were cutting (19%), site preparation (0.5%), artificial generation (0.3%), and other unspecified treatments (1.9%).
All the statistical analyses were performed with the JMP Pro (version 10, SAS, Cary, NC, USA) statistical software.We used two different types of statistical tests to evaluate the effect of the drivers (e.g., stand age, basal area, silvicultural treatments, ownership, site quality, and disturbance) on net C sequestration, timber volume, and water yield interactions among these ecosystem services.First, to evaluate the effect of drivers on provisioning rates, we used a multiple regression analysis with a backward selection.Homoscedasticity of the residuals were examined using scatter plots and the visual relationship between residual-predicted values.Normality was examined with the normal quantile plot and the Shapiro-Wilks goodness of fit test.When normality assumptions were not met, data for the three variables were log transformed.A p < 0.05 and parameter estimates were used to determine statistical significance of a particular driver on net C sequestration, timber volume, and water yield.
Second, using the coded ecosystem services provision level as the dependent variables, we used a multiple logistic regression to test the effect of the same drivers (e.g., stand age, basal area, silvicultural treatments, ownership, site quality, and disturbance) on the interactions between net C sequestration and water yield, as well as between timber volume and water yield.Using Timilsina et al.'s [25] approach that analyzed the effects of drivers on forest C stock hotspots in Florida, two dummy variables (i.e., 0 and 1) were analyzed as response variables to represent trade-offs and synergies among services.Specifically, trade-off plots were "0" and the synergistic plots were "1".The results of the multiple logistic regression were interpreted using a chi square of the Likelihood Ratio tests, parameter estimates and/or odds ratios, and all were tested for significance using an α = 0.05.

Ecosystem Services Provision and Interactions
Table 2 provides the mean and range for all 3 ecosystem services across the 377 plots.Table 3 lists the descriptive statistics and categorical levels of the three studied ecosystem services.Overall, the classification framework indicated that for all the variables, the data histograms have a positive skew, indicating that plots provided low levels for all ecosystem services (Figure 2).However, as shown in Figures 4-6, the spatial distribution of high water yield values mostly overlap with low and medium values for net C sequestration and timber volume.Categorically classifying provisioning of ecosystem services allowed interactions to be calculated (i.e., IC categories) and these were used accordingly to evaluate synergies and trade-offs among services as presented in Table 1.Mapping plot-level ecosystem service provision (Figures 4-6) displays the spatial distribution of the three provision levels across the study area, and Figure 7 displays the plots where trade-offs (i.e., crosses) and synergies (i.e., pins) occurred.

Effects of Drivers on Individual Ecosystem Service
Results indicate that stand age, treatment and site quality were significant drivers of C sequestration (Table 4).Carbon sequestration decreased with increases in site age and as a result of implementing silvicultural treatments.Furthermore, site quality was also a significant predictor of C sequestration, as more productive sites (i.e., lower class codes) were associated with higher C sequestration rates.Although ownership was not a statistically significant driver, C sequestration was higher in private forests than in public forests.Finally, disturbance regime was not a statistically significant driver, but disturbed plots tended to have higher C sequestration.
Stand age, silvicultural treatment, and site quality all had a statistically significant effect on timber volume (Table 4).In contrast to C sequestration, timber volume was greater in older stands (p = 0.0455), and stands that received silvicultural treatments had higher merchantable timber volume (p < 0.0001).Overall, in plots with silvicultural activities, merchantable timber volume was more likely to increase, compared to plots with no silvicultural treatments.Site quality also had a significant effect (p < 0.0001) on timber volume, but ownership and disturbance regime did not.
Finally, all drivers, except disturbance, had a significant effect on water yield (Table 4).Overall, water yield significantly decreased as stands aged (p < 0.0001), and greater water yield values were associated with treated plots (p < 0.0001).Ownership was also a significant predictor, as stands managed under private ownership yielded more water, compared to public owned forests.Finally, in contrast to net C sequestration and timber volume, higher water yield values were more associated with lower site quality stands.Our multiple logistic regression also corroborated the significant effect of these same drivers (e.g., stand age, silvicultural treatments, ownership, site quality, and disturbance) on the interactions between net C sequestration and water yield, as well as between timber volume and water yield as shown in Table 5.Our study used available plot-level forest inventory data, ecological function models, and a ranking-classification framework to map and better understand C sequestration, timber production and water yield interactions in Pinus elliottii forests in North Florida, USA.It also used commonly measured forest management and ecological disturbance metrics as means to better understand the biotic and abiotic drivers affecting these ecosystem service interactions.
During the analysis period, stands sequestered between −11.7 and 9.2 (mean = 0.6) Mg ha −1 yr −1 of aboveground tree C.These findings likely underestimate the total value as the belowground portions and the fate of harvested products, which delay C emissions [25], were not accounted for in this study.Using our ranking classification framework, 72% of the plots studied were classified as providing low C sequestration levels, with an average of −0.51 Mg ha −1 yr −1 .Nearly 22% and 6% of the plots were categorized as medium and high provision, with averages of 2.7 and 5.7 Mg ha −1 yr −1 , respectively.While this might be a result of the classification method used, these C sequestration rates are consistent with the findings of Timilsina et al. [25], who suggest that slash pine forests present a lower probability of being a hotspot for C storage compared to other forest types (e.g., upland hardwood and oak hickory).In addition, our study's mean (0.6 Mg ha −1 yr −1 ) and maximum (9.2 Mg ha −1 yr −1 ) values were slightly lower than C sequestration rates (10-12 Mg ha −1 yr −1 ) reported by Shan et al. (2001) [14] for a 17-year slash pine stand.However, the mean age for the plots analyzed in our study was 33 years, and it has been reported that that older forest sequester less C than younger ones [48].
Given the importance of timber and water to the regional economy, our results present a means for managing forest for timber while providing other services.We found that merchantable timber volume averaged approximately 55.6 m 3 ha −1 , and our framework identified 67% of the plots in the low provision level with an average of 22.6 m 3 ha −1 .The classification categorized 28% and 5% of the plots as having medium (101.5 m 3 ha −1 ) and high (243.5 m 3 ha −1 ) provision levels.Since slash pine is considered the dominant softwood species in Florida [13,14], these values are within range of estimates reported by [13] for softwoods in Florida, which averaged 82.5 cubic meter per hectare (m 3 ha −1 ).
During the analysis period there was positive water yield for all analyzed FIA plots, with values ranging from 461.0 to 6298.7 m 3 ha −1 yr −1 .This range of water yield values is similar to findings by McLaughlin et al. [21], who report annual water values ranging from 500 to over 6000 m 3 ha −1 yr −1 .However, water yield did decrease as a result of biomass growth.Increased forest water use (i.e., ET) and resulting reductions in water yield with stand age in forest plantations is well documented [17,49,50].However, potential multiple use management practices (e.g., recreation infrastructure, wildlife habitat development) that directly affect drivers and enhance water use efficiency from forest plantations are often not fully considered and adopted into management [20], representing a potential area to improve synergies among multiple ecosystem services (e.g., timber and water yield).
Although most of the three ecosystem services provision levels were classified as low, only 20% of the plots were categorized in the 111 category (i.e., low synergy).However, synergy or trade-offs can be defined a priori based on defined management objectives or societal values such as specific timber targets or pre-defined CO 2 reduction unit offsets [31].Most trade-off interactions were a result of the dominance of water yield over the other two services (38% of all plots), while timber volume dominated in 13% and net C sequestration dominated in12% of plots over the other two services.However, our framework did identify a few areas (only 1%) where the three ecosystem services in the bundle were in synergy at the intermediate provision level (i.e., 222; Table 1).Finally, in a few cases there was synergy between timber and water (around 4% of plots), between C and water (2% of plots), and between C and timber (over 11% of plots).

Management Implications
An important part of this study was to analyze the influence of drivers on the likelihood of a plot being associated with a trade-off or a synergy.As shown in Table 5, stand age, treatment, and site quality were the most significant drivers of interactions, as well as of services individually.Thus, this information can be useful for identifying areas where biomass management activities can satisfy multiple predetermined objectives (e.g., optimize C sequestration offsets while maintaining targeted water yields).
The results reported in this study suggest that site quality, silvicultural treatments, and stand age are all significant drivers of forest C sequestration and merchantable timber volume.The significant effect of site quality on C sequestration and timber volume is plausible as this driver indicates the increased capacity of a land to grow biomass on an annual basis [25,42].Silvicultural activities that control competition from other vegetation (e.g., herbicide, mechanical and/or fire removal of understory) also benefit tree growth and subsequent tree C sequestration in forest stands [14].The lack of a significant effect of disturbance on timber is in-line with another study that found that timber production forests in north Florida are under low risk of hurricane damage [26], although both studies did not explicitly test the same drivers.
It should be noted, however, that our C sequestration rates were for aboveground tree C, meaning that understory contributions to total stand C sequestration were not considered.The potential effect of forest understory on total C sequestration and other services needs further study.Moreover, given FIA's method for measuring silvicultural treatments, we have no way of discerning the effects of different and disparate silvicultural activities (e.g., fertilization versus thinning), which could have different and even opposite influences on ecosystem service provisioning.
Our study, while not proposing specific management strategies for certain biophysical drivers of water yield, found that stand age, silvicultural treatment, ownership, and site quality were significantly associated with water yield.The results suggested that higher stand age is associated with lower water yield.This can be explained by basal area and LAI, which affect water yield, since both increase as a forest becomes older [17,21,49,51].Our findings also show that slash pine stands managed under private ownership tended to have higher water yield than their counterparts managed under public land tenure.This may be due to increased biomass management goals and less silvicultural activities (e.g., clearcuts) on public lands [14,52], as many forest lands in Florida's public areas are managed for biodiversity and other regulation and cultural ecosystem services while forests under private land tenure in the US southern region are primarily managed for timber production [52][53][54].
A limitation of this study was the classification algorithm used to rank and classify the data.Given the continuous nature of the variables, the classes generated by the natural breaks algorithm may hide underlying trends in the dataset.Further research using our data set in its continuous form in multivariate regression, multi-criteria decision analyses or optimization modelling is warranted.Indeed, knowing the actual value the different segments of society place on each of the services could be used to weight the importance of each service-digit in our IC [52].Additionally, the 377 plots used in the study were distributed across nearly 50,000 km 2 suggesting a relatively small sample size.In addition, modeled water yield estimates are modeled predictions and do not account for actually observed water yield [21] or temporal changes (e.g., seasonal climate fluctuations) that occurred during the analysis period.Finally, another limitation is related to plot size used in this analysis, which makes the contribution of a plot, in terms of increase in water yield at a watershed scale, appear to be very minor.

Conclusions
While our study focused on understanding management influences on ecosystem service interactions at the plot scale, our modeling approach and the availability of nation-wide FIA plot-level data means our framework can be easily scaled-up and used to guide forest management at the state, ecoregion, or biome level.Indeed, the use of FIA forest inventory data for regional and state-level applications is one of the program's key missions [40,42].Regional goals could determine if management should focus on synergies (e.g., Type 111) or maximization of other highly valued services (e.g., C sequestration at the expense of water yield).Accordingly, forest stand prescriptions and landscape level management activities could focus on the most influential drivers to best meet these regional goals [31].Most importantly, using field measurements and site-specific informationas opposed to the use of land cover-based proxies-to better understand how land management influences the relative provisioning among ecosystem services is critical [33], and this study provides a framework to do so.
There has been increased interest in forest C sequestration and timber production because of their influential effects on climate change and in providing energy alternatives and other market and profit opportunities [52,54].At the same time, water yield is a highly valued service by society [52].Thus, the framework developed in this study can be used to assess, map and manage subtropical ecosystems for optimal provision of these three and other services.The framework could also provide the necessary information needed for production functions and metrics related to the economic valuation of trade-offs among these ecosystem services [54].This study also provides a repeatable and simplified approach to identify specific areas where synergies occur among different ecosystems services provided by a forest stand dominated by a single tree species (i.e., pine plantations).Our simplified approach can facilitate the quantification of not only ecosystem services but also their different provision levels, management drivers, and locations of synergistic interactions, providing an important tool to better guide multiple-use forest management objectives.

Figure 1 .
Figure 1.The northeastern and northwestern Forest Inventory Analysis (FIA) survey units and study area in the northern portion of the State of Florida, USA.

Figure 2 .
Figure 2. Histograms of the normalized values on a 0-1 scale for net carbon sequestration (a), timber volume (b), and water yield (c) using a 3-class classification scheme generated by the natural breaks algorithm using ARCGIS.Negative values were assumed to be 0. Note that the Y axes are counts and X axes are normalized values.Insert tables indicate the minimum (Min) and maximum (Max) values used to assign each ecosystem service provision level where 1 denotes low, 2 denotes medium, and 3 denotes high according to Jenk's natural breaks algorithm.

Figure 3 .
Figure 3. Diagram of the ecosystem service interaction classification (IC) framework.The order, or sequence, of the three digits in the IC is always carbon sequestration as the first digit, timber the second digit, and water yield as the third.

Figure 4 .
Figure 4. Net carbon sequestration rate provision levels for forest inventory analysis (FIA) Pinus elliotti plots in Florida's northeastern and northwestern survey units, 2002-2011.

Figure 6 .
Figure 6.Water yield provision levels for forest inventory analysis (FIA) Pinus elliotti plots in Florida's northeastern and northwestern survey units, 2002-2011.

Figure 7 .
Figure 7. Ecosystem services interactions among net carbon sequestration, timber volume and water yield for Forest Inventory Analysis (FIA) Pinus elliotti plots in Florida's northeastern and northwestern survey units, 2002-2011.

Table 2 .
Descriptive statistics for three ecosystem services and drivers for Florida Pinus elliotti stands.
Note: N = 377 plots for all variables.

Table 3 .
Descriptive statistics and categorical levels of provision of each ecosystem service for Florida Pinus elliotti stands.

Table 4 .
Parameter estimates of the predictors of net carbon sequestration, timber volume, and water yield in Florida Pinus elliotti stands.

Table 5 .
Effect likelihood ratio tests and parameter estimates for synergy (1) and trade-off (0) interactions between net carbon sequestration timber volume and water yield in Florida Pinus elliotti stands.