Examining Naturogenic Processes and Anthropogenic Influences on Tree Growth and Development via Stem Analysis: Data Processing and Computational Analytics

: The objective of this study was to develop a stem analysis data processing and computational algorithm and associated software suite that was (1) applicable to temperate and boreal forest tree species, (2) mathematically consistent with excurrent tree stem geometric and allometric principles, (3) compatible with data structures obtained using proprietary and non-proprietary imaging systems, and (4) executable on Windows ® -based operating systems. Computationally, the suite denoted SAP (Stem Analysis Program), deployed sectional-specific formulae that were in accord with the following geometric assumptions: (1) stump section was treated as a solid of revolution resembling a cylinder; (2) sections between the stump and the tip were treated as a solid of revolution resembling a frustum of a cone for sections with continuous annual increments, otherwise treated as a cone; and (3) tip section was treated as a solid of revolution resembling a cone. The algorithm also corrected for the slant-based sectional length measurements using Pythagorean Theorem and eliminated the requirement to predict age-specific apex height development through the use of a linear interpolation procedure. Based on input data structures consisting of annual ring-width xylem sequences measured from cross-sectional disk samples acquired at multiple positions along the tree’s main stem, the suite produces a broad array of output, inclusive of radial and longitudinal ring-width sequences, apical growth increments, annual and cumulative sectional and cumulative volume production patterns, and historically reconstructed stem taper profiles. In total, the SAP creates six output data files for each tree analyzed: (1) input data reference summary (e.g., geometric mean ring-widths and resultant radii for each cross-section);


Introduction
Stem analysis, the demarcation and subsequent analysis of annual tree ring sequences that are encased within the main stem of temperate and boreal tree species, provides an analytical framework for examining naturogenic processes and anthropogenic effects on tree growth and development. As exemplified by numerous examples, stem analysis has been one of the principal approaches used in detecting and assessing historical growth and development patterns of excurrent tree species. Those examining naturogenic processes would include studies quantifying the effects of edaphic variation [1], resource competition [2], abiotic wind stress [3], insect defoliation [4], drought [5], and pathogens [6,7] on tree growth, whereas anthropogenic based investigations would include studies assessing tree growth responses to intensive forest management treatments, such as pre-commercial thinning [8] and commercial thinning [9], atmospheric pollution [10], and global warming [11]. Stem analysis has also been an essential tool for reconstructing ecological disturbance histories (e.g., [12,13,14]) and past climate variability within forest ecosystems (e.g., [15,16]). Furthermore, tree ring sequences have been used as terrestrial-based proxies in reconstructing past continental-scale temperatures and hydroclimatic changes in order to detect anthropogenic-generated climate change and quantifying the associated effects (e.g., [17,18]).
Contextually, stem analysis approaches could be categorized into one of three classes, based on the level of complexity deployed. The first level includes non-destructive approaches, where analyses are focused on the examination of radial annual-ring increment sequences extracted at a single stem position, e.g., at breast height (1.3 m), obtained from single or multiple sample tree(s), in order to draw inferences regarding temporal patterns of cross-sectional growth (e.g., diameter and basal area increments). The majority of research investigations within the fields of dendrochronology, dendroecology, dendroclimatology and dendrohydrology, which involve temperate and boreal forest tree species, could be classified as such (e.g., [19,20]). The standardization of data management protocols and methodologies in terms of cross-dating procedures and statistical approaches, inclusive of the provision of the prerequisite software analogues, are well developed for this analytical class (sensu Laboratory of Tree-Ring Research, University of Arizona, Tucson , USA). Furthermore, as a consequence, these approaches have received considerable acceptance and deployment, leading to the creation of well-populated open-access global repositories of standardized tree ring chronologies with pertinent metadata attached (e.g., International Tree-ring Data Bank [21]). Analysis of these data sets within the context of climate change have significantly advanced the understanding of the key drivers of tree growth, leading to more informed predictions of the consequences of anthropogenic-induced climate change across the globe (e.g., [22]). The second analytical class would include those approaches involving destructive stem analysis, where sample trees are felled and sectioned in order to obtain (1) radial annual-ring counts (cambial age) at multiple stem heights (e.g., branch nodes), so that that temporal height growth patterns of codominant and dominant trees can be reconstructed and subsequently quantified, yielding site index functions, or (2) outside cross-sectional diameters at multiple stem heights for examining and modeling stem taper and morphological-based allometric relationships. Numerous methodologies have also been established for this analytical class in terms of sampling designs, computational frameworks and biometrical modelling (e.g., [23,24,25,26]).
The third analytical class, which is the most informative and generally represents the "gold standard" in terms of analytical scope and inference potential, includes those approaches which attempt to reconstruct the entire growth and developmental tree-level histories through the acquisition of radial annual-ring increment sequences measured from multiple cross-sectional disks positioned at fixed or variable stem positions over the entire length of the main stem. Although this complete (whole tree) stem analysis approach is the most logistically and analytically challenging to apply, it does enable the 3-dimensional visualization of annual radial and vertical growth layer profiles, and consequently provides for an expanded inference space in terms of analyzing and interpreting the effects of intrinsic and extrinsic influences on temporal growth patterns, allometry and stem morphological development. The conceptual, analytical and inferential frameworks for the complete stem analysis approach were largely established through the efforts of a series of pioneering growth and morphological investigations of red pine (Pinus resinosa Ait.) in central Canada during the early 1950s to mid-1960s, by Duff and Nolan [27][28][29] and Forward and Nolan [30][31][32][33]. However, the computational burden of reconstructing temporal, radial and vertical patterns of tree growth and development from radial annual-ring increment sequences measured on multiple disks, can present a hindrance to its application in forest research investigations and operational forest management. Consequently, in order to potentially overcome this deterrent, the focus of this study was to present the prerequisite computational analytics and an associated software suite that was mathematically consistent with excurrent tree geometry, compatible with current software and hardware environments, and generated the response metrics that could facilitate comprehensive analyses of the effects of naturogenic processes and anthropogenic influences on tree growth and development.

Stem Analysis Computational Algorithm: Structure and Development
The VisualBasic.NET software suite, denoted SAP (Stem Analysis Program) herewith, was developed within the Visual Studio ® development environment (VS.NET (Version 2012), Microsoft Inc. (MS), Redmond, WA, USA) and utilized a number of third party graphical and tabular reporting tools (i.e., ProEssentials (Version 6), GigaSoft Inc., Keller, TX, USA, and Input Pro for Windows Forms (Version 2.5), Charlotte, NC, USA, respectively). The main computational engine of the suite is the Fortran-based executable program, denoted SAPFOR.EXE (Stem Analysis Program -FORtran). SAPFOR is an improved and expanded version of the stem analysis algorithm previously presented by Newton [34,35]. Briefly, the SAPFOR program was written in Fortran using the Lahey/Fujitsu Fortran 95 compiler (www.lahey.com) and consisted of a main program and two subroutines. The primary function of the main program is to compute the geometric mean ring-widths and resultant radii, calendar date, and age from the raw data file, which consists of the annual ring width sequences from the multiple cross-sectional disk samples for a given sample tree. The first subroutine, denoted GDB for Growth Diameter and Basal area analysis, then computes radial growth estimates for the cross-section at breast-height, specifically. Subsequently, the second subroutine, denoted GVA for Growth Volume and Allometric analysis, computes the sectional and cumulative growth estimates, and stem developmental profiles. The underlying structure, computational assumptions, and mathematical formulae are detailed explicitly in Appendix A.
In terms of data processing, the suite also converts WinDendro TM -based (www.regentinstruments.com) Excel ® input files (W3*.xlsx (Version 3 format) or W4*.xlsx (Version 4 format) into the required Fortran-based configuration for use within the SAPFOR. Specifically, the program processes the selected WinDendro TM -based Excel ® input file by transcribing it into a 2dimensional rectangle array format after removing (parsing) all extraneous data columns and infilling blank cells with null zero values, resulting in a correctly structured input file. Generically, the structure of the input file (denoted SWijklmn for Stem Width increments, where ijklmn is a numeric code for identifying the sample tree, and the plot and (or) stand it was sampled from) consists of the each disk's height position (m) and its total age (years), followed by the cumulative length of the each radial sequence (mm or cm) and each annual ring-width increment (mm or cm) along each radii, starting from the pith and progressing outward to the outmost annual ring. Note that, as long as the increment data is presented in this structure, non-Windendro TM generated input files may also be processed by the program. Refer to input file SW11002.xlsx in Supplementary: Material S1 for an applicable example.
More specifically, employing radial xylem sequences derived from cross-sectional disks located at multiple sample points along the stem of excurrent (single-stem) forest tree species, the program calculates: (1) annual increments and one periodic increment for diameter and basal area growth at breast-height (1.3 m) in both absolute and relative terms; (2) annual volume increments and one periodic volume increment for each section in both absolute and relative terms; (3) annual increments and one periodic increment for volume growth in both absolute and relative terms for all sections combined; and (4) annual diameters and corresponding heights for each year. Note that the time interval pertaining to the periodic measurement is defaulted to unity, which in essence treats the periodic computation as an annual one. Currently, this option is included to enable further program expansion for addressing situations where (1) biotic or abiotic damage has resulted in partial disk samples, where all annual increment measurements along the radial sequence are not available (e.g., pathogenic incipient decay within the inner core) or (2) logistical limitations or study objectives dictate the annual measurement of only a fixed outer portion of the radial sequence (e.g., analyzing only post-treatment annual growth response patterns arising after the implementation of a midrotation silvicultural treatment).
Relative measures of growth included relative growth rate [36], growth percentage [37], specific volume increment [28], and relative production rate [38]. Relative growth rate and growth percentage express absolute growth as a proportion of initial size. Specific volume increment is defined as the rate of stem volume production per unit of stem surface area and, hence, indirectly reflects the ability of a tree to supply photosynthate (i.e., assimilates produced through photosynthesis) to the vascular cambium for xylem production. Relative production rate reflects the rate of change in stem volume growth between successive growth years. Conceptually, relative growth rate, growth percentage, specific volume increment and relative production rate, attempt to remove the effect of cumulative tree size, thus providing an unbiased method to assess growth performance among sample trees, irrespective of their size differences.

Exemplification of the Utility of the Resultant SAP Suite
In order to demonstrate the utility of the SAP suite, the processing of the annual ring-width increment sequences obtained from cross-sectional disks obtained from a jack pine (Pinus banksiana Lamb.) sample tree, and the associated analysis of the derived growth layer and specific volume increment profiles, are examined within the context of (1) identifying intrinsic developmental patterns generic to excurrent tree species, and (2) examining longitudinal growth responses to thinning treatments. Silviculturally, the selected sample tree had grown in a natural-origin, evenaged, monospecific (> 90% jack pine by basal area) stand that had been established via natural seeding in the year following a wildfire in 1945. The stand was located in the Tyrol Research Forest, which falls within the Nipigon Forest Section (B.10) of the Canadian Boreal Forest [39]. The stand was situated on a dry-to-moderately fresh sandy-textured conifer site type (ES 13 [40]) and was aboveaverage in terms of its productivity (site index = 19 m [41]). In 1962, the stand was subjected to an operational precommercial thinning treatment when the sample tree was 17 years of age. This treatment reduced the stand density from approximately 8500 stems/ha to 2419 stems/ha (nominal 2.0 × 2.0 m spacing). Thirty-five years later, the stand was subjected to an experimental thinning treatment in 1997, when the sample tree was 52 years of age. This treatment reduced the stand density from approximately 1669 to 748 stems/ha (3.7 × 3.7 m nominal spacing). At the individual tree level, this treatment reduced the number of immediate competitors (i.e., trees within a 3 m radius of the sample tree) from five relatively smaller and larger sized (height-based) competitors down to two smaller sized competitors. Stand-level mensurational parameters at the time of the stem analysis sampling (autumn 2015) for stand density, quadratic mean diameter, mean dominant height, basal area, total volume, merchantable volume and relative density were, respectively: 622 stems/ha In-forest stem analysis sampling procedures consisted of felling, delimbing, measuring for crown length and total stem length, and subsequently sectioning the main stem at the following height positions: stump (0.3 m), breast-height, and at 10% height intervals (10, 20,…,90). A crosssectional disk was then removed from each of these positions, yielding a total of 11 disks. The section lengths were measured and the disks were then returned to cold storage until laboratory processing. In-laboratory procedures for each disk consisted of (1) sanding its cross-sectional surface, and (2) determining and marking its geometric mean diameter from which two annual ring width radial sequences were measured via the Windendro TM imaging system (Version 2012b; www.regentinstruments.com). The raw data files were saved in the Version 4 format (denoted W411002.xlsx; see Supplementary Material S1). In order to demonstrate the SAP suite employing all three possible input data file configurations, a Windendro TM Version 3 file and a SW data file, were also created (denoted W311002.xlsx and SW11002.xlsx, respectively; see Supplementary Material S1). Note, the file name (identifier) used for this example was 110022 (ijklmn) which indicated that the sample tree number was 022, sampled from plot number 110.

Overview of the Resultant SAP Suite
Contextually, the SAP suite is a data-base driven MS Access application, consisting of a VBbased program in which input data files are assessed for their size parameters and subsequently manipulated in order to comply with the input requirements of the Fortran-based computational algorithm (SAPFOR). Consistent with widely accepted bole shape principles for excurrent tree species (sensu [42]), the SAPFOR-based computational formulae were based on the following geometric identities: (1) stump section was treated as a solid of revolution resembling a cylinder; (2) sections between the stump and the tip were treated as a solid of revolution resembling a frustum of a cone in cases where the annual increment was continuous throughout the section, or a cone in cases where the annual increment was not continuous throughout the section; and (3) tip section was treated as a solid of revolution resembling a cone. The algorithm also corrected for the slant-based sectional length measurements using Pythagorean Theorem and eliminated the need to predict height at a given age through the employment of a linear interpolation procedure. Based on input data structures consisting of annual ring-width xylem sequences measured from cross-sectional samples acquired at multiple positions along the main stem, the resultant suite produced a broad array of output, including radial and longitudinal ring-width sequences, apical growth increments, annual and cumulative sectional and cumulative volume production patterns, and historical stem developmental profiles. Graphics were also generated for data verification (e.g., evaluating the measured annual ring-width data by radii and disk) and interpreting temporal growth and developmental patterns (e.g., analyzing vertical growth layers and specific volume increment profiles by age or year).

Input File Format and Structure, Program Execution, and Generated Output
The input file containing the ring-width sequences must be formatted as a complete rectangular matrix with all cells populated if the input file is manually created from an imaging system other than Windendro TM . The required structure for manually inputted annual-ring sequences for each cross-sectional disk is the following: age of the kth cross-section (year); height of the kth cross-section (m); radial-specific ring-width sequences, each consisting of the total radial length (cm), and annual ring widths (mm) arranged in a pith-to-bark order. The order of the disks must also follow a specified pattern: sequences obtained from the stump cross-section, then the breast-height cross-section, then the sections between the breast-height and the tip ordered by increasing height, terminating with the tip section (refer to file SW11002.xlsx in the Supplementary Material S1 for an example of the correct structure). Otherwise, in the case of Windendro TM generated Excel ® files, either in Version 3 or 4 format, the suite automatically completes this data management step. Specifically, it parses and rearranges the raw Windendro TM data files and then completes the necessary preliminary calculations resulting in the correctly structured SGijklmn.DAT (Stem Growth increment) input file. In terms of analytical scope, one to eight radial sequences per disk, and up to 16 cross-sectional disks per tree, can be accommodated.
Upon execution, the SAP suite creates six output data files which can be readily exported for further analysis: (1) GFijklmn.DAT (Geometric Formatted file), consisting of the geometric mean ringwidths and resultant radii for each cross-section; (2) BIijklmn.xlsx (Breast-height Increment data file), consisting of the diameter and basal area growth metrics for the breast-height cross-section only; (3) SSijklmn.xlsx (Stem Sectional data file), consisting of volume growth metrics for each section; (4) STijklmn.xlsx (Stem Total data file), consisting of the cumulative volume growth metrics based on all sections combined; (5) SFijklmn.xlsx (Stem Form data file), consisting of the historical reconstructed stem taper profiles by year; and (6) REijklmn.xlsx (REference data file), consisting of a labeled summary of all five files collated together. Graphics are also produced for data verification and preliminary growth analysis. Specifically, in terms of data verification, the actual increment widths by year for each radii on each disk are presented to assist the analyst in verifying the input data derived from the imaging system (e.g., graphical detection of outliers or incorrect measurements). An average trend line is also produced for each disk, in order to visually evaluate the general growth trends within and among disks. In terms of preliminary analyses and to aid interpretation, the growth layer and the specific volume increment profiles for the entire tree are also graphically illustrated.
By way of example, in reference to the sample tree demonstration, the input files were placed in a separate computer directory and, upon execution of the SAP suite, loaded into the suite using the "Data Loader" tool bar ( Figure 1). Note, the end-user must specify the (1) sampling year, (2) number of years for the periodic increment (currently the value is set to a default value of unity), (3) annual increment measurement units (cm or mm), and (4) identify the source folder (C:\Users\pnewton\Documents\SAP-BETA-EXEC\TESTTREE) and associated increment file (e.g., W4110022.xlsx). Following the selection of the input data file and setting the above parameters, the SAP suite provides the original source input file from the imaging system (e.g., W411002.xlsx) and its associated SAP-compatible reformatted input file, followed by the generated six output data files. Note, the RE110022.xlsx file which includes all 5 of the individual labelled output files is presented in the Supplementary Material S2. Within the SAP graphical user interface (GUI), all of these files are presented as a set of tables under the "Data" tool bar ( Figure 2). The data verification, growth layer and specific volume increment profiles graphics are presented as a set of figures under the "Graph" tool bar ( Figure 3). Additionally, for the profile graphs, end-users have the option of selecting increment profile illustrations for one or more of the years analyzed.    Figure 4 illustrates the actual increment widths by year for each of the two radii on each of the 11 disks obtained from the jack pine sample tree. Visual examination of the radius-specific increment values, along with the mean increment trend line by disk, revealed no discernable outliers in the raw increment data. The resultant growth layer profile for the tree is presented in Figure 5. The growth layer profiles were consistent with expectations: (1) for any given growth layer, the radial increment exhibited a general increase from the stump upwards and then declined rapidly as the apex was approached; empirically the maximum increment occurs at the 80% height percentile position, which corresponded to a point approximately 1/3 into the basal portion of the live crown, (2) as the tree aged, the point at which the radial increment was minimal occurred at a height of approximately 20%; and (3) although this increment profile pattern became more pronounced with age, presumably due to increasing intraspecific competition, the width of the radial increments temporally increased following the thinning treatments, particularly after the second thinning ( Figure 6). These longitudinal increment patterns follow those espoused by Farrar [43]: i.e., (1) the point of the maximum increment occurring at the stem position where the crown density was the highest (assumed to occur at an 80% height percentile for the sample tree analyzed in this study), (2) the increment increasing from the apex to the mid-crown position, but thereafter declining with decreasing height until a minimum is approached (assumed to occur at a 20% height percentile for the sample tree analyzed in this study), and (3) increments then slightly increasing from that point until the base of the tree is reached. This increase in the lower bole (stump region) has been associated with mechanical stability requirements of excurrent tree species (i.e., butt swell; sensu [43]).More specifically, the growth layer profile indicated a period of rapid early growth, reflecting generally, the open-growing conditions lasting approximately 12 years in duration (1946)(1947)(1948)(1949)(1950)(1951)(1952)(1953)(1954)(1955)(1956)(1957)(1958). Examination of the within-tree specific volume increment profiles (Figures 7-9) indicates an increasing apex-to-base growth compression during the years preceding the precommercial thinning treatment (1959)(1960)(1961), immediately followed by a pattern consistent with thinning release effects. Specifically, the thinning treatment elicited an increase in the increment widths towards the upper part of the stem (upper two thirds) which was largely maintained until the effects of intraspecific competition took hold at approximately 15 years after the PCT treatment (1978). A period of continuous growth reduction throughout the entire stem continued until the time of the second thinning (1997).

Interpretation Guidance and Related Inferences Derived
Year Figure 4. Data verification for the 70 year-old jack pine example tree: illustration of the annual increment widths (mm) by year, radii (denoted R1 and R2 in legend) and cross-section. Note, the (1) section denotations (S: 1-11; y-axis) refer to the sections sampled at stump height, breast-height, and 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80% and 90% relative height positions, respectively, (2) mean trend line for each cross-section is superimposed, and (3)   . Selected growth layer profiles for the example tree illustrating the longitudinal distribution of radial annual increment widths with increasing stem height, for the three years prior to the 1st and 2nd thinnings (1959-1961 and 1994-1996, respectively), the year of thinning treatments (1962 and 1997), the following three years after the 1st and 2nd thinnings (1963-1965 and 1998-2000, respectively).  Similar to the first treatment, this second treatment elicited an increase in the increment widths towards the upper part of the stem (upper two thirds) until the effects of intraspecific competition took hold at approximately 7 years post-treatment (2004) and continued for the following 11 years (2005-2015). Note, these patterns before and after the second thinning are also evident in the treelevel tabular responses, as shown in Figure 2, in terms of relative growth rate and specific volume increment. The specific volume increment profiles for each 3 year period before and after the two thinning treatments, in addition to the last 3 years of growth, as shown in Figure 7, exhibited a consistent apex-to-base C-shaped pattern. More precisely, based on the vertical distribution of growth efficiency, maximum values occurred within the uppermost portion of the stem but then rapidly declined with decreasing stem height until reaching a minimum value near the mid-lower portion of the stem, and then slightly increasing as the stem base was approached. Although lacking the distinctive reverse S-shape reported for jack pine trees following thinning treatments by Goudiaby [9], the observed patterns are in general accord with the conceptual constructs advanced for conifers [43]: (1) specific volume increment maximizes near the point at which the concentration of foliage is the greatest (e.g., 33% live crown percentile (vertical position)) due to the proximity of the photosynthate production centre; (2) the distribution of photosyntate then declines with decreasing height, resulting in declines in specific volume increments; and (3) photosynthate is allocated in a manner which maximizes survival probabilities. Examples with respect to the last construct could include (1) increased allocation to the mid-portion of the stem in order to increase sapwood area for servicing the higher hydrological requirements arising from increased foliage production following thinning, and (2) increased allocation to the stump region in order for the tree to adjust its taper for addressing the additional stability requirements arising from the increased crown mass (as plausibly evident in this study from the slight increase in specific volume increment in the lowest portion of the stem). This latter effect results in a decreased height/diameter ratio [44] thus yielding increased stability, and is similar in theory, to that proposed for the development of butt swell (sensu [43]).

The SAP Suite
Computationally, the accuracy of the growth estimates and developmental profiles obtained through stem analysis are partially dependent not only on the sampling strategy, in terms of the number of cross-sections sampled and the number and type of radial sequences measured (e.g., [4,25,45]), but also on the applicability of the underlying geometric assumptions embedded in the computational algorithms used. For the SAP suite, the geometric assumptions deployed within the core computational algorithm (SAPFOR) were consistent with those commonly accepted for excurrent tree forms [42]: geometric solids of revolution resembling a (1) cylinder for the stump section, (2) frustum of cone for central sections, where annual radial increments were continuous throughout, or, alternatively, a cone, when annual radial increments were not continuous, and (3) cone for the apex section (tip). The algorithm employed a tangential height growth projection method to estimate the location of embedded terminal height growth positions. Specifically, deploying Carmean's linear interpolation procedure for all non-terminal sections (Equation (A8); [46]) and Newberry's [47] modification of Carmean's method for the apex section (Equation (A9)). This approach has been shown to be acceptable in terms of the tenability of the underlying assumption (tangential linear growth) and predictive precision when applied to excurrent forest tree species (e.g., [26,48,49,50]). Additionally, in the field, section lengths are measured along the outer periderm surface, which geometrically represents the slant length. However, this is not equivalent to the vertical length measurement required for the conic computations. The SAPFOR utilizes Pythagorean Theorem to estimate vertical lengths from the measured slant length and radial measurements before commencing with the volumetric calculations (e.g., Equation (A7)). Likewise, based on this theorem, the vertical length estimates along with radial measurements are used to calculate the unobserved slant lengths of embedded growth years, which is a prerequisite for calculating internal lateral surface areas (e.g., Equation (A13)). Collectively, these advancements provide a mathematically sound foundation for deriving growth estimates and reconstructing historical developmental profiles when deploying the stem analysis method.
Furthermore, the inclusion of innovative relative growth output measures, such as the specific volume increment, relative growth rate, and relative production rate, increases the comprehensiveness of the response metrics available to the analyst when evaluating anthropogenic influences and (or) examining naturogenic processes. Notably, however, the growth percentage and relative growth rate express growth as a proportion of initial absolute size, and hence caution should be exercised in their interpretation, given that they can intrinsically decline with increasing size irrespective of external causes (sensu [51]). Conversely, relative production rate reflects the rate of change in absolute growth between successive years or periods, and thereby largely eliminates the potentially confounding effect of initial size. In terms of year-to-year temporal variability, specific volume increment and relative growth rate inherently exhibit more stability than relative production rate when analyzing patterns over time, given that the cumulative memory of size is carried forward within these indices. Consequently, at the analytical stage, consideration of the deploying functional growth analysis methods may be a feasible option since this approach would largely mitigate the year-to-year temporal variation commonly observed when using relative production rate (sensu [38]).

Analytics of Whole Tree Stem Analysis
Logistically, three analytical steps are required to implement whole tree stem analysis once the sample trees are selected: (1) extracting cross-sectional disks at fixed or variable height positions along the main stem, (2) given (1), measuring their radial annual-ring width sequences along selected radii, and (3) given (2), processing the resultant data streams and calculating the growth and developmental metrics of interest. Specifically, with respect to the first step, the actual sampling design deployed is largely driven by the underlying objective of the analysis. For example, studies examining volumetric growth patterns have found that a variable-based, percent-height sampling approach to be optimal in terms of minimizing estimation error (e.g., sampling disks at approximately 10% height intervals [4,35]). Apart from the specifics of the overall sampling design employed, most stem analysis investigations require the measurement of one or more radial annual ring-width sequences from each of the cross-sectional disks acquired. A common practice is that, when disks arrive in the laboratory from the field, they are placed in refrigerated storage in order to avoid pathogenic-caused decomposition or excess drying and associated wood checking. The first processing step commonly involves power sanding or planning the cross-sectional surfaces in order to enable the demarcation of the annual ring-width increments along the selected radii. Various elliptical-based radial selection strategies have been devised, given that cross-sectional surfaces are rarely perfect circles. These commonly include randomly selecting two semiaxes along the axis corresponding to the geometric mean diameter, two semiaxes in which the minimum and maximum semiaxes are measured, two semiaxes along the major axis, or four semiaxes along both the minor and major axes (sensu [35]).
The second step consists of measuring the selected radial axes, employing one of the various semi-automated systems, and obtaining the annual-ring width sequences. Historically, advances in instrumentation and computer science over the last four decades have led to a continuous pattern of innovation in the development of automated systems for tree ring measurement and analysis. Representative systems along this innovation pathway include the Digital Positometer [52], ADDO-X System [53], Holman-Digimicrometer [54,55], Tree Ring Scanner [56], Tree Ring Measurement System [57], TRIM [58,59], DendroScan [60], MacDendro TM and WinDendro TM Image Analysis Systems ( [61] and [62], respectively), and Advanced Tree-Ring Image Capturing System (ATRICS [63]). Although this fragmented and discrete pattern of hardware development has hindered the establishment of common standards across systems in terms of input sample requirements, measurement protocols, data structures, computational analytics and software-based algorithms.
The third step involves developing and (or) deploying applicable data processing and computational algorithms which yield the required indices for drawing inferences regarding the effect of naturogenic processes and anthropogenic impacts on tree growth and development. Although a number of very useful software computational algorithms have been developed in the past (e.g., [34,35,58,61,[63][64][65][66][67][68]), they have been generally non-generic in their applicability and are specific in terms of their input data structure requirements, underlying computations and associated assumptions employed, growth and morphological metrics produced, and tabular and graphical outputs generated. Furthermore, even though only a few system-specific processing and analysis software programs are still currently available and actively supported (e.g., XLSTEM part of the WinDendro TM Image Analysis Systems (Regent Instruments Inc, Quebec City, Quebec, Canada) and TSAP-Win TM within Digital Positiometer System (Rinntech ® , Heidelberg, Germany)), the accessibility of comprehensive non-proprietary computational suites, that are able to readily produce output metrics that are consistent with the conceptual foundation and analytical stem analysis framework advanced by Duff, Nolan and Forward [27][28][29][30][31][32][33], and compatible with current computer hardware and software operating systems, are generally lacking. Hence the SAP software suite, which is (1) applicable to a wide range of temperate and boreal forest tree species, (2) mathematically consistent with excurrent tree geometric and allometric principles, (3) compatible with data structures obtained using both proprietary and non-proprietary imaging systems, and (4) executable on the dominant microcomputer operating system, represents a consequential contribution in fulling this knowledge gap. Specifically, the development of a software suite for processing and analyzing stem analysis data that (1) in addition to the commonly computed growth and developmental metrics, provides sectional and cumulative estimates of cambial stem surface area, specific volume increment [28], relative growth [36] and relative production rate [38], which are increasingly used to analyze naturogenic processes (e.g., [2]) and anthropogenic influences (e.g., [9]), (2) includes the computational innovations introduced for excurrent forest tree species (e.g., correcting for slantbased sectional length measurements and deploying an internal linear interpolation algorithm for height estimation), and (3) is compatible with current and foreseeable computer hardware and software environments, represents an incremental step in advancing stem analysis methodology in forest research and management investigations.
Moreover, the SAP suite was developed by (1) generalizing the input data requirement so that the core computational algorithm could accept a diverse range of input data structures, (2) expanding the computational structure in order to generate enhanced tabular and graphic outputs, and (3), given (1) and (2), integrating the resultant algorithm into a hybrid VB-Fortran programming platform. Furthermore, the employment of the revised Fortran-based program (i.e., SAPFOR) eliminated the requirement to transcribe and recode the core computational engine into the VB.NET programming language. Consequently, the analytical innovations originally introduced in stem analysis algorithm [34,35] were preserved, including (1) correction for slant-based sectional length measurements using the Pythagorean theorem, which was commonly overlooked in previous stem analysis programs, and (2) eliminating the requirement to externally predict terminal heights using an embedded linear interpolation procedure. In terms of deployment, the scope of the SAP suite's applicability currently includes temperate and boreal excurrent forest tree species. The program is able to accept input data structures derived from a variety of imaging systems, provided that they are formatted correctly. The executable variant of the SAP suite requires no additional software programs other than Windows ® 7 SP1 or newer and Microsoft's NET Framework 4.6 to run. All the other executable software required for the GUI and graphical and tabular reporting are embedded. Further suite developments will focus on enabling the optional processing of periodic-based increment sequences and establishing a webbased knowledge transfer site for facilitating the distribution of the executable variant of the suite to non-commercial end-users.

Conclusions
The conceptual, computational and analytical foundation of the stem analysis approach is well documented within the scientific literature. The approach has been widely used as the primary investigative tool for quantifying the growth, development and allometric responses of excurrent trees to a broad array of anthropogenic influences and naturogenic processes. However, the evolution of instrumentation, along with the processing and analytical algorithms required for implementing the stem analysis approach has been largely fragmented over the last 4 decades, as evident from the infrequent short-term periods of rapid activity interspersed by long periods of inactivity. As a result, the computational side of the stem analysis method has generally not kept pace with technological advancements in imaging systems, computer hardware or software design, nor in terms of incorporating new growth and developmental response metrics useful in forest research and management. Consequently, the SAP suite presented in this study which (1) overcomes the significant computational challenges involved in processing radial-based sequences extracted from multiple cross-sectional disks obtained throughout the entire tree stem, (2) analytically addresses the complex dimensionality of excurrent tree species via the deployment of a mathematical approach that is consistent with underlying geometric and allometric principles, (3) is compatible with data structures obtained using both proprietary and non-proprietary imaging systems, and (4) is executable on the globally dominant microcomputer operating system, represents an incremental contribution in terms of advancing the utility of the stem analysis approach in forest research and management investigations.

Conflicts of Interest:
The author declares no conflict of interest. The funder had no role in the design of the study, in the collection, analyses, or interpretation of data, in the writing of the manuscript, or in the decision to publish the results.

Appendix A: Analytical Structure, Computation Formulae and Underlying Geometric Assumptions Deployed in the SAPFOR Algorithm
The SAPFOR algorithm computes: (1) annual increments and one periodic increment for diameter and basal area growth at breast-height (1.30 m) in absolute and relative terms; (2) annual volume increments and one periodic volume increment for each section in absolute and relative terms; (3) annual increments and one periodic increment for volume growth in absolute and relative terms for all sections combined; and (4) annual diameters and corresponding heights for each section. This program was written in Fortran using the Lahey/Fujitsu Fortran 95 compiler and consisted of a main program and two subroutines, denoted the Growth Diameter and Basal Area Analysis Subroutine (GDB) and the Growth Volume and Allometric Analysis Subroutine (VGA). The primary function of the main program is to compute the geometric mean ring-widths and resultant radii, calendar date and age from the input data file, and store the results in a transitional input data file (denoted GFijklmn; Growth Formatted data file). The GDB subroutine computes radial growth estimates for the cross-sectional disk located at breast-height (1.3 m), utilizing the growth data file (GFijklmn) and deposits the resultant diameter and basal area growth estimates into the Breast-height Increment data output file (denoted BIijklmn). The VGA subroutine computes the sectional and cumulative growth estimates, and stem taper developmental profiles, similarly using the growth data file values (GFijklmn). The resultant output is deposited into one of three output files: (1) SSijklmn (Stem Sectional data file) consisting of volume growth analysis for each section; (2) STijklmn (Stem Total data file) consisting of the cumulative volume growth analysis based on all sections combined; and (3) SFijklmn (Stem Form data file) consisting of the reconstructed annual-based taper profiles.

Growth Diameter and Basal Area Analysis Subroutine (GDB)
The GDB subroutine computes radial growth estimates for the cross-section at breast-height for sample trees having non-zero cross-sectional data. These estimates include diameter increment (cm/(year or period), basal area increment (m 2 /(year or period)), relative diameter increment (loge(cm)/(year or period) [36]), relative basal area increment (loge(m 2 )/(year or period) [36]), diameter and basal area growth percentage (%·10 −2 /(year or period) [37]), and relative diameter and basal area production rate [38]. More specifically, absolute and relative diameter growth estimates were computed via Equations (A1-A3): is the relative diameter increment (loge(cm)/(year or period)) at time t [36], and ( ) GP t DI is the diameter increment expressed as a growth percentage (%·10 −2 /(year or period)) at time t [37]. Similarly, absolute and relative basal area growth estimates were computed via Equations (A4-A6): BAI is the basal area increment (m 2 /(year or period)) at time t, is the relative basal area increment (loge(m 2 )/(year or period)) at time t [36], and is the basal area increment expressed as a growth percentage (%·10 −2 /(year or period)) at time t [37].

Growth Volume and Allometric Analysis Subroutine (GVA)
The GVA subroutine computes the sectional and cumulative growth estimates, and stem developmental metrics. The computations are based on the following assumptions: (1) the stump section was considered a geometric solid of revolution resembling a cylinder; (2) sections between the stump and the tip were considered geometric solids of revolution resembling a (i) frustum of a cone when annual increments were continuous throughout the section, or (ii) cone, in the case where annual increments were not continuous throughout a section; (3) the tip section was considered a geometric solid of revolution resembling a cone; and (4) vertical section lengths were based on the Pythagorean theorem (Equation (A7); [34,35]): where V L is the vertical sectional length (m), S L is the slant (field measured) sectional length (m), and B R and T R are the radii (m) at the bottom and top of the section, respectively. Carmean's [46] linear interpolation procedure was used to estimate height at a given age for all non-terminal sections (Equation (A8); [34,35]), whereas Newberry's [47] modification was used for the terminal (top) section (Equation (A9) [34,35]): where Hij is the estimated total height (m) for the jth growth ring based on the sectioning point i, hi is the measured height (m) at the ith sectioning point, hi+1 is the measured height (m) at the ith+1 sectioning point, ri is the number of growth rings at the ith sectioning point, ri+1 is the number of growth rings at the ith+1 sectioning point, and j is the growth ring number based on the pith being the starting point (j=1,…,ri).

Sectional Computations (1) Stump Section
Volume and lateral surface area computations for the stump section were based on a geometric solid of revolution resembling a cylinder: Equation (A10) ( [34,35]) and Equation (A11), respectively.
is the volume (m 3 ) of the stump section at time t, values were then used to calculate absolute and relative growth rate estimates specific to the stump section (i.e., absolute volume increment and specific volume increment by year or period).
(2) Sections between the Stump and the Tip The geometric assumptions for the computations related to the sections between the stump and tip depended on the continuality of the increments throughout the section. Specifically, for sections in which increments are continuous throughout, volume and lateral surface area computations were based on a geometric solid of revolution resembling a frustum of a cone, Equation (A12) and Equation (A13), respectively [34,35].
Conversely, for sections in which increments are not continuous throughout, volume and lateral surface area calculations were based on a geometric solid of revolution resembling a cone, Equation (A14) and Equation (A15), respectively [34,35]. where is the volume increment (m 3 /(yr or period)) of the lth section at time t, and is the specific volume increment ((m 3 /m 2 )/(year or period)) for the lth section at time t [1,28].
(3) Apex Section Volume and lateral surface area computations for the tip section were based on a geometric solid of revolution resembling a cone, Equation (A18) and Equation (A19), respectively [34,35]. if increment is not continuous throughout section where 100 Tip volume increment was calculated by Equation (A20): where is the volume increment (m 3 /(year or period)) of the tip section at time t. These volume increments, combined with the lateral surface area estimates, yielded the corresponding specific volume increment values unique to the tip section.