Multi-objective biofuel feedstock optimization considering different land-cover scenarios and watershed impacts

1.   Introduction

The Energy Independence and Security Act of 2007 (EISA) has pushed for the expansion of biofuel target volumes and extended the ramp-up through 2022 ^[1]. The new Renewable Fuel Standard (RFS), expanded and extended in the EISA, required the use and production of 9 billion gallons of biofuels for 2008 and a target of 36 billion gallons in 2022. Of these 36 billion gallons, at least 16 billion should be developed from cellulosic biofuels, and no more than 15 billion gallons derived from corn ethanol. Additionally, it is becoming more of a policy priority to identify sustainable approaches to bio-energy production. The EISA also requests that federal agencies begin to identify, and report environmental concerns linked to biofuel production, for example, water and soil quality and land productivity.

In order to address the sustainability of feedstock production, assessment of environmental impacts is required. The Soil and Water Assessment Tool (SWAT) is perhaps the most comprehensive environmental modeling software available and most effective in predicting sediment, nutrient, and chemical yields to streams and rivers resulting from agricultural management practices in complex watersheds ^[2]. However, the SWAT model does not have an optimization method that will allow users to identify the optimal spatial placement of land-cover and management practices. Considering the mandated biofuel production levels in EISA, it becomes a challenge to select suitable feedstocks and locations for cultivation to achieve these production targets while preserving natural resources. Biofuel crop selection will not be uniform and will be based on regional factors, productivity, and sustainability, as crop yields respond differently to climate, soil, and management. To meet production demands and sustainability criteria, it is necessary to develop innovative strategies and optimization methods to assess production amounts and impacts. Biofuel production sustainability can be achieved by scientifically assessing potential effects of biofuel feedstock cultivation on water quality and quantity, sediment, and nutrient losses in runoff.

In recent decades, corn starch from kernels has been the primary input in ethanol production ^[3], which are grown with the highest fertilizer and pesticide use than any major crop in the U.S. ^[4]. Manure or fertilizers used on croplands may increase Phosphorus (P) and Nitrogen (N) concentrations in streams and lakes, increasing eutrophication of surface waters. Runoff and erosion are the primary means of P movement in surface runoff ^[5]. The extent of a hypoxic zone in the Gulf of Mexico, caused by excessive nitrogen concentration, is a result of these practices ^[6]. Perennial warm-season grasses (WSG's) are particularly effective in removing nutrients from runoff. For instance, a switchgrass barrier combined with a fescue filter strip can reduce organic N loss by 57~70%, and particulate P loss between 50% and 68% ^[7]. Loss of organic N and particulate P in runoff is correlated with sediment loss. As the barrier of switchgrass width increases, the nutrient loss in runoff decreases exponentially. Besides removing nutrients from runoff, WSG's lessens sediment transport in runoff. Sediment deposition is promoted with the slow movement and temporary ponding of runoff within grasses. McGregor et al. ^[8] reported that no-till soils under cotton lose 5.2 Mg/ha of sediment without switchgrass buffers, whereas only 2.2 Mg/ha is lost with switchgrass buffers. Furthermore, 91% of sediment was trapped with narrow switchgrass hedges ^[9,10].

WSG's have become a promising biofuel crop with many environmental benefits, including the displacement of fossil fuels, reduction in net CO₂ emissions through soil organic carbon (SOC) sequestration, and improvements in soil and water quality ^[11]. However, bioenergy crops' production may alter the hydrology and ecosystem services of a particular region, and the impacts may not always be the same. Impacts may be negative or neutral, depending on crop selection and management practices, while others may offer improvements in water quality and other positive benefits (e.g., afforestation and reforestation) ^[12]. Different agricultural management practices, such as the heavy use of nitrogen fertilizers, increased tillage, and the selection of crops for the given soil conditions and climate, have effects on watersheds and the environment. For instance, switchgrass and Miscanthus have the potential to reduce erosion and nutrient losses within the watershed ^[11]. On agricultural lands, erosion processes are the primary means of the movement of pollutants ^[11]. Excessive amounts of non-point source pollutants can contribute to eutrophication of the receiving water bodies and impair water quality.

2.   Background

Recent articles have addressed the importance of developing decision support tools that include ecosystem services to reduce the environmental impacts associated with agricultural systems ^{[13,14,15,16,17]}. However, there is a lack of research or information that quantifies water and soil quality, and other ecosystem services, despite the expected increases in second-generation biofuel feedstock production. Related research has not identified the soil and water impacts associated with the conversion and cultivation of biofuel feedstock crops or investigated the approach to managing agricultural landscapes. Additionally, very few have used SWAT with evolutionary optimization algorithms to improve agricultural management in watersheds. Those that have integrated SWAT with genetic algorithms have developed decision-making models mainly for Best Management Practices (BMP's) and detention basin locations to reduce pollutant loads and pesticide control. For example, Kaini et al. ^[18,19,20] and Artita et al. ^[21], developed a variety of methods for evaluating cost-effectiveness, optimum combination of detention ponds, parallel terraces, filter strips, and other BMP's to reduce pollutant loads. They proposed a model that can identify least-cost combinations based on size, type, and BMP location ^{[18,19,20,21]}. Maringanti et al. ^[22,23] developed a multi-objective tool for the selection and placement of BMP's for pesticide control by combining a genetic algorithm with a distributed parameter watershed model. According to Maringanti, other optimization models that had a dynamic linkage with water quality models could only analyze small-scale watersheds due to increased computational time. Instead of having a direct interface with the watershed model, Maringanti developed a database of pollution and cost information of the BMPs under analysis, allowing them to perform the optimization at a much larger scale.

Tallis and Polasky ^[19] presented the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) as part of a computer-based model designed for decision making for biofuel feedstock cultivation, considering a comprehensive set of ecosystem services. Herman et al. ^[24] developed a strategy based on a genetic algorithm to maximize steam health by coupling several hydrological models, including SWAT. However, the optimization of biofuel crops and soil effect was not considered. Gitau et al. ^[25,26] utilized a BMP database to optimize BMP placement and cost with a genetic algorithm and SWAT. Additionally, Muleta and Nicklow ^[27] developed an integrative modeling approach to simultaneously limit sediment yield and maximize farm-level profit by coupling SWAT with a Strength Pareto Evolutionary Algorithm (SPEA). However, the integrated model is simulated only at the approximate spatial scale of a farm, without addressing watershed management. Ng et al. ^[28] used SWAT to model four different land-cover change scenarios and their potential effects on riparian nitrate loads from cultivating Miscanthus instead of conventional crops in the Salt Creek watershed in East-Central Illinois. Nevertheless, the simultaneous optimization of multiple objectives was not considered.

Physically testing all possible land-cover change combinations to determine their optimal placement is not cost-effective nor feasible. Therefore, the main purpose of this study is to present a multi-objective optimization methodology that can be applied in properly calibrated and validated SWAT projects to assist researchers and stakeholders in understanding the environmental impacts of agricultural management practices, particularly when cultivating biofuel feedstocks to meet legislative benchmarks. The proposed MOEA uses outputs from the SWAT model to find the optimal spatial placement of land-cover changes to identify a possible landscape scenario that potentially would minimize environmental effects while maximizing biomass yields. The functionality of the proposed model is demonstrated with the analysis of three perennial grasses, namely switchgrass (Panicum virgatum L.) eastern gamagrass (Trip-sucum dactyloides L.) and johnsongrass (Sorghum halepense (L.), Pers.), each with three nitrogen application levels [33 29]. These grasses and nitrogen application levels were selected to create multiple solutions to execute the MOEA.

3.   Materials and methods

Multi-objective optimization problems require the simultaneous optimization of two or more objective functions that may conflict with each other. These problems usually find a set of Pareto Optimal solutions. The main idea behind the Pareto dominance criterion is to compare all solutions against each other, where the best-fitted solutions dominate the weaker ones, which is said to be dominated. The boundary defined by the set of non-dominated solutions is known as the Pareto-front that is presented as the solution set to the multiple objective optimization problems. The solutions that are part of the Pareto set belong to a category in which there is no mathematical foundation to discard any of those solutions over one another. Since they cannot be eliminated again, they are presented as a set of "options" that are all "equally" good at optimizing the desired objectives. The Pareto dominance criterion can be formulated as follows:

A solution s dominates solution s' if f_i (s) ≤ f_i (s′) ∀i, i ∈ {1, 2, … N} and there exists some i ∈ N for which f_i(s) < f_i(s′). In other words, solution s dominates another solution s' when solution s is no worse than solution s' in all objectives and solution s is strictly better than solution s' in at least one objective ^[30]. The next sections introduce the multi-objective optimization model developed for the existing Soil and Water Assessment Tool (SWAT).

3.1.   Soil and Water Assessment Tool (SWAT)

The multi-objective optimization model developed in the present work, uses SWAT version 2016/Rev 664, to quantify the environmental impacts associated with land-cover change at the watershed scale. SWAT was developed by the United States Department of Agriculture (USDA) and the Agricultural Research Services (ARS) with almost 30 years of effort and has gained international acceptance in hydrological and pollutant load assessments ^[31].

SWAT is a physically-based model that effectively predicts sediment, nutrient, and chemical yields to streams and rivers resulting from agricultural management practices in complex watersheds. The model includes major components such as weather, hydrology, soil properties and temperature, land management, plant growth nutrient parameters, pesticides, bacteria, and pathogens. SWAT divides watersheds into multiple sub-basins, which are further subdivided into Hydrologic Response Units (HRUs). HRUs are subdivisions with similar land cover, soil characteristics, and slope class scattered throughout the sub-basin. This component makes the SWAT execution faster than those with single fields by evaluating each HRU separately and then adding them together to define the loadings from the sub-basin ^[32].

The flowchart in Figure 1 shows the optimization modeling approach for the SWAT model. The first step is to simulate the same management scenario for all HRUs in the SWAT interface. Several output files are created in every SWAT simulation. The proposed optimization modeling framework uses two files created by SWAT, namely output.std and ouput.hru. While the summary output file (output.std) provides average crop values for each HRU and their corresponding yield (kg/ha) and biomass (kg/ha), the HRU output file (output.hru) contains summary information for each HRU in the watershed. For a full description of all variables see the SWAT Input/Output Documentation version 2012 ^[32].

After the simulation is executed, the output.std and output.hru files are used to create an optimization table. This is performed by running a script coded in MATLAB, called "Optimization Tables". This code requires to define four variables: (1) HRU numbers; (2) management code (MC); (3) HRUs included in the optimization study; and (4) three different response variables. The number of HRUs in the simulation allows the code to identify the location of average values for biomass in the output.std file. Since every management scenario simulated will create a new optimization table, the MC allows identifying the management scenario simulated in SWAT when possible solutions are evaluated by the MOEA.

The HRUs included in the optimization study allow excluding those HRUs that could not be optimized due to their natural characteristics (e.g., water or urban HRUs). The three different response variables allow the code to identify the location in the output.hru file of the SWAT output values desired to be minimized. The MP, biomass, and response variables are saved in an optimization table later used in the MOEA script. This process is repeated until all the management practices in the study are executed in SWAT, and their output variables are extracted and saved in the optimization tables. The next step is to run the MOEA script, also coded in MATLAB. The MOEA starts by creating a random set of land-cover scenarios using the optimization tables. This initial set is called initial population. Next, the iterative process of the MOEA continues until reaching a predefined number of iterations, called generations. The final solutions are presented as a set of Pareto optimal solutions, and one solution is suggested based on its proximity to the ideal point [1,0,0,0].

3.2.   Problem formulation

A general formulation to quantify all possible landscape scenarios is based on the total number of management practices and the total number of Hydrologic Response Units (HRUs). Mathematically, this can be formulated as follows (Eq 1):

where X = the total number of management practices, and HRUs = HRUs _total − HRUs _water − HRUs _urban (excluding water and urban HRUs). For instance, if two management practices were to be evaluated in a watershed divided into three HRUs, the total number of combinations would be 2³ = 8 and each combination should be evaluated to determine the combination with greater environmental performance. Testing these combinations may seem achievable. However, numerous management practices can be tested in a much broader space. In hydrological simulations, watersheds are divided into several HRUs, where for a simulation considering 50 HRUs, and six management practices, the total possible combinations would be 6⁵⁰ ≈ 8∗10³⁸. Manually testing these combinations to obtain better environmental performance scenarios, even with simulation technology will be extremely time consuming and practically unfeasible.

3.3.   Multi-objective evolutionary algorithm (MOEA)

Multi-objective evolutionary algorithms use a population-based search. A population is a collection of individuals representing a possible solution to the multi-objective problem, also called chromosomes. In order to demonstrate the functionality of the proposed multi-objective optimization model, two decision alternatives were considered. Such decision alternatives generate potential individuals or landscape scenarios, which are evaluated for optimal agricultural landscapes. These combinations are listed in Table 1. The land cover considers three different perennial grasses; switchgrass, eastern gamagrass, and johnsonsgrass. The amounts of nitrogen applications are 30, 60, and 90 Kg/ha.

In this approach, a chromosome represents a possible landscape scenario that includes all HRUs in the SWAT simulation. Each chromosome is encoded to evaluate three response variables, which are SWAT average output values extracted from the output.hru file. Furthermore, the objective function also maximizes biomass yields, whose average values are extracted from the output.std file. In this context, the multi-objective optimization problem is formulated as follows (Eqs 2 and 3):

where: BY = average annual biomass yield, R1 = average annual first response selected, R2 = average annual second response selected, and R3 = average annual third response selected.

The multiobjective optimization technique applied in the proposed optimization modeling framework was developed by Taboada and Coit ^[30], which adjusts various characteristics from numerous metaheuristic methods to achieve quality approximations to global optimal solutions. Figure 2 shows the MOEA flowchart. The algorithm starts by creating a random initial population composed of many individuals. These individuals are possible solutions representing different landscape scenarios. Every individual's objective function is evaluated [BY, R1, R2 and R3], then the three fitness functions based on distance, dominance count, and an aggregated fitness metric are obtained. The aggregated results rank the entire population allowing the best-fitted individuals to populate the succeeding generation. The algorithm stops with a predefined number of generations or when the quality of the solutions reaches a steady-state. Finally, a recommended solution is selected from the Pareto set based on the proximity to the ideal vector [1,0,0,0].

3.3.1.   Initialization

The initial population is randomly created to ensure a diverse population for an effective search space. The initial set of individuals (population size) is a fixed parameter that must be defined at the first stage of the algorithm and remains constant in every iteration.

The configuration of an individual's genes, or chromosome encoding, distinguishes one individual from others. The number of HRUs determines the number of genes of an individual. The management practice code is specified in each gene. Therefore, the information for each HRU corresponds to each of the individual's genes. An individual and possible solution, as well as its chromosome encoding, is illustrated in Figure 3. Every square represents a gene. The HRU number is represented by the first numerical figure in the upper center of each gene. As previously stated, the number of genes represents the number of HRUs. Therefore, in this example, there are 15 HRUs. The number in the lower center of each gene is the management practice code listed in Table 1. In this example, the land cover in HRU 1 will change to eastern gamagrass, and 30 Kg/ha of nitrogen will be applied. Successively, the land cover in HRU 2 will change to switchgrass, and 60 Kg/ha will be applied, etc.

A possible initial population is illustrated in Figure 4. In this example, eight individuals are randomly generated, and every individual contains 15 genes.

3.3.2.   Evaluation and fitness assignment

Every individual is evaluated according to the fitness functions, and the best fitted individuals will populate the following generation. The concepts of Pareto dominance and population diversity are considered in the algorithm. In this regard, the dominated individuals are not considered for reproduction, and only the non-dominated individuals are considered for the following generation. Additionally, population diversity is attained by assigning higher fitness to those solutions that are far away from other solutions. These two criteria are evaluated according to different fitness metrics: fitness metric 1-distance-based, f₁(i); and fitness metric-2: dominance count-based, f₂(i); Before calculating the fitness metrics, every objective's result is normalized using Equation 4, to avoid discrepancies in units using the following equation (Eq 4).

where:${f}_{i}\left(x\right)$ is the value in the nondominated set, ${f}_{i}^{min}$ is the minimum value in the non-dominated set and, ${f}_{i}^{max}$ is the maximum value in the nondominated set.

The normalized values are used to calculate fitness metrics 1 and 2. the fitness metric 1, uses the Euclidean distances between the solutions, the sum of the distances from each solution to the rest of the solutions is obtained, and the maximum and minimum value of all the sums is calculated. These values are normalized and used to rank the individuals in order from the largest to the smallest. The fitness metric 2 is based on the dominance count concept, and it aims to obtain proximity to the true Pareto front. This metric assures that those individuals with the highest dominance count will have a higher ranking than those with the lowest dominance. Finally, the third fitness metric used is the aggregated fitness metric, f_a(x). The aggregated fitness metric is the result of the sum of fitness metric 1 plus fitness metric 2 [${f}_{a}\left(x\right) = {f}_{1}\left(x\right)+{f}_{2}\left(x\right)]$. It aims to weigh both metrics equally. At this point it is assumed that solutions with higher aggregated fitness value are solutions that are closest to the true Pareto front and farther away from other solutions. These values are used to rank the individuals in order from the highest to the smallest value obtained and then, they go through selection and crossover (Taboada and Coit ^[30]).

3.3.3.   Selection and crossover

In every iteration, a fraction of the non-dominated solutions found in the current generation will survive into the succeeding generation. The remaining spots in the succeeding generation will be filled by mixing current non-dominated individuals creating new possible solutions, a process called reproduction. This research considers tournament selection for the reproduction process. Particularly, two random individuals are selected to compete against each other. The most fitted individual is chosen to be parent 1. The same concept is used to select parent 2. The first set of parents selected will produce two new individuals in a process called crossover.

While reproduction can be attained with numerous crossover types, employing a suitable method that satisfies the chromosome encoding of the problem will determine its effectiveness. For this research, one crossover point is selected at random to exchange genes between parents. The process in which the crossover is executed is shown in Figure 5. Two segments of the chromosomes are divided by one random point. Offspring 1 is created using the first segment of parent 1 and the second segment of parent 2. Similarly, offspring 2 is created using the first segment of parent 2 and the second segment of parent 1. This process ensures that the order of the HRUs is maintained in the chromosome. The unoccupied places after elitism will be filled using this reproduction methodology.

3.3.4.   Mutation and stopping criteria

A slight mutation chance prevents falling into a local optimum. When mutation occurs, four random genes switch places with each other. Figure 6 shows this process. In this example, offspring 1 and 2 exchange the genes in HRUs 3, 6, 12, and 14.

The iterative process can end in different ways, including reaching a satisfying solution, selecting a certain number of iterations, or detecting a steady-state system. However, this approach uses a variable number of generations to end the MOEA's iterative process.

4.   Results

The SWAT example dataset (Lake Fork Watershed) was used to test the MOEA. This dataset is included in the download files of the ArcSWAT interface ^[33]. The latest ArcSWAT version stores this dataset in the directory C:\SWAT\ArcSWAT\Databases\Example1. The functionality of the MOEA is demonstrated by executing the landcover scenarios (Table 1) in SWAT and extracting the output variables in each scenario. All parameters have been set without any calibration, as described in Winchell ^[33].

4.1.   Lake Fork Watershed

The input dataset provides raster datasets for the Lake Fork Watershed into the Albers Equal Area projection with a resolution of 100 × 100 m: a Digital Elevation Model (DEM), a DEM mask (amask), a land-cover grid differing between six classes, and a soil map (U.S. general soil map STATSGO). The subwatershed is divided into 21 sub-basins, which we further divided into 25 HRUs. These are past (PAST), range-grasses (RNGE), and water (WATR) HRUs. Weather data was simulated using one Cooperative Observer Program (COOP) weather station. This input data enables the simulation to run for a period between January 1st, 1902, and January 1st, 2100. In our demonstration, we set the simulation period from January 1st, 2020, to January 1st, 2025.

The multi-objective optimization methodology used to solve this example considers maximizing biomass production (BY, kg/ha) while minimizing sediment yield (SYLD, t/ha), organic nitrogen yield (ORGN, kg/ha), and organic phosphorous yields (ORGP, kg/ha). The Pareto-set of optimal solutions obtained using a population size of 1000 and 1000 generations consists of 11 non-dominated solutions, all of which are suitable compromise solutions between objectives without degrading any of them. Figure 7 displays the three-dimensional representation of these solutions. However, in the case of the simultaneous optimization of these four objectives, in which one objective is maximized, and the rest are minimized, the suggested solution is selected considering the proximity of the normalized objectives to the ideal vector ${Z}^{ideal} = \left[1 \ 0 \ 0 \ 0\right]$.

The algorithm was fully coded in MATLAB and run on an HP computer, with an Intel® Core™ i5 6200U CPU processor, operating at 2.30 GHz 2.40 GHz and 8 GB of RAM. The computational time to evaluate these landscape scenarios and to obtain a suggested solution was 1120.3 seconds.

The suggested solution obtained under these settings is displayed in Table 2. The management practice code identifies the landscape scenarios as described in Table 1. The succeeding columns correspond to the outputs obtained with the management scenario evaluated in each HRU. For instance, the land-cover change in HRU 1 to eastern gamagrass using 60 kg of N/ha, which corresponds to management practice code 2.6 in Table 1, produces 6,329.9 kg/ha of biomass, 3.449 metric tons/ha of sediment, 7.832 kg of N/ha of organic nitrogen yield, and 0.95 kg of P/ha of organic phosphorous yield.

Finally, the spatial location of the suggested solution in Table 2 is graphically represented in the Lake Fork Watershed map in Figure 8.

5.   Conclusions

This work presents a novel MOEA, which uses outputs from the SWAT model to determine the optimal spatial placement of management practices at the watershed scale. This model was coded in MATLAB, and it was scripted to be suitable for any SWAT project. The application of the proposed MOEA in properly calibrated and validated SWAT projects can assist researchers and stakeholders in understanding the environmental impacts of various management practices, particularly for those involving biofuel feedstock production. The functionality of this optimization framework was demonstrated for perennial grasses instead of regional crops in the Lake Fork Watershed (SWAT's example data set). Specifically, this demonstration evaluated switchgrass (Panicum virgatum L.) eastern gamagrass (Trip-sucum dactyloides L.) and johnsongrass (Sorghum halepense (L.), Pers.) each with three levels of nitrogen application considering the minimization of sediment yield, organic nitrogen yield, and organic phosphorus yield while maximizing biomass production.

Nevertheless, there are significant limitations in our approach. First, while our model maximizes biomass yields, most landowners select crops according to market prices and crop cultivation experience. Even when our aim is focused on environmental impacts, a market price optimization variable should be considered in future research to make the suggested solutions more attractive for stakeholders. Since market prices are not steady, and they are not considered in the SWAT model, a new approach is needed to satisfy the characteristics of the proposed MOEA. Second, to properly compare and select optimal management systems, a ranking weight system should be included to assign weights to intervals of crop yields based on market prices and optimize these ranks instead of yields. Third, this framework requires the output.std and output.hru files after management scenarios are evaluated in the SWAT interface; it may be more efficient if management scenarios were developed directly in MATLAB. These limitations can be approached with a Graphical User Interface (GUI) in which SWAT users are allowed to develop their landscape scenarios, management practices, and select optimization criteria. However, this approach was not considered in this work as the main purpose of this paper is the establishment of a multi-objective optimization framework for optimal spatial placement of agricultural management practices using the SWAT model. Future goals are to develop an open-source GUI that can handle these current limitations to exploit the capabilities of the proposed MOEA framework.

Acknowledgments

This work has been supported by the National Institute of Food and Agriculture, under the Hispanic Serving Institutions Program, under award numbers 2015-38422-24112 and 2016-38422-25542. Any opinions, findings, conclusions, or recommendations expressed are those of the author(s) and do not necessarily reflect the view of the U.S. Department of Agriculture.

Conflict of interest

The authors declare no conflict of interest.