Coarse Woody Debris Decomposition Assessment Tool: Model validation and application

Coarse woody debris (CWD) is a significant component of the forest biomass pool; hence a model is warranted to predict CWD decomposition and its role in forest carbon (C) and nutrient cycling under varying management and climatic conditions. A process-based model, CWDDAT (Coarse Woody Debris Decomposition Assessment Tool) was calibrated and validated using data from the FACE (Free Air Carbon Dioxide Enrichment) Wood Decomposition Experiment utilizing pine (Pinus taeda), aspen (Populous tremuloides) and birch (Betula papyrifera) on nine Experimental Forests (EF) covering a range of climate, hydrology, and soil conditions across the continental USA. The model predictions were evaluated against measured FACE log mass loss over 6 years. Four widely applied metrics of model performance demonstrated that the CWDDAT model can accurately predict CWD decomposition. The R2 (squared Pearson’s correlation coefficient) between the simulation and measurement was 0.80 for the model calibration and 0.82 for the model validation (P<0.01). The predicted mean mass loss from all logs was 5.4% lower than the measured mass loss and 1.4% lower than the calculated loss. The model was also used to assess the decomposition of mixed pine-hardwood CWD produced by Hurricane Hugo in 1989 on the Santee Experimental Forest in South Carolina, USA. The simulation reflected rapid CWD decomposition of the forest in this subtropical setting. The predicted dissolved organic carbon (DOC) derived from the CWD decomposition and incorporated into the mineral soil averaged 1.01 g C m-2 y-1 over the 30 years. The main agents for CWD mass loss were fungi (72.0%) and termites (24.5%), the remainder was attributed to a mix of other wood decomposers. These findings demonstrate the applicability of CWDDAT for large-scale assessments of CWD dynamics, and fine-scale considerations regarding the fate of CWD carbon.

A variety of empirical models, including single-exponential [18][19][20], multiple-exponential [21], linear [3], and lag-time [2] have been developed in attempts to characterize CWD decomposition using interactions among wood properties and abiotic factors. However, broad inferences from predictions based on empirical models are constrained by the measurements taken at study sites or regions that were used to develop the predictive equations. In contrast, mechanistic models can integrate fundamental processes using forcing functions based on long-term measurements and experiments, hence providing more robust applicability.
Russell et al. [22] utilized CWD surveys conducted in conjunction with forest inventory plots to effectively estimate CWD decomposition by assessing transitions among decay classes in the eastern United States. Further, several models have been developed to assess CWD decomposition without relying entirely on empirical relationships. Yin [23] developed a numerical model, based on the methodology used to analyze the long-term dynamics of C and nitrogen (N) in forest soils, as suggested by Agren and Bosatta [24]. That model incorporated substrate quality and microbial activity that was sensitive to abiotic conditions. Zell et al. [25] incorporated meta-analysis and a mixed effects model to predict CWD decomposition for 42 species that was sensitive to abiotic conditions. Although these models are not fully mechanistic, they are better than simple empirical models. The process-based soil C model, Yasso, added functionality to assess woody litter decomposition in soils and was tested using litterbag data in Canada [26,27], and subsequently was used as the foundation for assessing CWD decomposition in boreal forests [28].
When considering the role of CWD decomposition with respect to the various forest C pools, environmental fluxes across habitats, and a changing climate, it is necessary to use a mechanistic model that reflects inherent biogeochemical processes. The CWDDAT (Coarse Woody Debris Decomposition Assessment Tool) has been developed to simulate CWD decomposition in forests by incorporating organic matter decomposition processes and ecological drivers regulating CWD mass loss as well as the fate of the wood C for climatic conditions that range from the tropics to the boreal zone [29].
The main functionalities of this model have been tested using different ambient conditions from 89 sites (14˚N to 65˚N latitude and 58 to 139˚W longitude) distributed from the tropics to boreal zones, with a large range in mean annual temperature (-11.8˚to 26.5˚C), annual precipitation (181 to 6,143 mm), annual snowfall (0 to 612 kg m -2 ), and elevation (3 m to 2,824 m above mean sea level). The model has also been tested with respect to wood properties, including species group (hardwood and softwood), wood density (0.1 to 1.0 g cm -3 ) and size (4.0 to 45 cm in diameter). The results from the model assessment showed that CWD decomposition was highly sensitive to climatic factors, elevation, CWD properties and position (standing vs downed) [29]. Accordingly, the assessment affirms that CWDDAT should be a potentially effective tool for estimating CWD dynamics in forests.
Here we report on: (1) model calibration using observations from four sites to affirm whether or not the model performs within the design objectives [29]; (2) model validation using observations from five sites to assess whether or not the model performs as well with different observational inputs as it did under testing and calibration conditions, and (3) model application by assessing decomposition of CWD left by Hurricane Hugo on the Santee Experimental Forest (Santee EF) in South Carolina in 1989. To achieve the first two objectives the observations from nine sites on the FACE (Free Air Carbon Dioxide Enrichment) Wood Decomposition Experiment (FWDE) are used to provide decomposition response data (6 years) for three tree species. The assessment of CWD following Hurricane Hugo is predicated on pre-and post-storm forest assessments which encompass eight decades [30][31][32][33].

Study sites
The FWDE utilized nine Experimental Forests (EF) representing different bio-climatic zones to study CWD decomposition across the continental USA (Fig 1, Table 1) [34]. The study sites range in mean annual precipitation from 576 mm at the Fraser EF in Colorado to 2,259 mm at the Coweeta EF in North Carolina based on local climatic observations from 2000 to 2017 (S1 Table). The lowest annual mean temperature of 2.2˚C is at the Tenderfoot EF in Montana and the highest temperature of 18.6˚C is at the Santee EF site in South Carolina ( Table 1). The elevation ranges from 8 m above mean sea level at Santee EF to 2,710 m at Fraser EF. There are substantial differences in vegetation coverage among sites, ranging from sparse shrub cover at the San Dimas EF in California to intact forests at Coweeta in North Caronia and Caspar Creek EF in California, leaf area index (LAI) ranges from 0.1 to 5.5 m 2 m -2 (Table 1), with an average of 2.6 m 2 m -2 .
The Santee EF is located in South Carolina of USA (33.15˚N, 79.8˚W); it includes a paired watershed system [33], comprising a reference watershed (WS80) and a treatment watershed (WS77). Hurricane Hugo, a Category IV hurricane (1989), severely impacted the Santee EF. Over 80% of trees within WS80 were broken or uprooted, about 130 Mg ha -1 of CWD were left in the watershed, based on the post hurricane inventory [30,35]. All CWD produced by the hurricane on WS80 were left in situ without salvage logging. The Santee EF has biotic and abiotic observations over the last 80 years, including climate, biomass, and hydrology, which are useful for evaluating CWDDAT model applications; the forest is also an area with the highest prevalence of subterranean termites [36]. This work utilizes observations provided by the FWDE [34], which is located on the USDA Forest Service, Experimental Forest and Range Network [37]. Details regarding the FWDE and individual sites are provided by Trettin et al. [34]. Several of the FWDE primary investors (PIs) are U.S. Forest Service employees (Drs. Trettin, Page-Dumroese and Lindner), and they thus are able to request use of the Experimental Forest lands. We do recognize the Experimental Forest managers we have worked with in the acknowledgements. Data supporting the analyses of the post Hurricane Hugo response of CWD on the Santee EF were obtained from the cited literature [30][31][32][33]35].

Observation setup
Three species were used for the FWDE: loblolly pine (Pinus taeda L.) from the Duke Forest FACE site in North Carolina, USA, and aspen (Populus tremuloides Michx.) and birch (Betula papyrifera Marshall) from the Rhinelander FACE site in Wisconsin, USA. Twenty-four aspen and birch logs (1 m in length) and 24 pine logs (2 m in length) were placed horizontally on the surface of the forest floor to represent downed CWD (Fig 2A). Also, 24 loblolly pine logs (2 m in length) were placed vertically to emulate standing dead snags by placing the logs upright on  Table); LAI, leaf area index, m 2 m -2 .

PLOS ONE
Coarse Woody Debris Decomposition Assessment Tool: Model validation and application a PVC pipe in a manner to preclude ground contact ( Fig 2B). The fresh logs following harvesting were deployed in 2011 [34], for this study we utilize response data following 6 years decomposition. Those FACE logs used for this experiment had intact bark.

Sample collection and analysis
A 3 to 4 cm thick disk was cut from the end of each log before they were installed in EFs to measure their initial wood density at time zero (T 0 ). The weight and length of each log and diameter at each end and in the middle of each log were measured, and the dry weight, volume, and surface area for each log were calculated. Initial wood density was measured by water immersion followed by drying (105˚C) and weighing [34]. In the sixth year (T 6 ) following installation, each log was weighed in the field, unless it was too rotten to be handled; two disks (3-4 cm) were also cut from each end of the log, weighed in situ, and immediately placed in a plastic bag and sealed for transport to the laboratory. The T 6 disk samples were used for density determination by immersion in water to determine volume, followed by drying at 105˚C. The moisture content of the disk samples was used to adjust the log field weight to a dry weight basis. The T 6 wood density and calculated log volume were also used to calculate a log weight. The moisture content and wood density averaged from the two disks collected from a FACE log in situ were used to calculate the wood mass loss for the log. The measured log mass loss (MLML) was the difference in log dry weight between the time T 0 and T 6 . The calculated log mass loss (CLML) was based on the difference in log mass based on the log volume and wood density of the disks between the time T 0 and T 6 .

Model calibration and validation
Calibration and validation are important steps in applying a numerical model [38], since (1) model calibration is used to determine whether or not the model is in line with the design and application objectives, and (2) model validation is used to determine whether or not the model is performing as expected based on measurements. Accordingly, a well-performing model will produce similar results in the calibration and validation steps, thereby affirming applicability of the model for the target study. However, the datasets for the model validation should be completely different from those used for calibration. Accordingly, the CWDDAT was calibrated using observations from four sites with a range in climate, elevation and latitude; and the model was validated using observations from the other five sites (Table 1) with varying ecological conditions among the sites. The simulation results were compared to the year 6 th FWDE observations of mass loss from the FACE logs.

Model parameterization
The CWDDAT was parameterized using climate, soil, and vegetation data from the EF sites (Table 1). Climate data were obtained from weather stations near the EFs (S1 Table) and the Daymet database [39] for fill-in missing climatic observations. The T 0 data for each EF site was were summarized for the model parameterization ( Table 2). The logs were divided into different size classes based on the measurements at the time T 0 ; The FACE logs were divided into two species groups: hardwoods-aspen and birch (coded as species 1, Table 2) and softwoodloblolly pine (coded as species 2, Table 2). Assumptions for the simulation were: (1) annual litter fall for each site was initialized based on current forest species and biomass, this is used because litter can cover downed logs thereby affecting moisture content and temperature, (2) LAI (leaf area index) was initialized for each site based on the current forest type because leaf area can affect soil water content and temperature by changes in throughfall and light radiation to the forest floor, and (3) aboveground biomass, litter fall, and LAI varied with time; the increments were estimated by species group, and stand age. The model was also parameterized for assessing CWD decomposition on WS80 at the Santee EF, where Hurricane Hugo left approximately 130 Mg ha -1 in CWD (Table 3) based on the damage inventory conducted by Hook et al. [30]. Data for vegetation, climate, and soil used to parameterize the model for this application are similar to those used for assessing the impact of the hurricane on C sequestration in this watershed [33].

Model performance evaluation
Model calibration and validation were evaluated using four widely employed quantitative methods [40]: the determination coefficient (R 2 , squared Pearson's correlation coefficient), where O i and P i are the observed values and the predicted values obtained from the modelling, respectively; � O is the observed mean. The evaluation criteria are: if E<0.0, the model is not applicable; 0.0�E<0.25, the model performs poorly; 0.25�E<0.5, the model performance is fair; 0.5�E<0.75, model performance is good; and if E�0.75, the model performs excellently.
The evaluation variables, PBIAS (criteria: between -25.0% and 25.0%) and RRS (criteria: 0.0-0.7), are calculated, respectively, as and where SD is the observed standard deviation; RMSE is the root mean squared error between observed and simulated values, calculated as

RMSE ¼
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi P where n is the number of samples or the pairs of the observed and simulated values. The time to achieve 50% mass loss of CWD from the original mass (half-life, T 50 ) is a useful metric for comparing decomposition responses; it was calculated as where k is the decomposition constant, y -1 , based on the single exponential model that is widely used to assess CWD decomposition as follows, where M t is the mass remaining at time t (years); M 0 is the initial mass; k is the CWD decomposition constant. Because the CWD decomposition may not follow a perfect exponential model [3,13,21,43], a function combining power and exponential components (Eq 7) was used to assess the CWD decomposition, i.e., where k 1 and k 2 are the CWD decomposition constants. Since Eq 7 has two decomposition constants, Eq 5 cannot be directly used to compute the half-life of CWD; accordingly, T 50 for this function was calculated using an iterative method [44].

Model calibration
Four metrics for evaluating model performance for the calibration step were used to determine if the model is in the line of the design objectives. The PBIAS resulted from model calibration was 10.3 for MLML, and 7.39 for the CLML (Table 4). RRS was 0.51 and 0.35 for calibration against the MLML and CLML respectively. These metrics were within acceptable ranges [PBIAS 2 (-25.0, 25.0) and RRS 2 (0.0, 0.7)]. The determination coefficient R 2 (Fig 3) showed that the simulated mass loss was highly correlated with both the MLML (R 2 = 0.80) and CLML (R 2 = 0.90). The simulated total mass loss fraction (0.331) over six years (T 0 to T 6 ) for model calibration was about 10.3% lower than the MLML (0.369), and the simulated total mass loss fraction (0.416) over six years was 7.4% lower than the CLML (0.449). The difference in mass loss fraction between the MLML and CLML reflect uncertainties in estimating log mass loss. Model performance efficiency [E = 0.71, E 2 (-1, 1.0)] for the calibration against the measured mass loss showed that the model performance was within the rating range of "Good" (0.50�E<0.75). The efficiency (E = 0.87) of model calibration for the calculated mass loss using the disks collected in situ showed that the model performance was excellent (E�0.75), and reflects a difference in model performance efficiency for model calibration against data obtained from different methods (see Discussion below). These performance metrics from model calibration step indicate that the model is applicable for assessing CWD decomposition across the range conditions on the EFs, and it is in line with design objectives [29].

Model validation
Similar to the model performance evaluation for the model calibration, the four model evaluation metrics provide a basis for assessing the validation step, which affirms the applicability of  (Fig 4; Table 4), respectively. These metrics were within acceptable ranges for the metrics

Application
To further test the model and demonstrate a practical application, we simulated CWD decomposition following the massive biomass blow-down following Hurricane Hugo in 1989. Our simulations tracked the rapid decomposition of the CWD with an estimated 1.6% of the original CWD mass (130 Mg ha -1 ) following the Hurricane Hugo remaining 30 years after the storm. The decomposition pattern followed an exponential function (Fig 5), with a decomposition constant of 0.132 y -1 based on the single exponential model with an intercept forced to 100% of the initial CWD mass (100%), and about 0.139 y -1 based on the exponential model of the simulated decomposition without forcing the intercept. However, the R 2 (squared Person's correlation coefficient) was high (>0.99) for both exponential equations with and without forcing the intercept (Fig 5).
The CWDDAT also accounts for transfer of CWD-C into the forest floor and mineral soil. A small amount of fine woody debris (0.056% of original mass) derived from CWD decomposition remained on the forest floor, and a small amount of C was incorporated into the mineral soil due to leaching. The predicted DOC incorporated into soils over the 30 years was 30.3 g C m -2 in total, accounting for 0.474% of the total CWD resulting from the hurricane. The predicted particulate organic carbon (POC) incorporated into the mineral topsoil was 0.13 g C m -2 in total, only about 0.002% of total C loss to the CWD decomposition over 30 years.
Decomposition of CWD in WS80 was primarily mediated by fungi, accounting for approximately 72.0% of the total CWD mass loss over the last 30 years. Termites, another key decomposer, accounted for approximately 24.5% of the total CWD mass loss. The contributions of other decomposers (including bacteria and beetles) to CWD decomposition were small, accounting for approximately 3.5% of the total mass loss.

Model calibration and validation
The four model evaluation metrics consistently showed that the CWDDAT model is applicable for assessing CWD decomposition in forests across a wide bioclimatic range. However, those metrics varied depending on the basis for comparison; for example, the calibration performance efficiency with the measured CWD mass loss (E = 0.71) was less than the calculated CWD mass loss (E = 0.87). This difference in model performance might be related to the errors between measurement and calculation (Fig 6). The predicted means for model calibration were 10.8 and 8.1% lower than the MLML and CLML, respectively. The simulated mean for validation was 2.1% lower than the measurement, but 3.6% higher than the calculated value. The predicted average wood mass loss from all samples, including those used to calibrate and validate the model, was 5.4% lower than the measured value, and 1.4% lower than the calculated mean, reflecting that the model slightly underestimates CWD decomposition rates.
Since the simulation results are compared against two different estimates of mass loss, the differences in model performance metrics as compared to the "measured" or "calculated" mass loss essentially represent uncertainty in the estimate of mass loss. The important aspect to convey is that the model performance was good, regardless of the source of the mass loss estimate. Fig 6 is important in showing the uncertainties between the "measured" and "calculated" mass loss estimates; it shows that generally mass loss "measured" > "calculated". However, importantly, the model did a fine job predicting decomposition across a wide range of mass loss over the 6 years (differences in mass loss fraction among the logs: <0.02-0.78).
The differences in model performance relative to the basis for comparison reflects both uncertainties within the model algorithms and in the variables used for comparison. The challenge in assessing wood decomposition is to obtain a metric (e.g., mass, wood density) that is representative of the log at each time-step assessed. Since the FWDE is a long-term study, it was not possible to destructively sample entire logs, so weighing the entire log was used to overcome some of the limitations associated with only using samples from the end-of-log disks to represent the log decomposition. However, the assumption that the moisture content from the disk samples represented the log is also problematic, since the disk has one surface https://doi.org/10.1371/journal.pone.0254408.g006 exposed. Accordingly, if the sample moisture content were less than the mean for the entire log, the measured log mass would be overestimated. Analogous issues arise when results from field incubation studies are used to infer decomposition of CWD [2,11,22]. However, it's highlighted in this assessment when to two, albeit not independent, metrics are used as the basis for model comparison.
Reassuringly, both the model calibration and validation showed that the model performed well for assessing CWD decomposition across a wide range in climatic conditions, forest types and soil in the continental USA. However further assessment of the CWDDAT is warranted to affirm that it is robust across all forest lands. The current assessment utilizes relatively small logs, especially in comparison to the large trees typical of the tropics or pacific north west region. Accordingly, trials that focus on large diameter material is suggested. Similarly, it well established that wood decomposition rates vary among species [19,22,[45][46][47], hence testing the model against other species is needed, especially for the dense tropical species. Similarly, testing across the highly variable boreal zone is warranted [5,48,49]. Overall, this work demonstrates that the CWDDAT [29] is a useful tool for assessing CWD decomposition.

Model application
CWD decomposes quickly in warm and humid locations. Accordingly, the fast CWD decomposition at WS80 reflects the subtropical forest in South Carolina, and the presence of subterranean termites. As a result, over 50% of CWD produced by Hurricane Hugo in 1989 had decomposed within 6 years (T 50 � 6). The predicted main contributor of CWD decomposition on Santee is fungi, accounting for about 72% of total CWD mass loss, as the warm and humid climate regime is good for fungal growth [2]. Termites contributed about 24.5% of the total wood mass loss because Santee is located in the area with the highest risk of termites in USA [36].
Although the estimated mean CWD mass loss from the wood decomposition on WS80 followed an exponential decomposition function, there was a small difference in the decomposition constant (k) with and without forcing the intercept of the regressed function to the initial wood mass; forcing the intercept to 100% of the initial mass at T0 seems appropriate to have the function reflect the actual starting condition. The effect of altering the intercept changes the half-life (T 50 ) from 5.0 to 5.2 years based on the decomposition constants from regressed exponential equations with and without forcing the intercept, respectively, reflecting that there is a small time-lag in this subtropical climate regime. However, when Eq 7 was used to assess CWD decomposition, the dual constants k 1 and k 2 were 0.1535 (k 1 ) and 0.1303 (k 2 ), respectively, and T 50 was approximately 6.2 years. Accordingly, the simulated CWD decomposition constant on WS80 was within the global range of CWD decomposition and within the range in North America [25,50]. The CWD decomposition constant from this study was also similar to that reported for other studies from similar climatic conditions [50][51][52].
The temporal contributions of decomposers to the mass loss of the CWD in WS80 were substantially different; the estimated contribution of termites was approximately 72% within the second year, and their contribution decreased rapidly to 28.0% of the total annual decomposition 2 years after the hurricane and their proportion of total decomposition continued to decline after year 2. Similar to termites, the contribution of beetles to the CWD decomposition decreased over time ranging from 5.3% to null. In contrast, the contributions of fungi to decomposition increased with time. Fungal respiration increased from 21.6% of the total annual mass loss within the first year to �70% after one year, and remained over 70% of the total annual wood mass loss until the end of the simulation time period (30 years). Bacterial respiration increased from 0.4% in the first year to 3.5% of the total annual mass loss in 2017, reflected that bacterial contribution to CWD decomposition was small overall [53].
The decayed logs remained the primary pool of CWD (average 98.1 g C m -2 ) after 30 years. The predicted amount of fragmentation during CWD decomposition was approximately 23.5 g m -2 y -1 , with approximately 7.3 g m -2 remaining in the forest floor at the end of 30 years. The large difference in the mass of fragments produced and that remaining on the floor is due to decomposition of the fragmented material within the forest floor. There was only a small amount of wood C as POC (0.13 g C m -2 ) from CWD decomposition, which was incorporated into the soil as a result of leaching.
DOC from CWD decomposition was the primary pathway for C incorporation into the soil [54]. The 30-year predicted total DOC from the CWD produced by Hurricane Hugo was about 35.5 g C m -2 , and the net incorporation was about 30.3 g C m -2 , a rate of approximately 1.01 g C m -2 y -1 on average, reflecting that DOC is the main contributor of wood C into the underlaying mineral soil. DOC incorporation into soil changed non-linearly over time, increasing from 0.16 g C m -2 y -1 in the first year to 5.9 g C m -2 y -1 in the 10 th year after the hurricane. Transfer by DOC decreased substantially after year 10, from 5.9 g C m -2 y -1 in the 10 th year to 0.09 g C m -2 y -1 at the 26 th year due to a reduction in CWD decomposition. The DOC incorporated into soils can be absorbed or aggregated and therefore maintained in soils for a longer time [55,56]. These simulation results provide context for suggesting the magnitude of the transfers of C from CWD into the soil. Those predictions are predicated on functional relationships [29], but those functions warrant testing through field studies.
This application demonstrates the utility of the model for assessing CWD dynamics, in this case following a large episodic event. Similarly, the model could be used to assess the annual turnover of CWD in forest stands affected by different management regimes or climate scenarios. The fine performance of the model across contrasting forest landscapes suggests the opportunity to conduct large-scale assessments of CWD dynamics. The availability of national forest inventory data that includes periodic assessments of CWD stocks (e.g., Forest Inventory and Analysis data in the U.S.) suggests an opportunity to assess the model performance among multiple forest types within or among regions. Such a large-scale assessment would further test the model, and if successful, provide a tool that could be used in conjunction with forest biomass simulators to fully assess the dynamics of forest biomass.

Conclusions
The CWDDAT model used to assess CWD decomposition in forests was evaluated against wood mass loss observed from nine sites with different ecological conditions across the continental USA. Four quantitative model evaluation methods consistently demonstrated that the model is an effective tool for assessing CWD decomposition and the associated C fluxes in forests. Simulating CWD decomposition following Hurricane Hugo also demonstrated that this model is applicable for simulating responses at the watershed-scale. CWDDAT was built to reflect biogeochemical processes mediating wood decomposition; the calibration, validation and application simulation results demonstrate that those processes are sensitive to the primary factors known to regulate wood decomposition. The CWDDAT model incorporates both micro-organisms and arthropods as mediators of CWD decomposition, providing a basis for assessing the contributions of decomposer communities.
The CWDDAT model provides the basis for simulating the C cycle associated with forest CWD. The tool has considerable utility for assessing functional linkages between CWD decomposition and soil C, thereby providing a basis to the interactions of forest management, disturbance regimes (e.g., tropical storms), or climate change on this important component of the forest C cycle. Further evaluations of the model are warranted, and we hope that the findings reported here stimulate further studies to further strengthen the basis for model development and applications.