Long-Term Ecosystem Nutrient Pool Status for Aspen Forest Harvest Simulations on Glacial Till and Sandy Outwash Soils †

: Sandy outwash and glacial till soils compose large amounts of public forestland due to historically poor agricultural yields. The outwash soils have low fertility, poor nutrient retention and are restricted from whole-tree harvesting (WTH) in some states, whereas the glacial till has medium nutrient retention and fertility, and is unrestricted from WTH. To assess the long-term sustainability of harvesting, a nutrient budget was constructed from ﬁeld measurements, the National Cooperative Soil Survey (NCSS) database, and literature values for stem-only harvesting (SOH) and WTH at a 45-year rotation length and 11 rotations were simulated. The budgets showed that SOH and WTH recovery years, or the time necessary for the inputs to match outputs through leaching and one harvest, exceeded common rotation lengths for both soil types under all weathering scenarios, and the average WTH reduced the total available rotations by one harvest. The large variation in soil nutrient pools and harvest removals complicated the ability to identify the difference between SOH and WTH early in the model, but differences became apparent with sequential harvests. The recovery years were 2–20 times the 45-year rotation length under all weathering rates. Taken together, models in this study bridge the gap between short- and long-term studies and bring into question the sustainability of WTH and SOH practices on nutrient-poor soils.


Introduction
Nutrient-poor and medium fertility soils, such as glacial outwash (Typic Udipsamments, Entic Haplorthods) and tills (Alfic Haplorthods), are abundant throughout the Great Lakes region, USA. The soils compose large percentages of managed public forestlands after being abandoned as the "lands nobody wanted" [1] and are commonly managed in a regime of even-aged aspen clearcut harvesting. Bigtooth aspen (Populus grandidentata Michx.) and quaking aspen (Populus tremuloides Michx.) are the most widely distributed trees in North America [2] and contain large amounts of Ca and Mg in woody tissues [3]. However, there is still much uncertainty as to the frequency of harvests that can be conducted in lower to medium fertility soils and the effects of harvest on long-term available nutrient stocks [4][5][6][7][8][9][10][11].
The shift towards including woody biomaterials in the economy has increased the demand for woody biomass yield in forest harvesting practices [12]. Intensive forest harvesting is a concern for nutrient sustainability in low fertility or nutrient-poor soils and examine seasonal changes in the effects of WTH to determine if there are seasons of operability wherein nutrient removals can be minimized due to observed differences in residuals from breakage [19]. To meet these objectives, seasonal input-output budget models of forest tissue and soil elemental pools were constructed. The harvest models considered inputs and outputs on a 45-year rotation length and compared SOH to winter, spring, summer, and fall WTHs. Each harvest was repeated over time to determine the long-term effects of the harvest. Recovery years, defined as how long it would take for inputs to replenish values lost to one harvest and leaching, were calculated for comparison with typical rotation lengths.

Overview
This study examined nutrient inputs and outputs in aspen forests occurring on harvestrestricted soil complexes in Wisconsin, USA. We leveraged soil data collected from 14 new soil pedons (Figure 1), national soil databases, and literature values (see appendix for details). A nutrient budget was constructed for the soils considering inputs as weathering and deposition and outputs of soil leaching and harvest removals in comparison to available nutrient pools similar to previous research [5,6]. Soil weathering has previously been measured on the sandy outwash and till soils considered in this study, and the depletion method rates were used for the annual weathering inputs [24]. The total deposition was derived from deposition maps ( [25]; Figure A1). Soil and tree tissue information was collected across the northern portion of the state of Wisconsin (Figure 1). The nutrient budget considered SOH and WTH during the winter, spring, summer, or fall, under three regional deposition budgets and three weathering scenarios. seasonal changes in the effects of WTH to determine if there are seasons of operability wherein nutrient removals can be minimized due to observed differences in residuals from breakage [19]. To meet these objectives, seasonal input-output budget models of forest tissue and soil elemental pools were constructed. The harvest models considered inputs and outputs on a 45-year rotation length and compared SOH to winter, spring, summer, and fall WTHs. Each harvest was repeated over time to determine the long-term effects of the harvest. Recovery years, defined as how long it would take for inputs to replenish values lost to one harvest and leaching, were calculated for comparison with typical rotation lengths.

Overview
This study examined nutrient inputs and outputs in aspen forests occurring on harvest-restricted soil complexes in Wisconsin, USA. We leveraged soil data collected from 14 new soil pedons (Figure 1), national soil databases, and literature values (see appendix for details). A nutrient budget was constructed for the soils considering inputs as weathering and deposition and outputs of soil leaching and harvest removals in comparison to available nutrient pools similar to previous research [5,6]. Soil weathering has previously been measured on the sandy outwash and till soils considered in this study, and the depletion method rates were used for the annual weathering inputs [24]. The total deposition was derived from deposition maps ( [25]; Figure A1). Soil and tree tissue information was collected across the northern portion of the state of Wisconsin ( Figure 1). The nutrient budget considered SOH and WTH during the winter, spring, summer, or fall, under three regional deposition budgets and three weathering scenarios. The major nutrient-poor outwash (Typic Udipsamments-Entic Haplorthods), glacial till (Alfic Haplorthods), and low silt + clay soils from the gridded Soil Survey Geographic (gSSURGO) database. The sample pedons were collected with NRCS soil scientists to increase available pedons in outwash soils, and NCSS pedons were retrieved from the database for the outwash and glacial till soils. The aspen species within this study have large ranges across North America and are shown within the United States (derived from [26]). The major nutrient-poor outwash (Typic Udipsamments-Entic Haplorthods), glacial till (Alfic Haplorthods), and low silt + clay soils from the gridded Soil Survey Geographic (gSSURGO) database. The sample pedons were collected with NRCS soil scientists to increase available pedons in outwash soils, and NCSS pedons were retrieved from the database for the outwash and glacial till soils. The aspen species within this study have large ranges across North America and are shown within the United States (derived from [26]).

Field Sampling
Site selection was limited to aspen forest stands from the Wisconsin Forest Inventory and Reporting System (WisFIRS) within the 40-55-year harvest window. Sites were located on Entic Haplorthods or Typic Udipsamments soil map units from the gridded Soil Survey Geographic (gSSURGO) that are currently restricted from WTH [8]. Soil pedons were sampled cooperatively with Natural Resource Conservation Service (NRCS) soil scientists during the summers of 2016 and 2018 (Figure 1), following standard procedures, with bulk densities taken in triplicate [27]. The 2018 pedons were also sampled for leachates from 150-160 cm during the summer (12 July 2018) and resampled by bucket auger from 150-160 cm during winter (pre-frost layer, 8 January 2019) and spring (snow-on conditions, 12 March 2019) to capture the seasonal variation of water-soluble nutrients below the rooting depth. The below-rooting depth was chosen as previous research has defined the rooting depth of aspen on outwash soils at 150 cm [24]. The nearest dominant or codominant aspen tree to the soil pit was sampled for nutrient content from 18-25 June 2018 and on 11 September 2018. Each tree was felled and sampled for a basal disc, a disc at diameter at breast height (DBH), bolewood, barkwood, live branch, dead branch, twig, and leaf.

Laboratory Analysis
The soil samples were transported to Michigan Technological University, Houghton, MI. Leachate concentrations were obtained by shaking 20 g of field moist soil with 100 mL of deionized water for one hour and filtered through a Whatman 42 filter and a 0.45 µm filter [28,29]. The P, K, Ca, and Mg concentrations were determined with a Perkin Elmer Optima 7000 DV Inductively Coupled Plasma Optical Emission Spectrometer ICP-OES, and the N was obtained by Shimadzu DOC/TN analyzer. Soil pool values for the outwash soil sampled with the NRCS were obtained from the results of the USDA-NRCS NSSC National Soil Survey Laboratory. Soil pool values were also obtained from the National Cooperative Soil Survey (NCSS) database. Total N was determined by method 4H2a (dry oxidation), extractable P by method 4G4, and extractable Ca, Mg, and K were determined by method 4B1a1a (NH4OAc and 2 M KCl rinse) [30]. Total P, Ca, Mg, and K were obtained by ashing in a muffle furnace to remove carbon from the sample for 8 h at 500 • C and then acid digested by EPA 3052 method and processed by a Perkin Elmer Optima 7000 DV Inductively Coupled Plasma Optical Emission Spectrometer (ICP-OES). Missing values were imputed through kNN imputation in the Caret package [31] in the R statistical platform at the horizon level, and the weighted average by horizon depth was calculated for a depth of 150 cm.
The tree tissue samples were transported to the Wisconsin Department of Natural Resources, Forestry Headquarters, Rhinelander, WI, and dried at 55 • C to a constant mass. Tissues for N analysis were ground in a ball mill for 5 min, equipment cleaned with ethanol between samples, and 5 g were analyzed on a Costech 4010 Elemental Analyzer. For nutrients P, K, Ca, and Mg, tissues were ground in a Wiley mill, acid digested by EPA 3052 method, and processed by a Perkin Elmer Optima 7000 DV Inductively Coupled Plasma Optical Emission Spectrometer (ICP-OES). The DBH disc was used for both bolewood and barkwood samples. The barkwood was sampled from the outer to the inner bark, and bolewood was sampled from the inner bark to the center ring. The branches were sampled with a cross-section taken from the center of the branch, with care taken to ensure the cross-section remained intact. The twigs and leaves were sampled in their entirety. Roots were grouped into coarse (≥5 mm) [32] and fine roots (<5 mm) and pooled by soil genetic horizon type (A-AE, E-EA, Bhs, Bs1, Bs2, Bs3, BC, C1, C2) to obtain enough volume for nutrient analysis.

Ecosystem Inputs
The two main nutrient inputs for forest ecosystems are soil weathering and atmospheric deposition. Soil weathering has been estimated by several methods, and the depletion method is recommended for use as there is more confidence in the assumptions and it has the least uncertainty [24]. The depletion method measured the texture classes of the soils, conducted elemental analysis by X-ray diffraction, converted to elemental percent, and used multiple regression to calculate depletion [24]. The minimum, average, and maximum weathering rates for Ca, Mg, and K from the depletion method located on outwash soils [24,33] were used to create linear regression equations between elemental concentrations and sand percentages, except for the minimum K equation, which was non-linear and required exponential regression. The weathering rates of outwash soils have been related to the aeolian mass; hence the equations were necessary to adjust the weathering rates from the depletion method for Ca, Mg, and K to the average soil texture observed within the soil pits (Table A1). The texture-adjusted minimum, average, and maximum weathering rates were used as inputs into the nutrient budget for the observed soils sampled with the NRCS. The average weathering values from the outwash soils were used for the NCSS outwash soils and the Warba series for the NCSS Alfic Haplorthods [24]. The weathering rates from both the texture-adjusted and non-texture-adjusted outwash soils were used as model parameters to calculate recovery years after a harvest. Since the textural variation has been observed previously on the outwash soils [24], we treated the weathering rate as a parameter and considered the minimum, average, and maximum weathering rates per element.
Atmospheric deposition was comprised of wet and dry deposition. Wet atmospheric deposition has been measured since 1978, and total atmospheric deposition has been estimated from 2000 to the present [25]. The deposition amounts vary spatially but generally show a pattern of lower deposition in the northeastern part of the state of Wisconsin and higher in the south and southwestern portion within an example 45-year rotation length ( Figure A1). Given the differences in deposition amounts, the nutrient model was calculated for three regions. The regions were northeast, northcentral, and northwest, corresponding to the dominant advancement and retreat patterns from the last glaciation [34]. The deposition values for N, K, Mg, and Ca were extracted from the total deposition raster files [25] by the sampling points for years 2000-2017, and an average value per region was calculated per year. The deposition time series were transformed, detrended, and forecasted annually to the end of the rotation length (45 years) using the R libraries t-series [35], TSA: Time Series Analysis [36], and forecast [37,38]. The augmented Dickey-Fuller test was performed to ensure no unit root existed after transformation and detrending, which would result in a skewed forecast. The total deposition values were best detrended by differencing and log transformation ( Figure A2; Table A2). Typically, the series passed or nearly passed the adjusted Dickey-Fuller test (α = 0.05), except for the K in the northeast and northcentral areas (Table A2). However, across all nutrients, there was a spike in the 2017 values, which appeared to cause the K time series to fail since the removal of the 2017 value results in generally white-noise or detrended processes ( Figure A2). Since the transformed K series appeared to be white-noise and other transformations also did not pass the Dickey-Fuller tests, the differenced and log-transformed values were used to forecast the K time series. The annual observed values from 2000-2017 combined with the forecasted values for one rotation were used for repeated rotations.

Ecosystem Outputs
The two main nutrient losses from the forest ecosystem are leaching and harvest removals. Leaching rapidly increases after harvest and slows to the undisturbed state after approximately 5-10 years [5,39]. To approximate this post-harvest effect as a time series, the undisturbed leaching rate must first be provided by the field observed values. The undisturbed leaching rate was established by multiplying the observed seasonal nutrient concentrations by the outflow of a local water balance. The water balance was calculated monthly and was comprised of inflow as precipitation minus output from evapotranspiration. The change in water storage, as determined by the water holding capacity of the soil pedons, was used to complete the water balance [28,40]. Evapotranspiration was calculated by monthly average temperature and hours of daylight, as determined by latitude [41]. The undisturbed leaching values could then be modified by leaching ratios at 1, 2, and 5 years post-harvest [5], with the weighted average being applied between the available times. The control (CON) assumed no biomass removal and used the undisturbed leaching rate.
Decisions about harvest removals must consider differences between the nutrients removed from SOH and WTH, which determine the volume of removals by tissue. The starting forest volume available for harvest removals was calculated by scaling the sampled tissues per tree on a per hectare basis. The sampled tree tissue volumes were normalized to 45 years and scaled allometrically by tissue type; bolewood, barkwood, live branch, leaf, twig [42], dead branch [43], coarse roots [32] and stump, and large and medium roots [44]. The volume per tree was scaled to trees per hectare by the average basal area and site index of aspen stands located within the rotation length window (40-55 years) and on restricted soils in the stocking guide [45]. The available forest volume for harvest was then modified by harvest type, with SOH consisting of the merchantable volume of bolewood and barkwood [44] and WTH removing approximately 30% more fine woody debris (FWD) and 60% more coarse woody debris (CWD) [19]. The foliar and twig volumes were considered as FWD, and the live and dead branches were considered as CWD. The winter WTH was modified as it was found to leave 33% more FWD on-site than the remaining WTHs due to increased breakage [19]. Seasonal changes in tissue nutrient concentrations were applied to these volume estimates in determining how the harvest type affected changes in nutrient removal in different harvest times.
The observed tissue concentrations were multiplied by the volume of each tree component to provide the number of macronutrients removed by harvest type. The FWD and CWD concentrations were also modified by the ratios between the seasonal differences. Seasonal differences were determined by the average of the season for foliar (leaf + bud) and twig concentrations [13,[46][47][48][49]. The seasonal fluctuations and the seasonal averages across studies are presented for the foliar and twig nutrient concentrations (Appendix A Figures A2 and A3).

Nutrient Capital and Modeling Assumptions
The harvest intensity simulations considered losses to the system from leaching (OL) and harvest removals (OH) and inputs to the system from total deposition (ID) and weathering (IW) on an annual basis.
Six harvest types were considered: CON, SOH, and seasonal WTHs (Winter, Spring, Summer, Fall). The CON did not remove biomass, and the SOH scenario did not include seasonal changes because the branches were left on site. The WTH harvest type required the twig and leaf-bud values to be adjusted to the seasonal averages by the seasonal ratios and sampling time [13,[46][47][48][49] and differences in volume removed [19]. The initial soil value is the starting nutrient capital from which the inputs and outputs can be compared annually. The soil pool was defined as the top 150 cm, and soluble nutrients below this depth are considered below the primary rooting depth [24]. The unharvested volume composed of root and stump values was then combined with the soil values to constitute the nutrient capital post-harvest. The nutrient budget was calculated on a rotation length of 45 years which is a common rotation length of aspen in the region [50]. Each rotation included one harvest year, and the harvests were repeated until the capital pool reached depletion or up to 11 rotations (495 years). Confidence intervals were calculated for the harvest removals and soil capital (α = 0.05). The maximum soil and minimum harvest bounds were used to create an ecosystem upper bound, and the minimum soil and maximum harvest bounds were used to create a lower bound for the ecosystem. Although a forest stand will likely experience different types and seasons of harvest, the goal of this study was to evaluate sustainable harvest scenarios from a soil element "recovery" perspective. Recovery years were calculated by the time it took for the annual outputs and inputs of an element to replace the values lost to one harvest. The average atmospheric deposition per year across the three modeled locations was used for the deposition input, and the undisturbed leaching value was used for the output 45 years and later. The recovery years were calculated for comparison with typical rotation lengths to observe the harvesting rate under several weathering scenarios.

Ecosystem Nutrient Distribution
The soil contained the largest nutrient stocks in these forest ecosystems ( Table 1). The sampled outwash pedons held a larger amount of K and lower amounts of Ca and Mg than the NCSS pedons ( Table 1). The NCSS Alfic Haplorthods pedons, by comparison to the outwash soils, held much larger P, Ca, and Mg, similar N, and smaller K values ( Table 1). The vegetation contained the second largest ecosystem nutrient values ( Table 1). The vegetation on outwash soils contained a large proportion of ecosystem K, Ca, and Mg, whereas only the vegetation K values on the Alfic Haplorthods contained a large percentage of this cation ( Table 1). The N and P within the vegetation pools were small in comparison to outwash and Alfic Haplorthods soils ( Table 1). The largest pools of nutrients within the vegetation tended to be in bolewood, barkwood, leaf, live branches, and stump combined with large and medium roots across all reported nutrients. In comparison, the coarse roots, dead branches, and twigs tended to harbor much less of the overall nutrients within the ecosystem (Table 1). Foliar tissues showed a spike in N, P, and K during the growing season, whereas Ca and Mg showed a peak during the dormant season ( Figure A3). Twig nutrient concentrations generally decreased during the growing season and had larger concentrations during the dormant season ( Figure A4). Winter P values were not available, and the twig concentrations are more limited than the foliar, particularly in the winter, where only one study presented values.

Modeled Harvest Removals
Not surprisingly, WTH removed a larger portion of all nutrients than the SOH due to more volume being removed (Table 2, Figure 2a). The winter season removed the lowest amounts of macronutrients for WTH ( Table 2). The remaining seasons had different impacts by element, with the spring harvest removing larger amounts of N and K, the summer harvest removing larger amounts of K and Ca, and the fall harvest removing the largest amount of Ca ( Table 2). The P and Mg pools did not differ seasonally among WTHs, but the WTH removals are larger than the SOH ( Table 2). The harvest removals were one or two orders of magnitude larger than the leaching outputs depending on the nutrient, showing the degree to which the harvest removals drive the scenario outputs ( Table 2). The inputs did not vary greatly since the deposition values for K and Mg are similar to the average weathering rates across regions, but the Ca deposition values were larger than the weathering rates ( Table 2). The maximum weathering rates were larger than deposition inputs across all elements ( Table 2). outwash soils, the percentage of removal was approximately 5% or less of N and P, and between 20-35% of K, Ca, and Mg (Figure 2b). The harvest mass for the Alfic Haplorthods was approximately 5% or less for N, P, and Mg, and 10% or less for Ca, but remained within 30-40% for ecosystem K. In comparison to SOH, WTH removed an additional 5-10% of the ecosystem K, Ca, and Mg, and 1-2% of the N and P in the outwash ecosystems. On the Alfic Haplorthods ecosystems, the WTH removed an additional 1-2% N, P, Ca, and Mg and 7-10% K compared to SOH. The Alfic Haplorthods showed a much smaller percent of Ca and Mg per harvest, but K remained large. N and P constituted less than ten percent of the soil pools per harvest.

Harvest Time Series-Sampled Outwash Pedons
Although the deposition values varied spatially, the harvest removals drove the trend in the time series models. The time series from the northeastern model for the sampled NRCS outwash soils under the texture-adjusted average weathering scenario showed minimal impacts to N and P but rapid depletion of K, Ca, and Mg ( Figure 3). The sampled outwash pedons N and P pools showed negligible harvest impacts (Figure 3a,b), and the exchangeable K and Ca was depleted on average in 4 WTHs and 5 SOHs (Figure  3c,d) and Mg depleted in 5 WTHs and 6 SOHs (Figure 3e). In comparison to the harvests, the unharvested scenario (CON) showed gradual losses in P and K (Figure 3b,c) and accumulation in N, Ca, and Mg (Figure 3a,d,e). The total soil pool size was modeled to decrease through time as the weathering occurs (Figure 3a-e). The Alfic Haplorthods showed a much smaller percent of Ca and Mg per harvest, but K remained large. N and P constituted less than ten percent of the soil pools per harvest.
When the harvest mass was compared to the remaining ecosystem nutrients for the outwash soils, the percentage of removal was approximately 5% or less of N and P, and between 20-35% of K, Ca, and Mg ( Figure 2b). The harvest mass for the Alfic Haplorthods was approximately 5% or less for N, P, and Mg, and 10% or less for Ca, but remained within 30-40% for ecosystem K. In comparison to SOH, WTH removed an additional 5-10% of the ecosystem K, Ca, and Mg, and 1-2% of the N and P in the outwash ecosystems. On the Alfic Haplorthods ecosystems, the WTH removed an additional 1-2% N, P, Ca, and Mg and 7-10% K compared to SOH.

Harvest Time Series-Sampled Outwash Pedons
Although the deposition values varied spatially, the harvest removals drove the trend in the time series models. The time series from the northeastern model for the sampled NRCS outwash soils under the texture-adjusted average weathering scenario showed minimal impacts to N and P but rapid depletion of K, Ca, and Mg ( Figure 3). The sampled outwash pedons N and P pools showed negligible harvest impacts (Figure 3a,b), and the exchangeable K and Ca was depleted on average in 4 WTHs and 5 SOHs (Figure 3c,d) and Mg depleted in 5 WTHs and 6 SOHs (Figure 3e). In comparison to the harvests, the unharvested scenario (CON) showed gradual losses in P and K (Figure 3b,c) and accumulation in N, Ca, and Mg (Figure 3a,d,e). The total soil pool size was modeled to decrease through time as the weathering occurs (Figure 3a-e).

Harvest Time Series-NCSS Pedons
The seasonal time series from the northeastern model for the NCSS outwash soils under the average weathering scenario (Figure 4) and the NCSS Alfic Haplorthods ( Figure 5) showed different impacts by site. The NCSS outwash pedons showed minimal harvest impacts on N and P, although differences between SOH and WTH appeared after 5 rotations. The average remaining element pools were depleted between 3-10 rotations, varying by harvest treatment: K (rotations: 4 SOH; 3 WTH), Ca (rotations: 6 SOH; 5 WTH), Mg (rotations: 10 WTH) (Figure 4). The NCSS Alfic Haplorthods pedons showed only the average K was depleted (rotations: 4 SOH; 3 WTH), and N, P, Ca, and Mg manifested differences between SOH and WTH halfway through the model. Seasonal differences appeared but rarely changed the number of total harvests and only near the end of the model.
In addition to the time series, the recovery years also showed the harvest outputs were greater than inputs at the current rotation length under all scenarios and soil types ( Table 3). The K, Ca, and Mg recovery years exceed the average rotation length by a factor of two to four under the maximum weathering values on the non-texture-adjusted soils [24]. When the weathering rates were adjusted to textures observed in the field, the recovery years increased greatly and further demonstrated the non-sustainable removal rates for the outwash soils ( Table 3). The Alfic Haplorthods showed recovery years approximately double the average rotation length for SOH of Ca, SOH, and WTH of Mg, and triple of WTH for Ca (Table 3). SOH and WTH removals of K were not replaced through weathering and deposition under all scenarios except for the maximum outwash weathering rate, and even under that scenario vastly exceeding an average rotation length. 1 Figure 3. Simulated soil nutrients by scenario for the sampled NRCS pedons with texture-adjusted weathering rates. The N and P showed lower harvest impacts (a) N are totals (exchangeable unavailable) and (b) P total is unchanged (weathering unavailable). The depletion of K and Ca (c,d) is reached by 5 SOHs and 4 WTHs, and the Mg (e) is depleted in 6 SOHs and 5 WTHs. The control (CON) shows (a,d,e) accumulation in N, Ca, and Mg (b,c) gradual losses in P and K. Confidence intervals show large variability (α = 0.05).     Table 3. Recovery time in years for weathering and deposition inputs to replace outputs of a single harvest and leaching by harvest type. The number of recovery years demonstrates that the inputs do not replace outputs at current rotation lengths, besides K and Mg in the non-texture-adjusted ceiling and average of the maximum weathering rates. * indicates the weathering rates and deposition do not exceed the outputs indicating that the harvest nutrient pool is not replaced.

Harvest Removals
The annual inputs of K, Ca, and Mg to the aspen ecosystems over 45 years did not replace the nutrients lost to harvest. On the other hand, pools of nutrients associated with organic matter (mainly N and P) were minimally impacted by WTH in this study, which is consistent with extensive modeling efforts conducted in the same ecosystem [51]. This current study highlights the importance of base cations present in biomass removal. Even under maximum rates of base cation release from mineral weathering, WTH shortened the total number of rotations until depletion for the K and Ca (4 WTHs; 5 SOHs) and Mg pools (6 WTHs; 7 SOHs). The larger context of nutrient-poor soils represented by the NCSS restricted pedons under maximum weathering rates showed depletion for K (4 WTHs; 5 SOHs) and Ca (5 WTHs; 7 SOHs). Ca and Mg pools have long been reported as a potential concern using WTH in aspen stands [4,5] and especially on coarsetextured soils [7]. Calcium translocation from the deeper horizons has been proposed as an unmeasured input [52]. Although there are differences in bedrock types in these studies, if the relatively Ca-rich values from prior work [52] were incorporated into the harvest models in this study, the recovery years would be 47-135 for SOH and 60-173 for WTH, on soils 10 m from bedrock. This would indicate that only the highest values from the shallowest bedrock would approach satisfying removals on an average rotation length. In the absence of increased nutrient inputs to the soils, the recovery ages under the maximum weathering rates exceed the rotation lengths of aspen.
The harvest simulations call into question the practice of continual SOHs and WTHs on the sandy outwash soils in the absence of increasing the input terms of K, Ca, Mg, and also on the Alfic Haplorthods without increased inputs of K. Our approach in this assessment was conservative in two ways: (i) outputs to the system through harvest removals and leaching were smaller than previously reported (cf., [53]; Table A3), and (ii) inputs through weathering and deposition were parametrized across the range of possible values, and an additional input term of nutrient translocation was provided. There are some reasons for inconsistencies across studies evaluating base cation balances in response to harvesting. Intensively managed forests can be further along in the harvest rotation models or have lower starting soil pool values than the average field observed soil values (for example, in this study, 989 K, 3217 Ca, and 336 Mg kg ha −1 ). For instance, if the 30 cm pools from the observed WTH and SOH are extrapolated to 150 cm, the average exchangeable soil nutrients would be SOH (3650 K, 2100 Ca, and 1850 Mg kg ha −1 ) and WTH (3650 K, 550 Ca, and 1650 Mg kg ha −1 ) ( [10]; Table A3). Extrapolating the top 30 cm values to 150 cm will overestimate the soil macronutrient values since the 30 cm values are comprised of the A-B horizons, which have larger macronutrient pools than the lower soil horizons. The average Ca value reported after WTH has important consequences because the next harvest removal, regardless of the type, has the potential to deplete Ca. These results indicate that some outwash soils may be more at risk for negative harvest impacts depending on the management history and weathering status of the aeolian soil mass.

Long-Term Nutrient Budgets Comparison to Short-Term Harvest Intensity Experiments
The goal of this study was to model long-term patterns in the ecosystem and may not depict short-term phenomena well due to the absence of decomposition models. Although the outputs of the ecosystems through both harvest and leachate are smaller in this study than previously reported, the harvest time series show major losses to the outwash soils through SOH and WTH and shorter recovery years than rotation lengths for both soil types. In nutrient budget studies, WTH has been identified as having the potential to remove two to three times the nutrients as SOH [6], with larger impacts on coarse-textured soils [7], and Ca and Mg identified as at risk from WTH on more productive soils [4,5]. The short-term (≤1 rotation) observations of the soil pools post-harvest found non-significant or site-specific impacts between WTH and SOH [10,11]. However, the models did show patterns of larger WTH removals than SOH that decreased available nutrients, but significant differences may not be observed until after multiple sequential WTHs (Figures 3-5). The long-term effects of harvest intensity have been suspected since mature aspen has accumulated the most nutrients [17,54], and our models confirmed that pattern.
Complications of detecting differences in the exchangeable soil pools between SOH and WTH in the short term include seasonal water table fluctuations, nutrient translocation to foliage [10], and root decomposition [11]. The harvest confidence intervals show that significant differences may not be detected until multiple sequential harvests have occurred, and the variation within the ecosystem harvest intervals is large due mainly to the soil values. . Barring having any pre-harvest measurements of soil nutrient pools, differences in harvest type appear difficult to detect in soils after a single harvest, although differences in removal have been observed in harvest residues [19]. These data indicate that the detection of soil nutrient depletion in response to harvest systems may take several rotations, depending on soil type or harvest intensity.

Modeled Ecosystem Values
The harvest removals reported herein were consistent with previously reported aspen harvest values ( [5]; Table A3) and slightly above the reported harvest at 30-year estimation ( [4]; Table A3). In comparison to the most similar outwash soils, both the sampled and NCSS outwash pedons fall between extrapolated Entic Haplorthods [10] and the Typic Udipsamments [53], except for the Mg, which fell below both, but remained near a Typic Haplorthods ( [55]; Table A3). The differences in soil nutrient values within the Entic Haplorthods and Typic Udisamment studies are likely due to the textural variation within the outwash soils and the depletion status of the aeolian mass. Given the differences in soil type and potential overextrapolation of the literature soil values due to differences in sampling depths, the soil values reported in this study fall within historical values (Table A3) for the glacial till and outwash soils. Additionally, the leaching values from the water budget approach are more conservative than have been previously reported ( [53]; Table A3). The leaching rates measured at 60 cm [53] and 40 and 130 cm [5] may have over-estimated leaching losses in prior studies because these depths are still within the zone of active root uptake (150 cm) [24].

Desired Future Condition
It is important to remember many of the soils were not agriculturally productive and considered the "lands nobody wanted" [1], and exist now as publicly owned forest lands. The low fertility of the sandy soils in combination with the modeled outputs exceeding inputs at the current rotation length indicated that future yields could be diminished if the soil becomes depleted. The impacts from the increased removal of nutrients should manifest in decreased productivity when a single nutrient reaches the limiting level [56]. Declines in productivity have been attributed to WTH and removal of the forest floor after 20 years in coarse-textured soil [11], and yet differences in forest productivity were not observed after one WTH rotation [18]. These results indicate, in some locations, nutrients have not become limiting to the aspen, which is adapted to nutrient-poor soils. If such a point of depletion were reached, for example, after 3-6 rotations of intensive harvesting with no increases to the input terms, the sandy outwash ecosystem could be expected to convert to pine or oak savannahs and function more similarly to a barrens ecosystem. If the cations became limiting, the yields could be expected to decrease and reduce the harvest output terms sequentially. The modeling indicates that harvest removals outpace the inputs to the system. Since the recovery years are longer than the biological and economic rotation of aspen, a clear answer is to increase the inputs if harvest removals keep pace with current levels [6]. Together these analyses indicate that current harvest levels are not sustainable if the desired future condition of these forests is aspen. These findings can be used in evaluating the desired future condition of intensively harvested aspen forests (in transitioning to oak or pine barrens ecosystems).

Conclusions
Under all weathering scenarios, the recovery years from a SOH or WTH exceeded the average rotation length (45 years) of the aspen. Simulations showed that increased harvest intensity from SOH to WTH of aspen on nutrient-poor soils shortened the available nutrient pools on average by one full rotation length for Ca, Mg, and K, and K for the slightly richer Alfic Haplorthods. N and P were less impacted for both soils since the percent removed in each harvest was generally small in comparison to the residual tissues and soil pools. The harvest simulation models called into question the sustainability of the SOH practice on nutrient-poor soils given the impacts to Ca, Mg, and K, and for K in the Alfic Haplorthods. The seasonal differences between WTH removals were small in comparison to the difference between WTH and SOH. The large portion of base cations located within the vegetative pools implies that inputs by weathering and deposition could not keep pace with harvest removals within the sandy, nutrient-poor soils in northern Wisconsin. Although there was much uncertainty due to variation in aeolian mass, weathering depletion status, and harvest removals, our models indicate that the NRCS sampled soils would be depleted in K, Ca, and Mg in four or five WTHs and five or six SOHs. The NCSS outwash soils were depleted within three (K), five (Ca), and ten (Mg) WTHs and four (K) to six (Ca) SOHs, and the NCSS Alfic Haplorthods were depleted of K within three WTHs and four SOHs. The long-term models bridged the gap between short and long-term studies and showed the harvest intensity effects might not be apparent until multiple harvests have occurred.  Figure A1. Total atmospheric deposition (wet + dry) throughout one 45-year rotation (adapted from [25]) of nutrients across the state of Wisconsin. In general, the south and southwestern areas are receiving two to three times the amount in the northeastern area. Figure A1. Total atmospheric deposition (wet + dry) throughout one 45-year rotation (adapted from [25]) of nutrients across the state of Wisconsin. In general, the south and southwestern areas are receiving two to three times the amount in the northeastern area. Forests 2021, 12, x FOR PEER REVIEW 19 of 23 Figure A2. Time series of total atmospheric deposition time series for regional harvest model inputs after differencing and log transformation. When excluding the large spike in the 2017 deposition, seen as the second from right value due to differencing, the processes appear to generally be detrended or white-noise. Figure A2. Time series of total atmospheric deposition time series for regional harvest model inputs after differencing and log transformation. When excluding the large spike in the 2017 deposition, seen as the second from right value due to differencing, the processes appear to generally be detrended or white-noise.    Figure A4. Seasonal twig macro-nutrient concentrations by day of the year observed in aspen, as well as the average of the seasonal averages across studies (Derived from [13,[46][47][48][49]). The macro-nutrients all show increases during the dormant Figure A4. Seasonal twig macro-nutrient concentrations by day of the year observed in aspen, as well as the average of the seasonal averages across studies (Derived from [13,[46][47][48][49]). The macro-nutrients all show increases during the dormant season while decreasing during the growing season. Trembling aspen (Populus tremuloides) is represented by POTR, and bigtooth aspen (Populus grandidentata) is represented by POGR.