Powering sustainable development within planetary boundaries

The concept of planetary boundaries identifies a safe space for humanity. Current energy systems are primarily designed with a focus on total cost minimization and bounds on greenhouse gas emissions. Omitting planetary boundaries in energy systems design can lead to energy mixes unable to power our sustainable development. To overcome this conceptual limitation, we here incorporate planetary boundaries into energy systems models, explicitly linking energy generation with the Earth's ecological limits. Taking the United States as a testbed, we found that the least cost energy mix that would meet the Paris Agreement 2 degrees Celsius target still transgresses five out of eight planetary boundaries. It is possible to meet seven out of eight planetary boundaries concurrently by incurring a doubling of the cost compared to the least cost energy mix solution (1.3% of the United States gross domestic product in 2017). Due to the stringent downscaled planetary boundary on biogeochemical nitrogen flow, there is no energy mix in the United States capable of satisfying all planetary boundaries concurrently. Our work highlights the importance of considering planetary boundaries in energy systems design and paves the way for further research on how to effectively accomplish such integration in energy related studies.

Electricity transmission losses due to imports from state to state ' ' , Electricity exports from state to state ' , Electricity imports from Canadian region to state , Electricity transmission losses due to imports from Canadian region to state

Methodology and data sources: ERCOM-PB mathematical formulation
This section describes the modified mathematical formulation of the Emissions Reduction Cooperation Model (ERCOM) 1 that incorporates planetary boundaries, henceforth called ERCOM-PB.
We also highlight the data sources used in the original ERCOM (Supplementary Table 1). The work by Galán-Martín et al. 1 offers further details about ERCOM's assumptions and data sources.
Here, we extend the original ERCOM 1 by incorporating additional equations and constraints that include planetary boundaries into the model. In this section, we describe ERCOM-PB, which contains blocks of equations characterizing environmental burdens, planetary boundaries, loadmeeting constraints and equations required to assess the cost of electricity generation. ERCOM can be solved following a cooperative and non-cooperative mode between states 1 , yet for simplicity, this contribution deals only with the full cooperative mode. Thereby, the results shown in this paper represent an upper bound on the economic and environmental benefits that could be realized. Only the equations and constraints related to the full cooperative mode are described here; the interested reader is referred to the work by Galán-Martín et al. 1 for further details on cooperation-related equations.

Greenhouse gas emissions constraints
The United States (US) has pledged to reduce its greenhouse gas emissions by 26

Load meeting constraints
Electricity generation is limited by the generation potential available for each technology and state as follows: where is the electricity generation potential for technology in state . The constraint (SE3) is , applicable for all technologies except for those that share the same energy feedstock, namely coal and coal with Carbon Capture and Storage (CCS), natural gas with and without CCS as well as biomass and Bio-energy with CCS (BECCS). These technologies are modeled as follows: where the sets , and represent coal-, natural gas-and biomass-powered technologies, respectively. Constraints (SE4), (SE5) and (SE6) ensure that technologies that share the same tradable energy feedstock do not exceed certain sensible limits on the generation potential in a given state.
Such tradable energy resources are also bounded by their national availability, as stated in constraints (SE7-9).
Installed capacity and electricity generation, both standard and backup, are linked as follows: where is a continuous variable that represents the standard installed capacity for technology , in state , is also a continuous variable that provides the backup installed capacity for , technology in state , denotes the capacity factor used to link the capacity to the generation , with technology in state and is a scalar that represents the total number of hours of operation in a year. The capacity factor represents the ratio between the actual and potential electricity output.
This parameter limits the amount of electricity generated by each technology by considering those time periods in which the power plant is not operating. While the standard electricity generation is bounded by an inequality constraint (SE10), the backup electricity generation is modeled via an equality constraint (SE11) that ensures the grid reliability.
According to the US Energy Protection Agency, new nuclear power plants shall not be built 3 .
This policy is modeled as follows: where is a parameter denoting the current installed nuclear capacity.

,
The integration of more intermittent energy resources into the grid inevitably decreases its reliability. Many studies dealt with this modeling issue. One way to model this is to assign each technology an inertia potential, which is lower for intermittent technologies than for those that are not intermittent 4 . Thereafter, the grid reliability is maintained by setting a global system inertia potential as a lower bound in the model, which should be met by the grid 4 . Similarly, we use the backup generation methodology to ensure that for every unit of electricity generated from a nondispatchable resource, a certain proportion has to be generated with a dispatchable resource [5][6][7] . Both methodologies serve the same purpose, namely to maintain the reliability of the grid. Furthermore, both methodologies indirectly increase the marginal cost of production, particularly if the dispatchable technology chosen to support the grid is more expensive than the non-dispatchable one.
The reliability of the grid is modeled using equation (SE13), where is a scalar denoting the capacity share needed by dispatchable technologies for every unit of a non-dispatchable technology and is the set of non-dispatchable technologies. We also ensure that non-dispatchable technologies cannot generate backup electricity using equation (SE14).
Electricity can only be traded between neighboring states and Canadian regions as follows: , where represents electricity being imported by state from state , is the We also model the total electricity losses due to transmission as follows: where is the total amount of electricity received at the final destination from state and , ' ' is the total amount of electricity losses due to electricity transmission from state to state , , ' ' is the total amount of electricity received at the final destination from the neighboring , Canadian region and is the total amount of electricity losses due to transmission from , the neighboring Canadian region to state .
The total amount of electricity losses is modeled as a function of the distance and the amount of electricity being transmitted between regions as follows: where is the distance between state and neighboring state , is the distance , ' ' , between state and neighboring Canadian region and is a scalar that links the electricity losses to the transmission distance.
To be consistent with the current energy market structure, total electricity imports from neighboring Canadian regions and states cannot exceed a given share of the total demand as follows: where is a parameter that represents an upper bound on the total electricity demand that can be met by imports and is the total electricity demand in state .
To prevent states from acting as transmission nodes (i.e., import more electricity than the amount consumed to ultimately sell it), we constrain the total imports by each state by its total electricity demand (SE22).
The domestic electricity generation plus the imports and minus the total exports must satisfy the total demand for that state adjusted by a reliability factor as follows: where is the demand satisfaction factor required to maintain the grid reliability forcing the system to deliver a reserve margin above the electricity demand.

Economic objective function and constraints
The objective function seeks to minimize the total cost of electricity generation as follows: where is the total cost of electricity supply to the US and is the cost of electricity generation in state . The cost of electricity generation is structured as follows: where is the installed capacity cost in state , is the annual fixed cost in state , is the variable operating cost in state and is the electricity import cost from the neighboring Canadian regions to state . The installed capacities of both standard and backup technologies determine the state capital cost as follows: where is a parameter that represents the unitary capital cost of technology in state . The fixed , capacities of both standard and backup technologies determine the state annual fixed cost as follows: where is a parameter that represents the unitary annual fixed operating cost of technology in , state . The variable costs of both standard and backup technologies determine the state variable operating cost as follows: where is a parameter that represents the variable cost of technology in state . The total , electricity imports from neighboring Canadian regions provide the total electricity import cost as follows: where is a scalar that represents the Canadian unitary selling price of electricity.

Planetary boundaries constraints and objective function
Our work characterizes planetary boundaries and includes them into energy systems models. Each Earth-system process required to assess the performance of the grid in terms of planetary boundaries is modeled as follows: where represents the total performance of each Earth-system process linked to planetary boundary and is a parameter denoting the total environmental burden linked to planetary The US power sector share of the safe operating space linked to each planetary boundary connected to the Earth-system processes considered in this work are obtained as follows: where is the assigned share of the safe operating space to the US power sector, is the and only quantified in the other solutions (i.e., not optimized), is obtained as follows: where is the summation of the normalized transgression of each planetary boundary by its corresponding US power sector share of the safe operating space. Equation (SE37) is, therefore, the objective function of the minimization problem in the planetary boundaries solution (S3).
The following models are solved in the main manuscript to obtain three solutions namely S1, S2 and S3:  To compute solution S1, which corresponds to the US 2012 default developments in the power sector to meet the 2030 electricity demand, the share of each technology is fixed to its 2012 level and the demand is projected to 2030.
 To compute solution S2, which represents the least cost solution that meets the US commitment to the Paris Agreement, we minimize equation SE24, subject to constraints (SE1-SE29).
 To compute solution S3, which corresponds to the system minimizing the transgression of planetary boundaries at minimum cost, we minimize equation SE37, subject to constraints (SE3-SE37). We then take the value of the total normalized transgression, i.e., in equation To be more precise, we aim to quantify the likelihood of the Paris Agreement mix, solution (S2), being less expensive than the Business as Usual (BAU) mix, solution (S1). Here, we assume that LCOE values follow a uniform distribution where the ranges are tabulated in Supplementary Table 3 and perform a post-optimal analysis (i.e., optimal mixes are fixed to solutions S1 and S2 and then the LCOE values are varied between the current to future levels). The results are shown in Supplementary Figure 3 and discussed in Supplementary Note 7.

Uncertainty analysis approach corresponding to the life cycle inventory entries (results shown in the main body of the manuscript)
To assess the impact of the uncertainties involved in the life cycle inventory entries used to evaluate the performance of Earth-system processes on planetary boundaries, we generated 100 scenarios and solved the optimization model (i.e., ERCOM-PB) for each of them separately. Each scenario corresponds to a different materialization of all the uncertain parameters, i.e., life cycle inventory entries, simultaneously generated by applying Monte Carlo sampling on the probability distributions of the inventory parameters. This number of scenarios meets the Law and Kelton's test 9 , that is, it satisfies a confidence level of 95% with a relative error below 5% defined on the optimal objective function value. In Figures 2, 3, 4 and 5 in the manuscript, we show the results of the uncertainty analysis as error bars, where each error bar represents one standard deviation. We note that we apply the uncertainty analysis to the three solutions analyzed in the work (i.e., S1, S2 and S3).
where is the geometric standard deviation of the lognormal cumulative distribution function for The standard deviation of the uncertain parameter's natural logarithm is obtained as follows: while the mean of the uncertain parameter's natural logarithm is computed from its expected where is the expected value of the uncertain parameter (i.e., arithmetic mean), retrieved uncertain parameters, so a different seed is applied to sample each life cycle inventory entry separately. We applied the same approach to calculate the random life cycle inventory entries associated with electricity generation in neighboring Canadian regions.

Uncertainty analysis approach corresponding to the levelized cost of electricity (results shown only in this supplementary information and discussed in detail in Supplementary Notes 6 and 7)
We use two sets of LCOE values to (i) asses the uncertainty associated with future LCOE on meeting seven planetary boundaries concurrently, i.e., solution S3, (Supplementary Table 2) and (ii) quantify the impact of learning curves associated with the LCOE on the likelihood of the Paris Agreement, solution (S2), being more expensive than the BAU, solution (S1), (Supplementary Table 3 we performed a post-optimal analysis in analysis (ii) where we fix the mixes corresponding to solutions S1 and S2 and then vary the LCOE values reported in Supplementary Table 3 We assume that the LCOE follows a uniform distribution. To generate the Monte Carlo samples following a uniform distribution, two parameters are needed, namely a lower and upper bound for every electricity technology (Supplementary Table 2 and Supplementary Table 3). The probability density function for the uniform distribution can be defined as follows: where is the random value, is the lower bound (minimum LCOE value) and is the upper bound (maximum LCOE value). Supplementary where , and are the average capital, fixed operating and maintenance and variable operating and maintenance costs, respectively, per unit of electricity generated with technology , is the cost adjustment factor corresponding to state and is the cardinal (i.e., size) | | of set that represents the states considered in the model.

Supplementary Note 4 Limitations and future work
The following limitations and assumptions of our work are acknowledged:  We acknowledge that our allocation methodology assumes static shares of planetary boundaries.
Nonetheless, those shares could be dynamic in nature based on the 'elasticity' 13 of each sector, that is, on its ability to reduce its burdens while maintaining a given level of outputs. For example, the power sector has less technological resilience towards reducing the biogeochemical Nitrogen Moreover, population and gross value added vary over time and, therefore, the shares derived from such indicator should be updated periodically.
 According to the PB-LCIA methodology 14 , freshwater withdrawal (widely available in life cycle repositories, such as ecoinvent 10,15 ) is used to characterize the environmental flows related to the planetary boundary on global freshwater use. Nonetheless, there could be a slight mismatch between freshwater consumption and freshwater withdrawal, since part of the withdrawn freshwater could potentially be recycled 16 . Unfortunately, freshwater consumption data for each electricity technology are currently unavailable.   Figure 2). Therefore, our economic insights are robust within the assumptions highlighted in Supplementary Table 2.

Likelihood of the Paris Agreement mix being more expensive than the business as usual mix: postoptimal analysis focusing on learning curves (analysis (ii) in Supplementary Note 3)
When future energy mixes are studied, their future LCOE values should be considered to derive sound insights. Typically, learning curves are used to quantify the LCOE reductions for renewable technologies that currently might not be competitive when compared to conventional ones. Many studies concluded that the future LCOE of some renewable technologies (e.g., wind onshore and geothermal) would be more competitive than those of conventional ones (e.g., coal power plants) 12,19 .
Here, we are interested in evaluating the extent to which the Paris Agreement mix (solution S2) is better in economic terms than the BAU mix (solution S1). To this end, the outcome of learning curves is varied by performing a post-optimal analysis on the mixes found. First, we fix the mixes found with the average LCOE values (third column in Supplementary Table 2). Thereafter, we vary the LCOE values in the ranges reported in Supplementary Table 3  The Paris Agreement mix (S2) is found to be less costly than the BAU mix (S1) in 95 out of the 100 scenarios (i.e., runs) considered (Supplementary Figure 3). This is due to the competitiveness of some electricity technologies in the future, primarily wind onshore and geothermal -both combined representing only 4% of the BAU mix compared to 32% in the Paris Agreement solution. Therefore, despite the uncertainty associated with the LCOE of renewable technologies in the future, the Paris Agreement mix remains less expensive than the BAU mix with a probability of 95%.  Table 2). Whiskers represent the range from 5% to 95%. The dashed white line represents the mean and the solid white line represents the median. The 'x' sign represents the last 1% of the results and the '-' sign represents the maximum and minimum values in the sample. The x-axis shows the implications of the uncertainty in the future levelized cost of electricity values on the cost of US electricity generation and imports (primary y-axis), as well as on the performance relative to downscaled planetary boundaries (secondary y-axis).
Supplementary Figure 3 Box plot with whisker summarizing the likelihood of the Paris Agreement solution (S2) being more expensive than the business as usual solution (S1) due to the uncertainty in the levelized cost of electricity induced by the consideration of learning curves (Supplementary Table 3). Whiskers represent the range from 5% to 95%. The dashed white line represents the mean and the solid white line represents the median. The 'x' sign represents the last 1% of the results and the '-' sign represents the maximum and minimum values in the sample. Total transmission distance between states and Canadian regions Short et al. 25 Demand satisfaction factor, 1.05 , NERL 26 Potential of electricity generation Total United States potential of electricity generation United Nations 27 Gross value added for the United States total economy in 2016 United Nations 27 Gross value added for the United States power sector in 2016 100

Supplementary
100-year global warming potential per unit of electricity generated Remaining technologies: ecoinvent 10,15 Life cycle inventory entry per unit of electricity generated World Bank Group 32 United States population in 2016 World Bank Group 32 World population in 2016 Steffen et al. 17 and Ryberg et al. 33 Full safe operating space for every planetary boundary Paris Agreement 2 Paris Agreement target in 2030 Fripp 34 Losses due to electricity generation, 0.62% for every 100 km U.S. Army Corps of Engineers 35 Cost adjustment factor for every state Supplementary Table 2 Future average levelized cost of electricity for each technology and uncertainty ranges. a We use the share of the additional cost caused by adding a CCS unit to a biomass power plant for every cost component reported by Cuellar and Herzog 23 and multiply the shares by the biomass power plant cost components reported by EIA 12 . We assume 19% reduction from the resultant current levelized cost of electricity of BECCS to the one in 2030 in line with EIA projections 12, 20 for biomass plants. b We assume 30% of the cost corresponding to the capital cost of the CCS unit to follow the range reported by Koelbl et al. 18 (±25% of the mean) and the remaining 70% of the cost to follow the range reported by EIA 12