Spatially variable hydrologic impact and biomass production tradeoffs associated with Eucalyptus (E. grandis) cultivation for biofuel production in Entre Rios, Argentina

Climate change and energy security promote using renewable sources of energy such as biofuels. High woody biomass production achieved from short‐rotation intensive plantations is a strategy that is increasing in many parts of the world. However, broad expansion of bioenergy feedstock production may have significant environmental consequences. This study investigates the watershed‐scale hydrological impacts of Eucalyptus (E. grandis) plantations for energy production in a humid subtropical watershed in Entre Rios province, Argentina. A Soil and Water Assessment Tool (SWAT) model was calibrated and validated for streamflow, leaf area index (LAI), and biomass production cycles. The model was used to simulate various Eucalyptus plantation scenarios that followed physically based rules for land use conversion (in various extents and locations in the watershed) to study hydrological effects, biomass production, and the green water footprint of energy production. SWAT simulations indicated that the most limiting factor for plant growth was shallow soils causing seasonal water stress. This resulted in a wide range of biomass productivity throughout the watershed. An optimization algorithm was developed to find the best location for Eucalyptus development regarding highest productivity with least water impact. E. grandis plantations had higher evapotranspiration rates compared to existing terrestrial land cover classes; therefore, intensive land use conversion to E. grandis caused a decline in streamflow, with January through March being the most affected months. October was the least‐affected month hydrologically, since high rainfall rates overcame the canopy interception and higher ET rates of E. grandis in this month. Results indicate that, on average, producing 1 kg of biomass in this region uses 0.8 m3 of water, and the green water footprint of producing 1 m3 fuel is approximately 2150 m3 water, or 57 m3 water per GJ of energy, which is lower than reported values for wood‐based ethanol, sugar cane ethanol, and soybean biodiesel.


| INTRODUCTION
Using sources of renewable energy, such as biofuels, may result in cleaner, cost-competitive alternatives to fossil fuels (Sekoai et al., 2019;Winjobi et al., 2018). Cellulosic crops, crop residues, and woody biomass are promising bioenergy sources because they have shown to produce similar fuel yields per feedstock mass as first-generation biofuels such as corn-based ethanol (Lynd et al., 1991;Tilman et al., 2009). Short-rotation harvest of woody biomass is considered a major advance in bioenergy because of high rates of biomass production (Guerra et al., 2014). Eucalyptus is the most widely planted hardwood genus in the world, covering more than 19 million hectares . Eucalyptus trees are highly productive (e.g., >35 m 3 biomass/ha/year found by Albaugh et al., 2013), have a short-rotation length of 6-8 years, have used as lumber and pulp (Dougherty & Wright, 2012), and potential benefits for ecosystem restoration and carbon sequestration (Pazhavand & Sadeghi, 2020). Many parts of the world are experiencing a rapid increase in Eucalyptus plantations for biofuel (Gonzalez et al., 2011). There are several bioenergy products from Eucalyptus, including cellulosic biodiesel and ethanol (Gonzalez et al., 2011), as well as wood pellets for direct heating or electricity generation (Pirraglia et al., 2010).
Eucalyptus plantations have high water use efficiency (WUE; . Furthermore, fast-growing Eucalyptus are more efficient water users compared to slower growing trees (Otto et al., 2014). However, Eucalyptus plantations have been reported to have high water use compared to other species (Albaugh et al. 2013;Scott, 2005) and compared to the native plants that they replace (Farley et al., 2005;Ferraz et al., 2013). In fact, Eucalyptus has one of the highest ET rates among trees (Dye, 2013;Farley et al. 2005;Hubbard et al. 2010), due to morphological and physiological characteristics including high stomatal conductance, evergreen leaves, drought tolerance, and deep rooting (Whitehead & Beadle, 2004). Farley et al. (2005) observed a higher water use rate for Eucalyptus by converting grassland to Eucalyptus and pine plantations on a catchment scale. They concluded that converting to Eucalyptus would decrease the streamflow 25%, compared to conversion to pine. Maier et al. (2017) studied short-rotation Eucalyptus benthamii planting in South Carolina, USA at the plot scale and concluded that Eucalyptus had a 40% higher transpiration rate compared to pine (Pinus taeda). However, little is known about Eucalyptus cultivation impacts on specific hydrologic components, that is, baseflow versus surface runoff, and seasonality.
Proper site selection for biofuel-related land use conversion can be crucial for sustainably managing resources in a watershed (Cibin & Chaubey, 2015). An appropriately selected biofuel crop planted at a suitable location can reduce water quality impacts of biofuel development projects (Parish et al., 2012;Robertson et al., 2008). Spatial allocation of biofuel crops has been studied on different scales, from a national level in China using geographic information systems (Zhang et al., 2017) and at a watershed scale using optimization methods (Cibin & Chaubey, 2015;Femeena et al., 2018;Herman et al., 2016;Parish et al., 2012). However, the spatial variations in biomass production across the watershed are typically neglected. Biomass yield can vary significantly in cases where soil depth, soil quality, precipitation, or temperature change across the watershed.
Sustainable biofuel production requires assessing the hydrologic impacts-in terms of water use and water pollution-at the watershed scale Engel et al., 2010;Gopalakrishnan et al., 2009;Heidari, Mayer, Watkins, Propato, et al., 2019;Heidari et al., 2020;Watkins et al., 2015). Developing proper management practices to achieve high water use efficiency, while minimizing negative environmental impacts, requires quantification of Eucalyptus water demand. To fully understand the impacts of Eucalyptus development on water resources, their growth cycle and water use should be studied in more detail at sub-watershed scales.
Hydrological models have been used globally to study hydrological impacts of biofuels, especially for first-generation bioenergy crops (e.g., Lin et al., 2015;Love & Nejadhashemi, 2011;Schilling et al., 2008), but less so for second-generation bioenergy crops (Guo et al., 2018(Guo et al., , 2019Hillard, 2017). SWAT is a commonly used ecohydrological, physically based, spatially semi-distributed simulation model (Arnold et al., 2000). SWAT provides the opportunity for detailed simulations at scales ranging from tens of hectares up to watershed or river basin scales, including both hydrologic and plant growth sub-models. SWAT has been used to simulate biofuel development around the world for different crops (Babel et al., 2011;Cibin et al., 2016;Heidari, Mayer, Watkins, Propato, et al., 2019;Heidari et al., 2020;Sinnathamby et al., 2017).
The goal of this study is to determine how Eucalyptusbased biofuel feedstock cultivation will impact hydrological systems. Specific objectives are to assess the impacts of spatially varying patterns of Eucalyptus plantation (E. grandis), biomass productivity, and water use for biomass production on baseflow, surface flow, and evapotranspiration. The interannual variability of hydrologic impacts is also to be evaluated, along with the explicit tradeoff between biomass production and water use. These objectives are accomplished by calibrating and validating a SWAT model, using both hydrologic and plant growth data, for E. grandis plantations in a watershed located in Entre Rios, Argentina. Argentina is one of the largest biofuel producing countries in the world (Statista, 2019), and the Mesopotamia region of Argentina is an appealing candidate for continuing development with | 825 HEIDARI Et Al.
Eucalyptus plantations. Cultivation of E. grandis, which is considered to be one the most important Eucalyptus species globally (Dougherty & Wright, 2012), is expanding rapidly in the region (Phifer et al., 2017).
While SWAT has been used to study hydrologic processes in various watersheds in Argentina (Cisneros et al., 2011;Havrylenko et al., 2016;Romagnoli et al., 2017;Schwank et al., 2014;Troin et al., 2012), to the authors' knowledge, this is the first application of SWAT that focuses on improving E. grandis growth parametrization and investigating the hydrologic impacts of Eucalyptus plantations for biofuel development. Considering the rapid expansion of the E. grandis in this region of Argentina, there is a need for more studies of the water use, management, and productivity of the plantations. In this work, the SWAT hydrologic and biomass growth models are calibrated and used to assess the impacts of spatially varying patterns of E. grandis plantation, biomass productivity, and water use for biomass production. The SWAT model is used to determine the feedstock stage water demand for biomass, fuel, and energy production, as well as impacts on specific hydrologic components (baseflow, surface flow, and evapotranspiration) and the interannual variability of those impacts.

| Model setup, calibration, and validation
The selected watershed (see Figure 1) is representative of the Argentinian Mesopotamia region. Land cover in the region typically consists of rangelands, crops such as soybeans (Modernel et al. 2016), natural forests (Espinal), orange orchards, and rivers and wetlands draining into the Uruguay River to the east. The Yuqueri Grande-Concordia hydrologic station (Secretaría de Infraestructura y Política Hídrica [SIPH], 2015) was selected as the watershed outlet, and daily flow data for the period 1991-2013 were used for calibration and simulations. The contributing area to the gage was found to be 625 km 2 , using the 30-m resolution digital elevation model from USGS (U.S. Geological Survey [USGS], 2015) and ArcGIS 10.3. A customized streamline shapefile from the Argentina National Institute of Geography (Instituto Geografico Nacional, 2015) was used to improve the streamline delineation process.
Land use-land cover (LULC) maps from 2002 to 2003, 2005 to 2006, and 2013 to 2014 and soil maps were obtained from INTA. Land use-land cover classifications were made with high-resolution images including LANDSAT 5 and 8 with a spatial resolution of 30 m from USGS (U.S. Geological Survey [USGS], 2016). For each growing season, a majority voting approach was applied considering five supervised classification methods: Maximum Likelihood, Support Vector Machines, Random Forest, LOGIT regression, and Neural Networks (Waldner et al, 2016). Classes included orange orchards, agriculture, native forests, Eucalyptus plantations, and rangelands. Ground truth data for training and validation were compiled from different sources, including georeferenced photos, visual identification (in situ observation), georeferenced voice recordings, land owner's information, and visual interpretation of Very High-Resolution (VHR) images. The overall accuracy for each LULC maps were 0.89 for 2002-2003, 0.91 in 2005-2006, and 0.95 in 2013-2014. The series of LULC maps indicated a significant decline in natural forest land (−60%) and orchards (−76%), while the area planted with Eucalyptus expanded by slightly more than 100% over the 12-year period (see Table S1 for a summary of land cover change analysis). Preliminary assumptions for determining areas for biofuel development were that plantations would not compete with food crops (Paine et al., 1996) and no wetland areas would be converted. However, the land F I G U R E 1 (a) Land use/land cover map of the watershed, locations of the sub-basins, and precipitation gauges with average annual precipitation (in mm); (b) Land cover distribution in the watershed; and (c) Study site location within the state of Entre Rios, Argentina cover analysis indicated a large decrease in orange orchards and a slight variation in agriculture and rangelands over the time period evaluated.
Maximum and minimum daily temperature and daily precipitation, relative humidity, and wind speed data were compiled by INTA for the Aero Concordia weather station (see Figure 1). In addition, Climate Forecast System Reanalysis (CFSR) global weather data from four more nearby stations were added to the SWAT weather database. Spatial interpolation of the climate data indicates that the area receives an average of 1220 mm of precipitation annually, with the majority of rainfall occurring during October, November, and April (see Figure 2). Precipitation is higher in the eastern portion of the watershed (see Figure 1). The average annual temperature in the watershed is 19.4°C with slight variation across the watershed. Figure 2 shows the intra-annual variation of streamflow, temperature, precipitation, and SWAT-simulated estimates of evapotranspiration (ET) using the Penman-Monteith method. The highest monthly average streamflow is in October, as a result of heavy rain events and average ET. Runoff efficiency, the ratio of annual stream flow to annual precipitation in the watershed, was 0.22 over the study period. Even though the months of January and February receive around 100 mm precipitation on average, they are among the lowest streamflow months due to higher temperature and ET. The interannual variability in precipitation and temperature is shown in Figure S1, and seasonal climate variability is illustrated in Figure S2. In the simulation period of 1993-2013, annual precipitation ranged from a low of 833 mm in 1999 to a high of 1963 mm in 2002. The average annual temperature ranged from 18.4°C in 1998 to 20.9°C in 1997.
ArcSWAT version 2012.10_4.19 (Winchell et al., 2013) was used for setting up the model. The watershed was divided into eight sub-watersheds to assess the potential spatial variability of hydrologic impacts associated with Eucalyptus cultivation. SWAT further divides the sub-basins into noncontiguous hydrologic response units (HRUs), which represent homogeneous areas within each sub-basin with unique combinations of land use, soil type, and slope class. During the HRU definition, thresholds of 0%, 5%, and 15% were selected for LULC, soils, and slope classes, respectively, resulting in 185 HRUs. Rivers and wetlands comprise 18% of the watershed land cover, and thus wetlands were considered in the model. Wetland functionality is described in detail in SWAT theoretical documentation (Neitsch et al., 2011) and  ;Heidari, Mayer, Watkins, Propato, et al. (2019). The variable storage routing method was selected for this study.
Calibration and validation focused on both the hydrological and plant growth components of the model. The analysis was performed for 21 years (from 1993 to 2013) to include a combination of dry and wet years in both the calibration period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) and validation period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013). Periods with missing or unreliable data, attributed to a bridge construction project that impacted the stream gage measurements in some periods, were omitted from the goodness-of-fit calculations. The parameters controlling LAI were adjusted during the hydrologic calibration to affect the simulated ET. The calibration process included calculating the heat units for E. grandis in the region and changing the shape coefficients of the LAI growth curve. Parameters related to LAI development stages along with potential heat units were adjusted to reflect the evergreen nature of the tree, similar to Alemayehu et al. (2017). The maximum LAI, defined as the LAI that can be reached in the absence of water and nutrient stress, was adjusted based on literature values Maurice et al., 2010) and field measurements (J. Licata, personal communication, 30 Nov 2017).
The hydrologic calibration method was similar to Heidari, Mayer, and Watkins (2019); Heidari, Mayer, Watkins, Propato, et al. (2019), which included separating baseflow and surface flow (Arnold & Allen, 1999). This analysis resulted in the ratio of baseflow to total flow ranging from 0.33 to 0.51 on an annual basis. The next step was to conduct a sensitivity analysis using the p value and t-statistic sensitivity tests in SWATCUP SUFI2 (Abbaspour, 2013). Finally, the sensitive parameters were adjusted in groups. Final adjusted parameter values are presented in Table S2. Performance metrics for calibration included the Nash-Sutcliffe Efficiency (NSE; Nash & Sutcliffe, 1970), coefficient of determination (R 2 ), and percent bias (Pbias; Gupta et al., 1999). The ratio | 827 HEIDARI Et Al. of baseflow to total flow was also required to be within the historical range.
Biomass growth is dependent on LAI and solar radiation and does not directly influence the hydrologic cycle. Therefore, growth rates and the maximum biomass were calibrated after the LAI and hydrologic calibration. Using the final LAI parameters from the hydrologic calibration, the biomass growth was further calibrated by adjusting the radiation use efficiency and light extinction coefficient. Reported values for these parameters from De Costa and Jayaweera (1996) and  informed the biomass growth calibration process. It was assumed that E. grandis trees are planted as saplings, and the full growth cycle was simulated. In the simulations, LAI increases year by year until it reaches the specified maximum LAI, and biomass also increases each year until the trees are harvested ( Figure S3). The biomass growth calibration accounted for losses during the dormancy period, and simulated biomass at the time of harvest matched reported biomass yield in the area (INTA, 2016). Table S3 lists the adjusted plant parameters. Full descriptions of each parameter are presented in the SWAT theoretical documentation (Neitsch et al., 2011).

| Modeling scenarios
Biofuel development scenarios were formulated considering a number of variables, including the land cover types being replaced, locations of feedstock cultivation (e.g., in headwaters or downstream sub-basins), spatially variable soil fertility, and whether or not irrigation is applied. SWAT model simulations were performed for the period 1991-2013, using the corresponding hydroclimatic time series. This period included a 2-year warm-up period (1991)(1992)(1993) to establish initial conditions, followed by a 21-year period  for scenario evaluation. This period allows for a range of climate conditions to be represented, as well as several harvesting rotations. Specifically, the E. grandis trees were planted at the end of August and were harvested at the end of May, with two 7-year rotations and a 6-year rotation represented in the 21-year simulation (i.e., initial planting is towards the end of the first year of the SWAT simulations). The 6-to 7-year rotation length was selected based on previous studies Dougherty & Wright, 2012;INTA, 2016). Simulated plantations were fertilized (100 kg N/ha/year) to prevent nutrient stress. The scenarios for various land areas converted to Eucalyptus consider watershed, sub-basin, and HRU scales, as described in Table 1.
A bi-criteria optimization model was developed to determine Pareto-optimal combinations of sub-basins, that is, where B is cumulative biomass production over the simulation period, Q is total streamflow, s is the sub-basin index, and S is the total set of sub-basins. The optimization procedure was based on results from the one sub-basin-at-a-time scenarios, formulated as a knapsack problem to maximize total biomass production subject to a single constraint on the allowable change in total streamflow. The tradeoff curve was generated by starting with a low level of allowable change in streamflow and then incrementally relaxing the constraint to allow more conversion.

| Model evaluation
Comparison of simulated and observed monthly discharges, shown in Figure 3, demonstrates good performance of the hydrologic model. The Nash-Sutcliffe efficiency (NSE) of 0.6, R 2 of 0.55, and Pbias of less than 10% for the entire simulation period indicate satisfactory hydrologic model performance (Moriasi et al., 2007). Table S4 shows the goodness-of-fit statistics for the calibration and verification periods. The overall water balance matches measurements of the dominant hydrologic processes, including an average of 812 mm/year ET and 300 mm/year streamflow from an average precipitation of 1220 mm/year. The ratio of surface flow to baseflow is also maintained within the historical range (0.33-0.51) on an annual basis. The average annual deep percolation of 100 mm/year is less than 9% of the precipitation and is physically justified by the eastern portion of the watershed having deep soils that drain to the Uruguay River. Model calibration also resulted in LAI values matching the regional measured values and values reported in the literature, which ranged from 4.0 to 4.2 Maurice et al., 2010). Calibration of the biomass growth model in SWAT resulted in the most productive HRUs matching the highest reported yields for the area, approximately 100 tons/ha/rotation. The average simulated biomass yield was 75 tons/ha/rotation, also matching the average reported values for the region (INTA, 2016; see Figure S3 and Table  S5 for annual biomass production and LAI development simulated with SWAT). The range in biomass production is attributed to variability in water availability (Smethurst et al., 2003). The simulated N uptake rate of 65 kg/ha/year is within the medium-high range reported by . Successful LAI and biomass calibration was an important goal of this study attempting to address reported limitations of previous research on the hydrologic impacts of Eucalyptus plantations. These limitations have included utilizing plant parameters from trees other than the species in question (Brown et al., 2015), generalizations or application of empirical curves of water use based on vegetation type (Almeida et al., 2016), and simulating only portions of the life cycle (Brown et al., 2015). Hence, the plant growth parameters calibrated in this study (Table S3) can be a valuable reference for applications in other regions.
As shown in Figure S3 and Table S5, simulation results indicate that yields were generally resilient to droughts occurring in the hydrologic record, although some sensitivity to drought was noticed near the end of the rotation, when LAI is near a maximum (e.g., 1999 and 2006), resulting in a slightly lower yields for those rotations than the final rotation. This result is consistent with the fact that most of the E. grandis plantations are located in the eastern part of the watershed, where soils are deep (greater than 1250 mm; see Figure 7). The low rainfall in 1995 and 2008 also appears to have affected LAI growth, with lower values of incremental LAI in those years than in 2002, the wettest year in the record. Each of these years was the second year of a rotation. The effect was less pronounced in 2008, when the monthly rainfall distribution was more uniform across the growing season. Older Eucalyptus trees with well-established roots are able to access deeper water sources during the growing season while younger trees would be expected to be more sensitive to drought (Brando, 2018;Engel et al., 2005).

| Analysis of watershed-scale impacts
Simulation results from all scenarios are summarized in Table  2. The intensive scenario had an average yield of 77.1 ton/ ha/rotation (cumulative biomass = 9 × 10 6 ton). Under the intensive production scenario, streamflow was reduced at the watershed outlet on average by 28%. The surface flow was reduced by an average of 24%, with the greatest relative change in December through March (34% average decline). The average overall reduction in baseflow was 31%, with the months of January to April being the most impacted months, with an average baseflow decline of 39%. Figure 4 shows the changes in monthly average total, baseflow, and surface flow for the intensive scenario and the base case.
The E. grandis plantations had the highest annual average ET rate (842 mm/year) among the terrestrial LULC classes in the basin, which was 24% higher than the average of 638 mm/ year for all terrestrial LULC classes (the average over the watershed, including water bodies, was 812 mm/year). This E. grandis ET rate is similar to 880 mm/year as reported by  for high-productivity Eucalyptus for a location in Brazil with a climate similar to this study area. In the intensive scenario, the conversion increased the average annual ET over the newly converted land (319 km 2 ) by 32% (204 mm), corresponding to a 14% increase over the watershed (625 km 2 ). The large increase in ET in the converted area was due to large areas of rangelands being replaced. Conversely, converting orange orchards did not result in a large ET difference per unit area, as orange trees have similarly high ET rates and canopy interception. Figure 5 shows the monthly average ET for the intensive T A B L E 2 Summary of biomass production and hydrologic impacts for each scenario. All changes are relative to the Base case

F I G U R E 4
Average monthly baseflow, surface flow, and total flow for Base case and Intensive scenario scenario versus the base case during the simulation period. The substantially higher ET rate in January to April and September to December correlated to the greatest reductions in monthly streamflow shown in Figure 4. Interannual variation in climate produced some severe decreases in streamflow due to conversion to E. grandis. During the driest years (1999, 1995, and 2008), the average precipitation was 855 mm precipitation (compared to the mean annual precipitation of 1223 mm), and there was a 53% decline in the annual streamflow under the intensive production scenario. In wet years (25th percentile high annual precipitation), streamflow decreased only by 20%, on average. An annual precipitation of about 1200 mm was usually sufficient to saturate soils and fill the wetlands to capacity. The exceptions were 1997 (1182 mm) and 2009 (1332 mm), which had an average 38% decline in annual streamflow. These years both followed dry years, which caused large declines in soil moisture and wetland volume. Figure 6 shows the cumulative distribution (exceedance probability) of monthly streamflow values for the base case and intensive scenario over the entire simulation period. A significant shift downward in streamflow is observed as a result of replacing existing land cover with E. grandis plantations, especially for the lower flows. The shift was smaller for higher flows as they are associated with heavier rainfall.
The extreme scenario had an average biomass yield similar to the intensive scenario, with 77.2 ton/ha/rotation, but produced a higher cumulative biomass (12 × 10 6 ton) as a result of converting 83% of the watershed to E. grandis. This conversion increased the average annual ET by 18%, causing a 37% decline in the average annual streamflow. When 434 mm/year of irrigation of E. grandis was added to the intensive scenario, the number of water stress days decreased by 85% and the cumulative biomass production of the watershed increased to 12.3 × 10 6 ton, an increase in 36% over the non-irrigated intensive scenario.

| Variability in biomass productivity due to spatial variations in soil properties and climate
The simulation results indicate a wide range in biomass production at the HRU scale (average area = 3 km 2 ), from 37 to 97 ton/ha/rotation for the intensive scenario. Figure 7 shows the variation in soil depth, precipitation, and yield across the sub-basins. The most critical spatially variable parameters for determining biomass productivity were precipitation and soil depth. The lowest productivities were associated with shallow soils (less than 50 cm deep), which reduce growth because the reservoir of available soil water is small, leading to water stress during low-rainfall or high-ET periods. Comparing HRUs with similar soil depths across sub-basins, different yields were mainly due to precipitation differences in these sub-basins. In wetter sub-basins, the relatively high precipitation maintained the water content of the soil, leading to a reduction in water stress. The results in Figure 7 allow comparison between the lowest biomass yield (subbasin B) and highest biomass yield (sub-basin H) sub-basins. Sub-basin H is typical of the downstream portions of the watershed, comprising the highest soil depths and precipitation. Sub-basin B is typical of the upstream portions, where soil depths are shallower and annual precipitation is about 200 mm less than in the downstream sub-basins. In the upstream sub-basins, the trees fail to reach the maximum LAI, which is reflected in their lower water use and biomass productivity. In the intensive irrigation case, the additional water increased biomass yield by 50% in the upstream sub-basins.
Using results from the intensive scenario, the HRUs were sorted from the highest biomass productivity to the lowest, and high-yield HRUs were defined as those having a productivity of more than 75 ton/ha/rotation (in the upper half of the reported range of 50-100 ton/ha/rotation). The highyield HRUs were then grouped so as to cover approximately one-third, two-thirds, and the total area of high-yield HRUs (a total area of 213 km 2 ). The simulation results in Table 2 show that converting two-thirds of the highest yield HRUs (HY2) resulted in the highest productivity (83.7 ton/ha/rotation) among all the non-irrigated scenarios simulated in this study. Table 2 also shows, however, the high water cost per biomass for the HY scenarios. This is due to most of the high-yield HRUs being located in sub-basins H and G, where incrementally higher yields were achieved but required significantly more water use, which reduced the overall water use efficiency. Figure 8 summarizes the tradeoffs between biomass production and streamflow impacts at the main outlet for the scenarios involving conversion of each sub-basin, one at a time. In the base case, the cumulative biomass yield was 1.6 × 10 6 ton, or an average yield of 75 ton/ha/rotation. Subbasins F, G, and H were inferior to the other sub-basins (i.e., they plotted below the Pareto optimal curve) because E. grandis planted in these sub-basins had relatively high water use and a greater impact on streamflow at the outlet compared to the other sub-basins. The relatively high impact is attributed to higher precipitation rates and deeper soils in these subbasins, which allowed for higher evapotranspiration rates. In contrast, planting in sub-basins A and B produced a considerable amount of biomass with a relatively small decrease in the streamflow. This was surprising as these sub-basins had high local impacts at a sub-basin level (Table 2). However, sub-basins A and B had small watershed-wide impacts because they received lower precipitation amounts and their contribution to the total streamflow at the outlet was relatively small (see Table 2). Further investigation of sub-basins A and B helped to understand how the hydrological impacts were dependent on which land cover was replaced, as well as the presence of water bodies. Sub-basin A experienced a high local impact (at the sub-basin level) since it was dominated by rangelands and it had a small area covered by rivers and wetlands. In sub-basin B, even though the total amount of converted area was greater, local hydrologic impacts were mitigated due to a larger area in this sub-basin being covered by water (Figure 1). Figure 9 shows the biomass production-water impact tradeoff analysis generated by the optimization model (see detailed results in Table S6). At low levels of allowable change in streamflow (less than 3%), only one sub-basin was converted at a time (i.e., E, F, C, A). As the streamflow constraint was relaxed, the model continued with the best combination of two or more sub-basins until all sub-basins were converted. During the optimization, sub-basins A and C were picked the most frequently (22 and 19 times, respectively, out of 29 solutions), even though sub-basins H and G had the highest productivity (selected 5 and 6 times each). This result is explained by the fact that the high biomass yield in those two basins came with the cost of high water consumption. In other words, the biomass production per unit of water consumption was highest for sub-basins A and C. A notable water-efficient solution, corresponding to ABC, produced 3.0 × 10 6 additional tons of biomass with only a 9.9% decrease in total streamflow relative to the base case. For comparison, the intensive scenario produced an additional 7.0 × 10 6 tons of biomass but resulted in a 28.8% decrease in total streamflow. Point ABCDE was also a critical point, as the slope of the tradeoff curve steepened beyond this point due to the optimization model being forced to select sub-basins F, G and H in the rest of the solutions. Sub-basins G and H were the ones least selected because of their low productivity per unit of water consumption.

| Green water footprint
Water footprints represent the total water consumption associated with a production system, with green water defined as precipitation that is stored in the soil and available for evapotranspiration, and blue water defined as water extracted from rivers, lakes, and aquifers (Chiu & Wu, 2013;Hoekstra & Chapagain, 2006). Recent research has highlighted the importance of focusing on green water resources in the assessment of water scarcity and potential ecosystem service impacts (Schyns et al., 2019). Table 3 summarizes the range of water requirements for different biofuel production T A B L E 3 Water requirement for biomass, fuel, and energy production estimated from this study and others scenarios estimated in this study and several others. The calculations for Table 3 were based on an assumption of using the total aboveground biomass (stems, branches, and leaves) with no losses (Guerra et al., 2016). Furthermore, this study only reports the water use at the farm gate, considering that total water use in the life cycle of biofuels is dominated by the feedstock production stage (Gerbens-Leenes et al., 2009).
Water use values in Table 3 mostly account for rainfall and soil moisture and can be considered green water, except for a few cases that include irrigation, which is categorized as blue water. On average, simulations conducted herein indicate a green water requirement of 0.79 m 3 to produce 1 kg of dry biomass, or 1.26 kg of dry biomass would be produced from 1 m 3 of water. This value is similar to that found by Maier et al. (2017), 0.69 m 3 water/kg dry biomass for shortrotation Eucalyptus benthamii in South Carolina, USA. Moreover,  reported a similar but slightly lower range (0.31-0.62 m 3 /kg wet biomass) for Eucalyptus production in Brazil.
Assuming each kilogram of biomass can produce 0.32 kg of fuel (GREET, 2016), and neglecting the water used at the refinery, an average of 2150 m 3 of water would be used to produce 1 m 3 of biodiesel. Further assuming this liquid fuel would have an energy content that is similar to conventional diesel fuels, 43 MJ/kg (Lemmon et al., 2019), results in a water footprint of 57.1 m 3 water/GJ, or 206,000 l/MWh. The "additional water" shown in Table 3 is defined as the difference in water use (ET) between the intensive case and the base case. This represents the direct hydrologic impact of converting land to E. grandis, which is just 0.17 m 3 /kg of biomass, or 14.6 m 3 water/GJ of energy. Notably, these water use estimates for E. grandis are orders of magnitude lower than what Dominguez-Faus et al. (2009) reported for irrigated corn and soybean. Moreover, the 57.1 m 3 water/ GJ found in this study is significantly lower than the reported values for wood-based ethanol (Schyns et al., 2017) and ethanol from sugarcane and soybean (Rodriguez et al. 2018). This indicates that planting E. grandis in the case study basin can be a water-efficient method for biofuel feedstock production, especially if plantations are located on deep fertile soils which, considering the region's high average annual precipitation, will eliminate the need for irrigation.

| DISCUSSION
Leaf area index is a key parameter for plant growth models. It is related to photosynthesis, water and nutrient use, rate of growth, and accumulation of dry matter (Ishak & Awal, 2007;Smethurst et al., 2003). Similarly for SWAT, LAI is a key parameter for simulating actual ET and biomass production (Neitsch et al., 2011). LAI measurements (Maire et al., Smethurst et al., 2003) indicate the nonlinear nature of LAI development over time in Eucalyptus trees. However, in SWAT, LAI increases at a constant rate until it reaches a specified maximum allowable LAI, and the annual growth rate is limited by a single parameter (number of years to maturity). This simplified model of LAI development can lead to inaccurate estimates of water use and annual incremental biomass production. LAI development in SWAT can be calibrated by changing the number of years to maturity, but this parameter also impacts the ratio of aboveground biomass to total biomass. Thus, there may be a tradeoff between accurate modeling of aboveground biomass and total biomass production.
The 6-to 7-year rotation length used in this study was selected based on previous studies Dougherty & Wright, 2012;INTA, 2016), with the assumption that trees are harvested before the LAI enters a plateau or descending phase (Smethurst et al., 2003). Forrester et al. (2012) and (2013) indicate that the LAI reduction can occur after the first 3 years of fast-growing E. grandis due to pruning or thinning, but this management practice was not considered due to a lack of measured values. Furthermore, SWAT does not simulate a descending phase of LAI due to age.
Another limitation in the growth model is the dormancy period. In SWAT, trees, perennials, and cool-season annuals go dormant as the day length nears the minimum for the year. Furthermore, the LAI starts to decrease to a specified minimum leaf area index during the dormant period. Both of these model assumptions are inaccurate for an evergreen tree such as Eucalyptus. Despite the improvements to modeling LAI in this study (see LAI-related parameters in Table S3), a dormant period was simulated for 2 weeks in winter (mid-July), when the minimum day length occurs. However, the biomass growth calibration accounted for losses during the dormancy period, and simulated biomass at the time of harvest matched reported values for the region.
This study reports the water footprint of biodiesel at the farm-gate level. Farm-gate level water use is commonly reported in biofuels research (e.g., Fazio & Monti, 2011). However, this study evaluated the gross production of bioenergy rather than the net production, meaning that the analysis did not account for energy inputs in the production chain. Neglecting energy inputs means that the water footprint will be underestimated, especially when bioenergy production systems have large energy inputs.
Potential water quality impacts of E. grandis plantation development were not considered in this study due to a lack of measured data for model calibration and validation. This is an important area for future research, as additional tradeoffs between water quality and quantity are likely to be relevant to biomass production and the maintenance of ecosystem services (Zhao et al., 2020).

| CONCLUSIONS
The main objectives of this work were to predict the hydrological impacts and evaluate the water-biomass tradeoffs associated with the development of E. grandis plantations for bioenergy production. Hydrologic model results indicated that the ET rates of E. grandis were the highest among the local terrestrial LULC classes in the watershed, which resulted in a decline in the streamflow the amount of which depended on the area and location of the plantations. For an intensive scenario of converting rangelands, orange orchards, and forest (62.5% of the watershed), an average annual decline of 28% in the total streamflow (including both surface and baseflow) was simulated. The greatest decline occurred during months of February, January, and March, with an average decrease of 37%.
Planting E. grandis in different parts of the watershed resulted in a wide range of biomass productivity (37-100 tons/ ha/rotation), due to variability in soil depth and precipitation across the watershed. Water-biomass tradeoffs resulted from the more productive plantations having higher ET rates and consequently greater impacts on water yield at the watershed outlet. To some extent, the tradeoffs could be mitigated by accounting for the land cover being replaced and the amount of water bodies in the area. The ET rate was higher for open water than all terrestrial LULC classes, making it a controlling hydrologic process for the sub-basins.
Based on model results, the average green water footprint of biodiesel produced from E. grandis was 0.79 m 3 water per kg of dry biomass, or 57.1 m 3 per GJ energy assuming conversion to biodiesel. The direct hydrologic impact of converting land to E. grandis was only 0.17 m 3 /kg of biomass, or 14.6 m 3 water/GJ of energy. These water intensity (m 3 water/ GJ) estimates for E. grandis are orders of magnitude lower than reported values for ethanol from irrigated corn and soybean; and they are significantly lower than reported values for ethanol from sugarcane and soybean, as well as other reported values for wood-based ethanol.