Integrated Stochastic Life Cycle Assessment and Techno-Economic Analysis for Shrub Willow Production in the Northeastern United States

: The refereed literature contains few studies that analyze life cycle assessment (LCA) and techno-economic analysis (TEA) methodologies together for lignocellulosic bioenergy systems, using a stochastic modeling approach. This study seeks to address this gap by developing an integrated framework to quantify the environmental and ﬁnancial impacts of producing and delivering shrub willow in the Northeastern United States. This study analyses four different scenarios from a combination of two different initial land cover types (grassland, cropland) prior to willow establishment, and two harvesting conditions (leaf-on, leaf-off). Monte Carlo simulations were performed to quantify the uncertainty of the results based on a range of ﬁnancial, logistical, and biophysical variable input parameters (e.g., land rental rates, transportation distance, biomass yield, etc.). Growing willow biomass on croplands resulted in net negative GHG emissions for both leaf on and leaf off scenarios for the baseline. The GHG emissions were lowest for the leaf-off harvest on cropland ( − 172.50 kg CO 2eq Mg − 1 ); this scenario also had the lowest MSP ($76.41 Mg − 1 ). The baseline grassland scenario with leaf-on harvest, results in the highest net GHG emissions (44.83 kg CO 2eq Mg − 1 ) and greatest MSP ($92.97 Mg − 1 ). The results of this analysis provide the bioenergy ﬁeld and other interested stakeholders with both environmental and ﬁnancial trade-offs of willow biomass to permit informed decisions about the future expansion of willow ﬁelds in the landscape, which have the potential to contribute to GHG reduction targets and conversion into fuels, energy, or bioproducts for carbon sequestration and ﬁnancial beneﬁts.


Introduction
In 2019, U.S. energy consumption consisted mainly of petroleum and natural gas at 36.7% and 32.1% of all energy sources consumed, respectively [1]. There has been an increased focus on reducing this high domestic fossil fuel consumption, and increasing consumption from renewable energy, as more attention is paid to the climate change impacts of greenhouse gas (GHG) emissions from fossil fuels. In 2020, the United States reached a new record in terms of renewable energy consumption, at 12% of total energy consumption. Of this total renewable energy consumption, approximately 39% came from biomass sources. Specifically, 18% was attributed to wood, while another 17% was attributed to biofuels [2]. An additional 4% came from biomass waste.
In terms of renewable energy consumption, using biomass, such as willow shortrotation coppice (SRC), has several benefits. It has a short harvest cycle of three years; can be grown on marginal lands; and has net GHG emissions that are negative or very low over its life cycle, depending on previous land use and transportation distances for the biomass [3,4]. Willow biomass grows particularly well in the northern, temperate regions of the United States, and can be used as a replacement for fossil fuels such as natural gas, or fuel oil for heating and other energy purposes [5]. To understand the impacts of replacing fossil fuel derived energy sources with renewable energy sources such as willow biomass, it is vital to determine the environmental and financial impacts of doing so. While many studies quantify the environmental [4] or financial impacts [6] of shrub willow production, there have not been previously published efforts to use an integrated life cycle assessment (LCA) and techno-economic analysis (TEA) approach. Unlike relying on a separate LCA or TEA methodology, an integrated approach measures both the environmental and financial impacts associated with certain feedstocks, technologies, or energy pathways by using consistent system boundaries, input data, and assumptions.
There are a limited number of refereed studies pertaining to lignocellulosic bioenergy systems that report both the LCA and TEA methodologies together (Table 1), and these studies face challenges in modeling complexities and data limitations. Consequently, these studies differ in their environmental and financial impacts reports, system boundaries, functional unit, and other key aspects (Table 1) [7][8][9][10][11][12][13][14][15]. Data availability and similar limitations often represent a complication for LCA and TEA of bioenergy system studies [16]. Such limitations can be associated with the lack of production-scale data, process equipment cost data, and dependable regional data. To overcome these challenges, experimental data can be used, and process equipment cost data can be obtained from several vendors [16]. Another challenge pertains to the incorporation of uncertainty in LCA and TEA studies [16]. Integrated bioenergy LCA and TEA studies are comprehensive, and can include multiple shared input parameters that exhibit substantial uncertainty.
Several bioenergy studies have utilized a comprehensive sensitivity analysis to handle uncertainty. One study calculated the GHG emissions and minimum selling price (MSP) of biofuel, biopower, and pellets produced from switchgrass, shrub willow, and Miscanthus grown in the Northeastern United States [12]. To account for uncertainty, a sensitivity analysis was performed on the yield, conversion rate, transportation distance, and internal rate of return (IRR). Results of this study concluded that biofuel production yielded the highest GHG emissions, and pellet production was responsible for the lowest GHG emissions along the life cycle. Additionally, the MSP range for bioenergy products was 7.7 to 47.9 $ GJ −1 . The lowest MSP was achieved by pellet fuel, and biopower production was responsible for the highest MSP. Another study used LCA and TEA methodology to model process designs that integrated ethanol and enzyme production, using spruce production logging residues for the feedstock [14]. Three sensitivity analyses were conducted and accounted for GHG emissions from ethanol production, primary energy balance, and MSP. This analysis concluded that it was possible to decrease GHG emissions from ethanol production through an integrated production process, and that the cost of obtaining enzymes (manufactured off-site) may be similar to the cost of the integrated process. Another recent study analyzed the GHG emission impacts and NPVs of pellet, biopower and biofuel from a variety of lignocellulosic feedstocks [17]. The lignocellulosic bioenergy LCA and TEA studies summarized in Table 1 differ in  their environmental and economic impacts modeled, system boundary, functional unit,  conversion technology, end product, feedstock, and region. The following three areas are  key to providing reliable LCA and TEA analysis results for bioenergy systems: integrated analysis, dependable data, and comprehensive uncertainty analysis. To the authors' knowledge, these three aspects have not been applied together in a single peer-reviewed study for shrub willow production. This study bridges this research gap by providing the first integrated cradle-to-gate LCA and TEA for shrub willow production in the Northeastern US. To accomplish this, an integrated analysis is performed with shared input parameters and reliable datasets. Additionally, feedstock (moisture content, root-to-shoot ratio, yield), soil quality (soil organic carbon), logistics (transportation distance), and financial (planting costs, diesel price, land rental rates) uncertainty is accounted for by using best fit analysis to develop probability distributions for these stochastic input parameters. Probability distributions are also used to represent the MSP values for each scenario, and reflect a range of potential prices. This study presents a comprehensive, integrated methodology to quantify the environmental and financial impacts of growing shrub willow.
An increasing number of renewable energy studies are utilizing stochastic methodology in their models [19]. Stochastic methodology is used to model uncertainty and often utilizes Monte Carlo Simulation (MCS) methods [19]. While this method is commonly employed in studies that assess either the life cycle environmental [20] or financial impacts [21] of biomass production and/or conversion systems, it is rarely employed in studies that report both impacts. This can be due to the already complex nature of analyses involving both LCA and TEA methodologies. Still, handling uncertainty in an integrated LCA and TEA would provide more impactful results to industry stakeholders and policymakers.
The objectives of this research were to: (1) Determine the GHG emissions, and minimum selling price (MSP) of shrub willow biomass production in the Northeastern U.S. under a combination of two land-use changes and two harvest seasons; (2) Assess the uncertainty associated with the GHG emissions and MSP based on variations of key parameters used in the integrated LCA and TEA. These variable parameters encompass willow biomass yield, moisture content, transportation distance from the field to the gate of a conversion facility, soil organic carbon change, diesel price, land rental price, planting price, and root-to-shoot ratio. The next section details the integrated LCA and TEA framework utilized, along with the system dynamics of the willow supply chain, and details regarding the individual LCA and TEA modeling. Section 3 depicts both the results and a discussion of the life cycle GHG emissions and MSPs, along with a discussion regarding carbon abatement cost (CAC) values.

Integrated LCA and TEA Framework of Willow Production
An integrated LCA and TEA approach was used to quantify the baseline climate change impact and MSP of a typical willow biomass production system in the Northeast United States, over a 24-year life cycle. The framework for this analysis (Figure 1) was created so the LCA model and TEA model utilize the same scenarios, system boundaries, and input parameters. The four modeled scenarios account for a combination of two landuse change scenarios, and two harvesting seasons: 1. Grassland + leaf-off (S1. GL + LOFF); 2. Grassland + leaf-on (S2. GL + LON); 3. Cropland + leaf-off (S3. CL + LOFF); and 4. Cropland + leaf-on (S4. CL + LON). The choice of these options was informed based on prior research findings [4,[22][23][24] and recent changes in practices among willow growers. Two types of direct land-use change were modeled due to their differences in soil organic carbon (SOC) changes, with SOC increasing in cropland and decreasing in grassland through conversion to willow cultivation [4,22,23]. Leaf-on harvest and leaf-off harvests were studied as separate scenarios, as recent studies have shown differences in harvester throughput and fuel consumption [24]. The LCA focuses on climate change impacts to produce one oven-dry Mg of willow chips, and follows the International Organization for Standardization (ISO) 14,040 and ISO 14,044 standards. The cradle-to-gate system boundary includes processes from the nursery, site preparation, planting, and crop management, harvesting on three-year rotations, transportation to the end-user, and removal of the willow stools after the 7th harvest ( Table 2). The same system boundary is used for the TEA. The financial and environmental impacts occurring after transportation to the end user are outside the scope of this study.
Environmental impacts originating from diesel refining, equipment manufacturing, and fuel combustion in the harvester and other equipment are included in the analysis. Direct land-use change impacts are considered through the incorporation of changes in SOC for grassland scenarios (S1 and S2), and cropland scenarios (S3 and S4). Certain activities associated with field preparation such as weed control (first growing season), plowing and tilling, planting of cover crops, and rolling were necessary only for the conversion of grasslands to willow. These steps were not needed for croplands because these fields were already in a cultivation cycle. Some activities or processes were unique to the LCA, while others were only applicable to the TEA. For example, carbon sequestration in roots, decomposition of litterfall, and change of soil organic carbon were accounted for only in the LCA model, and were not included in the determination of the MSP of willow biomass, as they have no current financial value. Similarly, because land rental rate was considered a financial parameter, it was not included in the LCA model.
The TEA updates the original EcoWillow 3.0 S model [6] to EcoWillow 4.0 S to calculate the baseline MSPs for shrub willow production in the Northeastern United States. The updated EcoWillow 4.0 S model incorporates new stochastic variables of moisture content, transportation (field to end-user), and New York State land rental rates, in addition to the original stochastic variables of willow yield, planting costs, and farm (base) diesel price. While the TEA and LCA models are run on separate platforms (Excel and Python, respectively), their biomass production processes and stochastic values for the uncertainty analysis are the same. All scenarios for both the LCA and TEA analysis incorporated the same probability distributions using a 95% confidence interval, and were used as input variables for the stochastic treatment of average annual willow yield (dry), moisture content (wet basis), and transportation distance (field to end-user) ( Table 3). Root-to-shoot ratio and SOC change are only included in the LCA, and financial variable parameters are incorporated into the TEA model.

System Dynamics of Willow Biomass Supply Chain
The willow biomass system was modeled using data collected from field trials and commercial-scale operations pertaining to the nursery, site preparation, planting, site management, harvesting operations, and transportation of the biomass to the end-user [6,20,22,24]. At the beginning of a growing season, willow was planted in a double-row system at a density of 13,500 plants ha −1 . It was then coppiced after the first growing season to foster the development of multiple stems per shoot and better yields. Willow was harvested on three-year rotations, with fertilizer applied at a rate of 100 kg N ha −1 at the beginning of each rotation.
A total of 210,778 hectares of land suitable for growing willow was identified in five counties (Jefferson, Lewis, Oneida, St. Lawrence, and Oswego) of New York State, by applying filtering criteria based on land cover classes, land use type, slope, hydrography, and property tax parcel in ArcGIS [4]. The mean transportation distance (one way) from these identified parcels to existing end-user facilities was 57.87 km for cropland, and 58.73 km for grassland scenarios. Willow biomass yield was estimated for each parcel, based on the National Commodity Crop Productivity Index (NCCPI), which considered soil, landscape, and climate factors [4,31]. Biomass yields from croplands and grasslands were considered separately. Cropland parcels produced a slightly higher average yield (11.65 Mg ha −1 yr −1 ) than grassland (10.76 Mg ha −1 yr −1 ), because of the generally higher NCCPI scores.
Material capacity (harvester throughput) and harvester fuel consumption were modeled for green standing biomass (delivered) for dry-weather conditions for leaf-on and leaf-off harvests [24]. These models were developed from nearly 700 wagon loads of chips from commercial-scale harvests, utilizing New Holland FR9080 and FR9090 forage harvesters, and a FB130 woody crop header. Both standing biomass and harvester efficiency were used to calculate the material capacity. The harvester fuel consumption (per hour) was impacted by fuel consumption (per hectare) and harvest rate. Please see Eisenbies et al., (2020) [24] for details regarding the equations used for LOFF and LON scenarios. Table 4 depicts the mean values obtained, and their respective units, for the harvester operation across all scenarios.

LCA Modeling
The primary focus of this LCA was to determine the life cycle GHG emissions of willow biomass production in the Northeastern United States. The United States Environmental Protection Agency (EPA) TRACI 2.1 method was used to calculate the life cycle impact. The LCA of the willow biomass system was modeled in Python 2.7 by linking equations with the variable input parameters and the life cycle impacts of materials, chemicals, and processes (e.g., fertilizer, diesel) that were calculated in SimaPro (version 8.2) (Amersfoort, The Netherlands).
Willow captures carbon dioxide from the atmosphere to sequester carbon in various components of the plant, such as leaves, stems, stool, coarse roots, and fine roots. A conservative approach was applied that considered only the long-term sequestration of carbon in the stool and coarse roots. The carbon stored in the harvested willow chips (shoot) was not included in the calculation of the carbon sequestration process, because the carbon will be released back to the atmosphere once the chips are delivered to the end-user and converted to an energy product. Carbon stored in fine roots and leaves was not included in the carbon sequestration because they have a fast turnover, some portion of which will contribute to SOC, with the rest returning back to the atmosphere [32,33]. A root-to-shoot ratio was used to determine the mass of coarse roots and stool from the standing biomass, which is the above-ground harvestable biomass. At the end of the third rotation, a rootto-shoot ratio of 0.67 indicates that, given a standing biomass of 35 Mg ha −1 (dry basis) after three years of growth, the estimated below ground biomass was 23 Mg ha −1 . The amount of root biomass was expected to level off, or only slightly increase, after the third rotation [34]. Thus, the most recent root-to-shoot data collected on the third rotation of the willow biomass system was used to estimate the amount of coarse root for the system [26].
SOC data for 30 cm soil depth resulting from direct land-use change of cropland or grassland to willow, across five counties, were informed from the literature [4,23,27]. We used the 30 cm depth values rather than 100 cm, because the marginal land where willow was grown in the region had soil depths that are often restricted by hardpans or perched water tables. The SOC model indicated that SOC decreased when grasslands were converted to willow fields, while the conversion of croplands to willow fields resulted in increased SOC [23,27].
Emissions of nitrous oxide from leaf decomposition were determined based on the International Panel on Climate Change (IPCC) 2006 standards [35]. Based on the leaf biomass and nitrogen content of willow, 1% of the leaf biomass is assumed to be released as nitrous oxide. For N fertilizer, 1% of the applied N was assumed to be released as nitrous oxide (i.e., 100 kg N applied resulted in 1 kg nitrous oxide being released). Nitrous oxide emissions were then converted into carbon dioxide equivalents using 100-year global warming potentials [36].

TEA Modeling: EcoWillow 4.0 S
The original EcoWillow 3.0 S model [6] was updated to the EcoWillow 4.0 S model. Similar to the original model, monthly cash flows were generated as a function of both deterministic and stochastic parameters to derive a 24-year (288-month) baseline MSP for shrub willow production. All site tasks, from nursery and site preparation, to transportation (field to end-user) were included in the model (Table 3). These were represented using seven, three-year rotations and accounted for all costs through the crop removal phase. Diesel fuel and planting cost data were from the original EcoWillow 3.0 S study [6]. The following aspects of our study were also kept the same: inflation rate (2%); New York State diesel tax ($0.02 L −1 ); and federal diesel tax ($0.06 L −1 ).
The EcoWillow 4.0 S model made five key modifications to the original EcoWillow 3.0 S model. First, land rental rates for willow biomass production on cropland and grasslands were included. The data came from confidential 2014 surveys from the United States Department of Agriculture's National Agricultural Statistics Service [30]. This analysis was based on 10,000 simulated data points that were derived from the original USDA datasets, to maintain their confidentiality. The original datasets included 10 random and unknown New York State counties. Downstate NY counties were excluded due to high rental rates and limited land area, which made these counties unsuitable for growing biomass for bioenergy systems, and highly unlikely to achieve positive financial returns. Second, additional site preparation tasks required by the grassland scenarios to prepare the site prior to crop establishment, and to manage potential soil erosion, were included (Table 3). Third, separate equations that reflect harvesting performance for leaf-on and leaf-off conditions were used, based on hundreds of monitored loads of harvested willow chips [24]. Fourth, the timing of the leaf-on scenario was adjusted to include the earlier harvesting and transportation site tasks (Table 3). Lastly, EcoWillow 4.0 S treated the following three additional input variables in a stochastic manner: moisture content, transportation distance, and land rental rates. Best fit analysis was used to develop probability distribution functions for these inputs. An updated field trial dataset containing 9718 tax parcels (88.7% were cropland/pasture and 11.3% were grassland) that were identified as suitable for willow biomass crops across the following five New York State counties (Jefferson, Lewis, Oneida, Oswego, and St. Lawrence) was used to develop the probability distributions for both the transportation distance and willow yield [4]. The moisture content dataset included moisture content from 205 samples taken from wagon loads of chips, during commercial-scale harvesting operations at two locations in upstate NY [25].

Monte Carlo Simulation for the LCA and TEA
The uncertainty was incorporated into the analysis via Monte Carlo simulation. The Monte Carlo analysis is conducted by selecting random values from an assigned probability distribution for each variable input parameter. This study applied the same approach used in [6]. Ten thousand modeling runs are conducted by Crystal Ball to develop the output probability distributions, and the software is set to run at normal speed and to stop on calculation errors, or when a 95% confidence interval precision control limit is reached [6]. For each run, the selected input parameter values are fed to the LCA and the TEA models. The list of the variable parameters and their probability distributions can be found in Table 2.

LCA Results
The life cycle GHG emissions are negative when willow is grown on cropland, and slightly positive when grown on grassland. In both land use scenarios, leaf-on harvesting results in higher GHG emissions. The baseline life cycle GHG emissions of willow biomass crops grown on cropland, from cradle-to-gate were −172.50 kg CO 2eq Mg −1 when harvests occurred after the leaves fell (S3. CL + LOFF), and −161.35 kg CO 2eq Mg −1 for leaf-on harvests (S4. CL + LON). A higher heating value of 18.9 MJ kg −1 corresponded with a −9.13 g CO 2eq MJ −1 impact for leaf-off harvest, and −8.54 g CO 2eq MJ −1 for leaf-on harvest for the cropland scenarios [37]. This means that carbon can be sequestered in this system, and sequestration may also be possible as this woody biomass is converted into bioenergy, biofuels, and/or bioproducts, depending on the added impacts from extending the system boundary. The baseline life cycle GHG emissions were slightly positive for willow grown on previous grasslands, with 33.20 kg CO 2eq Mg −1 for the leaf-off harvest system (S1. GL + LOFF), and 44.83 kg CO 2eq Mg −1 for leaf-on harvest system (S2. GL + LON) ( Figure 2). Harvesting during the dormant season when leaves are off the plant resulted in lower GHG emissions, because harvester throughput and ground speeds were higher, resulting in lower fuel consumption per hectare or per Mg of biomass [24]. Although the GHG emissions difference between leaf-on and leaf-off harvest systems was small, it reinforced the best management practice of harvesting after the leaves fall, which also results in the most vigorous regrowth of the plants, increased quality of the chips, and greater amounts of nutrients retained as litter cover [24,[38][39][40]. Nonetheless, willow growers may strategically decide to harvest with leaf-on because of the need for year-round supply of biomass and the challenges associated with harvesting on marginal lands, which in this region is often poorly drained and no longer freezes consistently each winter.
The net negative baseline life cycle GHG emissions for the cropland scenarios (S3 and S4) were mainly due to the long-term storage of carbon in the stool and coarse roots, and the carbon sequestration in the soil (SOC). The root and stool sequestered 160 kg CO 2eq Mg −1 , which is equivalent to 1.86 Mg CO 2eq ha −1 yr −1 (for both S3. CL + LOFF and S4. CL + LON). After carbon sequestration in the stool and coarse roots, in absolute values, SOC associated with land-use change was the second most important contributor to the baseline life cycle GHG emissions across all the scenarios considered. Considering an average, modeled changes in SOC within the geographic boundary of five counties in New York resulted in 0.29 Mg C ha −1 yr −1 (sequestration) in the first 30-cm depth soil for cropland, and −0.29 Mg C ha −1 yr −1 for grassland [4,23,27]; growing willow on previous cropland resulted in the sequestration of 96 kg CO 2eq Mg −1 or 1.10 Mg CO 2eq ha −1 yr −1 .
Among the variable parameters listed in Table 4, the life cycle GHG emissions of willow biomass was the most sensitive to root-to-shoot ratio and soil organic carbon (SOC) change (Figure 3). For a given aboveground biomass yield, increasing the root-to-shoot ratio from the baseline value (0.67) to 0.99 would reduce the net life cycle GHG emissions by 76.6 kg CO 2eq Mg −1 . Therefore, the breeding and deployment of high root-to-shoot willow cultivars that also have high yields, can play an important role in mitigating climate change while providing raw materials for bioenergy and bioproducts. Recent studies of willow growing in upstate NY reported root-to-shoot ratios of three cultivars at two different sites, ranging from 0.54 to 1.04 [41]. Considering the important role of root-to-shoot ratio on the net life cycle GHG emissions, further investigation is needed to better understand the complex relationship between root-to-shoot ratio, aboveground biomass growth, and environmental conditions (e.g., climate, soil conditions). Decreased sequestration of SOC would result in a substantial increase in the net GHG emissions of willow biomass. For instance, a reduction of the SOC from 0.29 to 0.19 Mg C ha −1 yr −1 , in the case of cropland, is equivalent to a net increase of 33.1 kg CO 2eq Mg −1 . Also, the sensitivity results depicted the opposite trends of the net life cycle GHG emissions for varying yield under both grassland and cropland scenarios. An increase in willow biomass yield resulted in a slight decrease in life cycle GHG emissions for both grassland scenarios (S1 and S2), and a slight increase for both cropland scenarios (S3 and S4). These diverging trends for the GHG emissions response to varying yield on cropland and grassland are counterintuitive. In fact, for all the scenarios, an increase in willow biomass yield would result in lower GHG emissions per Mg of willow biomass produced for the processes of field preparation, field maintenance, and harvesting because fewer hectares of land would be needed to produce the same amount of biomass. However, these fewer hectares of land mean that there would be less SOC sequestered (negative trend) for cropland scenarios, and less SOC released (positive trend) for grassland scenarios per Mg of willow biomass produced. As for root-to-shoot ratio, future research should aim to further understand the interaction between willow biomass yield and SOC, under a range of field and climatic conditions across multiple rotations. Previous LCAs of willow bioenergy systems reported results that were in the same order of magnitude as our study, but with substantial variation reflecting different selected scenarios, assumptions, and methodological choices [4,18,20]. A systematic review of data published in the literature showed GHG intensity, from cradle-to-gate, in the range of 0.6 to 10.6 g CO 2eq MJ −1 biomass for the production of short rotation willow and poplar biomass [18]. Most of the reviewed studies did not include the impact of SOC changes and sequestration in willow stool and coarse roots. Carbon sequestration in soils depends on many factors such as initial carbon level, climate, soil characteristics, and management practices [42]. However, as depicted in this study, the impact of SOC changes and carbon sequestered in roots should not be overlooked because they can shift the system from a net GHG emitter to a carbon sink. Nevertheless, when the impact of SOC was omitted, the life cycle GHG emissions for all four scenarios in this study, regardless of their prior land use (−4.0 to −3.1 g CO 2eq MJ −1 biomass ), were comparable to the range of −6.9 to −2.7 g CO 2eq MJ −1 biomass reported in the literature for short-rotation willow crop systems [20].

MSP Results
This study found that the lowest baseline MSPs occurred when willow was harvested leaf-off, while the highest MSPs occurred for the leaf-on harvests. The lowest baseline MSP of $76.41 Mg −1 was reported for the cropland + leaf-off scenario (S3), while the highest baseline MSP value of $92.97 Mg −1 was attributed to the grassland + leaf-on scenario (S2) (Figure 4). The Department of Energy's $84 dry Mg −1 threshold [43] was $7.59 Mg −1 higher than the lowest MSP modeled in this study (S3-CL + LOFF).
There were considerable differences in the baseline MSPs between the leaf-on and leaf-off values. The cropland scenarios produced baseline MSPs of $91.49 Mg −1 (leaf-on) and $76.41 Mg −1 (leaf-off). The difference of approximately $15 Mg −1 was primarily attributed to differences in harvester fuel consumption, ground speed, and throughput. Fuel consumption significantly differed between the harvesting scenarios and yielded 144.48 L ha −1 for the cropland + leaf-on scenario (S4), versus 74.89 L ha −1 for the cropland + leaf-off scenario (S3). Another key difference was the substantially lower harvest rate of 0.53 ha h −1 (S4. CL + LON) compared to the more than double harvest rate of 1.19 ha h −1 (S3. CL + LOFF). Material capacity (harvester throughput) varied considerably between the cropland scenarios with the leaf-on scenario handling 35.88 Mg h −1 , as opposed to the 74.57 Mg h −1 value for leaf-off scenario. Harvest speed was substantially higher for the cropland + leaf-off scenario at 5.93 km h −1 , while the leaf-on scenario's harvester speed was calculated at 2.67 km h −1 . Total harvest time also differed substantially between both scenarios, with the leaf-on scenario's harvest time being more than double that of the leafoff scenario. These harvest time differences resulted in a harvester rental cost difference of $26,127.99 (cropland + leaf-on), versus $12,409.92 (S3. CL + LOFF) per harvest. These trends were also notable for both grassland scenarios. Research on harvesting operations has noted that willow stems with leaves do not fall forward as easily when they are cut and are harder to feed into the harvester, both of which slow forward speed and throughput [24]. There were also slight financial differences between the grassland and cropland scenarios, albeit less substantial than between the leaf-on and leaf-off harvests. The baseline MSP calculated for the grassland + leaf-on scenario was $92.97 Mg −1 , while the baseline cropland + leaf-on scenario yielded $91.49 Mg −1 , a difference of $1.48 Mg −1 (Figure 4). The baseline leaf-off scenarios varied by $1 Mg −1 at $77.41 Mg −1 (S1. GL + LOFF) versus $76.41 Mg −1 (S3. CL + LOFF). These variations were attributed primarily to differences in land rental rates of $46 ha −1 yr −1 for grassland, versus $102.27 ha −1 yr −1 for cropland. These variations in MSP between original land uses were also due to differences in transportation distance, average annual biomass yield, and harvester and fuel consumption parameters (Table 4 and Figure 5). Figure 5 contains a breakdown of the costs and revenues associated with all four baseline scenarios. A majority of the costs came from both the harvesting process and from the total general data (land costs, internal administration costs, and stock removal). On the other hand, the total cost of crop maintenance (pre-emergent herbicide after planting, weed control (1st and 2nd season), coppice, and fertilizer) attributed the least amount to the total cost over the entire lifetime. The revenues were attributed to harvesting, both leaf-on scenarios yielded the highest revenues and the leaf-off scenarios were responsible for the lowest revenues.
In terms of the stochastic variables, the MSP was the most sensitive to changes in the willow biomass yield (Mg ha −1 yr −1 ), transportation distance to end-user (km), moisture content (% wet basis), and land rental rate ($ ha −1 yr −1 ). As yield increased, the MSP decreased across all scenarios ( Figure 6). On the other hand, as the transportation distance (km), moisture content (% wet basis), and land rental rates ($ ha −1 yr −1 ) increased, the MSPs followed the same trend. The MSPs increased considerably as land rental rates increased to at least 300 $ ha −1 yr −1 for both cropland scenarios. This highlights the importance of accounting for stochastic variables when assessing the financial viability of a feedstock.

Life Cycle Environmental and Financial Impacts
The scenario with the greatest negative GHG emissions was the scenario with the lowest MSP (S3. CL + LOFF) (Figure 7). The scenario that yielded the highest GHG emissions was responsible for the highest MSP (S2. GL + LON). Therefore, willow should be planted on less productive cropland, and a leaf-off harvest should be conducted to achieve the lowest MSP and GHG emissions. From an LCA perspective, the greatest GHG emission differences (~206 kg CO 2eq Mg −1 ) occurred between the grassland and cropland scenarios, which was driven by the increase in soil carbon when willow is grown on cropland, and a decrease in soil carbon on grassland converted to willow. The data that is available on changes in soil carbon under willow crops over multiple rotations is limited, and needs to be improved, and methods developed to verify these changes, especially since it is an important factor in the overall GHG balance of these systems. In terms of the current use of land that is suitable for willow in the region studied, the vast majority (88.7%) of it is in the cropland/pasture category. As a result, even using a mixture of cropland and grassland in proportion to its availability on the landscape would result in negative GHG emissions across the landscape. On the other hand, the greatest differences in baseline MSPs for the TEA were between the leaf-on and leaf-off scenarios, because harvesting occurs every three years over the life of the crops and makes up a large proportion of the MSP. While the lower harvester throughput in leaf-on scenarios had a large impact on MSP, the change in GHG emissions was smaller (~11.4 kg CO 2eq Mg −1 ). Results of this study are within the range of reported values in the literature, with significant variation capturing the particularity of each scenario. A study by Liu et al. [12] quantifies the environmental and financial impacts of producing bioenergy products such as biofuel, electricity, and pellet fuel from hybrid willow, switchgrass, and Miscanthus. The GHG emissions were 0.73 g CO 2eq MJ −1 (~11 kg CO 2eq Mg −1 ) for feedstock development, harvest, and delivery to a biomass processing facility, and the estimated 70 $ Mg −1 is closer to the reported value in our study for the leaf-off harvest. It appears that Liu et al. [12] did not incorporate the impact of leaves decomposition, SOC, and carbon sequestered in the roots in their LCA. A previous Ecowillow version (3.0 S) showed a median stochastic MSP of 75 $ Mg −1 [6]. Another LCA of willow biomass system in the Northeast U.S. [4] reported positive GHG emissions (28 CO 2eq Mg −1 ) for the grassland scenario, and negative GHG emissions (−127 kg CO 2eq Mg −1 ) for the cropland scenario for a 30-cm soil depth. Although the overall trends remain the same, the reported baseline GHG emissions by [4] are 15% to 37% lower (in absolute value) than the reported values in the current study, depending on the scenario. These differences are attributed to updated harvesting data that were partitioned for leaf-on and leaf-off harvests, and differences in modeling approaches. The current analysis uses stochastic modeling that inherently integrates the uncertainty associated with the variable input parameter instead of a simple average of the experimental data points. For instance, the average root-to-shoot ratio of the 10,000 data points generated based on the defined probability distribution was 12% higher than the average value from [4] for the same parameter. Using an integrated approach for the LCA and TEA, in addition to stochastic modeling, supports robust results.
This analysis showed that there were different key drivers for the baseline GHG emissions and MSPs from this system, which will have implications for future management and policy decisions related to the deployment of this cropping system. If crop management decisions are going to be made to improve the profitability of this system, the focus should be on factors like improving harvesting and yield, and using lower cost land. If the primary driver of policy decisions is on GHG balances, then factors that should be focused on are different and include the prior land use, root-to-shoot ratio, and, to a lesser extent, yield. Without a way to value the changes in carbon, these differences will create conflicts between potential policy drivers and decisions made by landowners who are growing the crop on their land.
While this study quantified both the life cycle and financial impacts of shrub willow production in the Northeastern United States, it did not include a value on carbon. Future research should focus on the inclusion of a price on carbon into the MSP of willow biomass. Specifically, the social cost of carbon (SCC) values provided by the Department of Environmental Conservation's (DEC) Establishing a Value of Carbon Guidelines for Use by State Agencies [44] can be utilized within this analysis. These values can be incorporated into the financial analysis as a revenue stream in terms of GHG emissions reduction from a baseline value, such as those from fuel oil or natural gas. The revenue stream can be used to calculate the MSP for shrub willow. Ultimately, this would provide a decrease in the MSP of willow, as the SCC revenue would provide a financial benefit. It is important to note that the SCC estimates provide an annual monetary estimate of the damage that adding a metric ton of CO 2 emissions has on the atmosphere.

CAC Values
Carbon abatement cost (CAC) values account for the decarbonization cost of switching fuel systems. The CAC is the cost of abating a unit of CO 2 and is not limited to renewable energy sources. Since combustion was not within the scope of this study, a calculation of the CAC values was not included within this study. However, CACs have been calculated in recent refereed bioenergy studies for the United States [43,45]. A study by Björnebo et al., 2018 [43] calculated the CAC values for several centralized district heating networks using combined heat and power (CHP) plants in the Northeastern U.S. that do not have access to natural gas infrastructure, and found that when biomass-fed district heating was implemented, negative CACs mainly within the range of −$250 to −$38 Mg CO 2eq −1 occurred. This was mainly due to the reduced GHG emissions and operating costs. This study also found that those district heating systems connected to the natural gas grid yielded negative GHG abatement costs from −$4500 to −$400 Mg CO 2eq −1 . When compared to the biomass CHP systems, the natural gas system resulted in a financial benefit due to the reduced operation and maintenance costs, along with combined CHP's increased electricity generation on a per unit basis. One Georgia-based study by Masum et al., 2020 [45] conducted a long-term analysis to understand the environmental and economic impacts of switching from a coal-based power plant to biomass-based feedstocks, for the purpose of electricity generation. This study found that to be competitive against a coal powerplant for the electricity sector, a $40 t CO 2eq −1 carbon tax was necessary to produce competitive bioenergy feedstocks, in order to reduce carbon emissions. In the future, research related to dedicated energy crops, such as shrub willow, can focus on the incorporation of the fuel combustion process into an integrated LCA and TEA analysis, in order to calculate the CAC values for biomass-based heating systems versus fossil fuel systems.

Conclusions
This study performed an integrated LCA and TEA of shrub willow production in the Northeastern U.S. under uncertainty. The baseline life cycle GHG emissions and MSPs were quantified and analyzed across two land use and two harvesting scenarios. The results depicted large differences between the baseline GHG emissions produced between the scenarios. In terms of baseline results, the cropland leaf-off scenario produced net negative GHG emissions at −172.50 kg CO 2eq Mg −1 and the lowest MSP of $76.41 Mg −1 . On the other hand, the grassland leaf-on scenario yielded the highest net GHG emissions and MSP results at 44.83 kg CO 2eq Mg −1 and $92.97 Mg −1 , respectively. The baseline cropland scenarios depicted negative emissions due to larger amounts of atmospheric carbon being sequestered in the roots and the soil, than generated through the willow life cycle. It is important for scientists, farmers, policymakers, and other interested stakeholders to understand both the environmental and financial impacts of biomass feedstocks biofuels, bioproducts, and biomaterials production versus other non-resources. Acknowledgments: Thank you to the National Agricultural Statistics Service Northeastern Regional Field Office for their time. Therasme is grateful for the start-up funding received from SUNY ESF which covers some of his time and the article publication fee.

Conflicts of Interest:
The authors declare no conflict of interest. The funders 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.