Process design within planetary boundaries: Application to CO2 based methanol production

(cid:1) Incorporation of absolute sustainability criteria in process design. (cid:1) Integration of process simulation, surrogate modelling and planetary boundaries. (cid:1) Methodology to design processes consistent with the planet’s ecological capacity.


Introduction
The ongoing transition towards a more sustainable chemical sector calls for advanced decision-support tools to incorporate environmental criteria in the design of emerging processes and retrofit of existing ones.While designing greener processes is now receiving increasing interest, it is indeed a topic long considered in the chemical engineering literature, which covers multiple approaches based on a wide range of metrics.Notably, in a pioneering work, Young et al. introduced the WAR algorithm (Young et al., 2000), which quantifies the potential environmental impacts of a chemical process.Sikdar (2003) proposed a hierarchical approach based on 3-dimensional (3-D) indicators covering different sustainability aspects.This framework was further refined by Martins et al. (2007), who included 3-D and economic indicators into process design.Aliff Radzuan et al. (2019) developed the SUI metric, which combines economic, sustainability, and environmental aspects into a weighted index, and used it to compare design alternatives for the production of cyclohexane from benzene and hydrogen.Benoit et al. (2019) developed a tool for assessing the eco-efficiency of a process using different impact metrics.With a similar spirit, Pereira et al. (2018) developed the eco-efficiency comparison index to enable comparisons between alternative designs.
In recent decades, it has become clear that identifying environmentally superior technologies requires enlarging the scope of the analysis beyond the chemical plant.Hence, a proper environmental analysis of chemical technologies requires embracing impacts over the whole life cycle to avoid shifting burdens across echelons in the chemical supply chain (Algunaibet and Guillén-Gosálbez, 2019).In this spirit, life cycle assessment (LCA) (Kleinekorte, 2020) has become the prevalent method in the environmental assessment of chemical processes, either as a standalone tool or coupled with process optimization within the life cycle optimization framework (LCO) (Azapagic and Clift, 1999).The latter approach has found multiple applications in Process System Engineering (Pieragostini et al., 2012;Grossmann and Guillén-Gosálbez, 2010;Guillén-Gosálbez et al., 2019), including the areas of process design (Guillén-Gosálbez et al., 2008;Hugo et al., 2004), and supply chain optimization (Mota et al., 2015;You et al., 2012).
The main limitation of the plant-based methods and the LCA and LCO approaches employed in green engineering is that they rely on metrics that are hardly interpretable from a worldwide sustainability perspective.Notably, LCA indicators enable comparisons between process alternatives but provide limited insight into whether technologies are sustainable in absolute terms.LCO inherits this critical limitation, as optimal solutions minimizing LCA impacts are not guaranteed to be sustainable from an absolute sustainability viewpoint, i.e., when evaluated considering the finite capacity of the Earth.
Among the approaches available, the planetary boundaries (PB) concept, proposed by Röckstrom (Rockström, 2009;Steffen, et al., 2015), provides a robust framework to quantify absolute sustainability accurately.In essence, the PBs define a set of ecological thresholds on nine Earth system processes critical for the resilience of the planet.PBs on nine Earth system processes were put forward initially (Rockström, 2009), i.e., climate change, rate of biodiversity loss (terrestrial and marine), interference with the nitrogen and phosphorus cycle, stratospheric ozone depletion, ocean acidification, global freshwater use, change in land use, atmospheric aerosol loading and chemical pollution.These bounds were later refined in (Steffen, et al., 2015), which further developed the biogeochemical flows, biodiversity integrity, and novel entities PBs.
Transgressing any of these boundaries could trigger a set of deleterious events that could shift the current state of the Earth to a new unknown, and very likely less friendly for humanity, equilibrium state.Ongoing research on the PBs aims to define these critical thresholds more accurately, while their application to the assessment of chemical processes is still in its infancy.
The chemical sector is witnessing a paradigm shift toward more sustainable technologies (Guillén-Gosálbez et al., 2019;Bakshi, Jun. 2019), including the transition to renewable carbon, e.g., biogenic carbon and captured CO 2 , less energy-and water-intensive processes, and technologies closing the materials loop through the adoption of circular economy principles.In this context, assessing the absolute sustainability level of emerging technologies becomes essential to guide future experimental work and regulations more sensibly (González-Garay et al., 2019).
Bearing the above in mind, here we incorporate the PBs in the design of chemical processes to identify process flowsheets that are entirely consistent -environmentally speaking-with sustainable development.Process design is a fundamental problem in chemical engineering, often addressed by maximizing economic performance and imposing bounds on environmental impacts, quantified via simple metrics and/or constraints based on current regulations (e.g., CO 2 emissions, waste, etc.).Some works applied the LCO framework to minimize the life cycle impact of chemical processes (Guillén-Gosálbez et al., Feb. 2008;Hugo et al., 2004), yet even this approach cannot ensure that the design identified is environmentally sustainable.Similarly, standard pollution prevention methods for process design (El-Halwagi and Prevention, 1997) can help to reduce the emissions and waste of a chemical plant, yet they cannot guarantee that the final design is sustainable from a worldwide sustainability perspective.
In the derivation of our process design method, we capitalize on the application of the PBs framework to other engineering problems, namely the design of sustainable power mixes (Algunaibet et al., 2019), and fuels supply chains (Wheeler et al., 2020;Ehrenstein et al., 2020), as well as the evaluation of land-use strategies (Heck et al., 2018), laundry washing (Ryberg, 2018), ammonia production (Samaroo et al., 2020) and the chemical sector (Galán-Martín et al., 2021), among others.We note that we share with other works the overall aim of embracing multiple sustainability criteria in the analysis, as done in the recently developed SOFTSCAPES framework (Mangili et al., 2019  Going well beyond standard green process design methods, here we link the performance of a chemical flowsheet with its impact on the planet using a recently proposed LCA-PB methodology coupled with process design tools.Notably, we formulate the process design task as an optimization problem that incorporates explicit LCA-based constraints linking the mass and energy flows of the process with the control variables of the PBs.The resulting bi-criteria optimization model minimizes then the total cost and the transgression of the PBs to identify process flowsheets that are entirely sustainable in absolute (environmental) terms.
To illustrate the capabilities of our approach, we apply it to green methanol production based on CO 2 hydrogenation (González-Garay et al., 2019;Hank, 2018).Green methanol has attracted increasing interest in the process modeling and catalysis communities, as it could be used as a green precursor for a wide range of chemicals via methanol-to-olefins and methanol-toaromatics processes (Tian et al., 2015;Samimi et al., 2018;Yang et al., 2019).However, green methanol requires large amounts of renewable energy to activate the inert CO 2 molecule, which raises concerns about its sustainability level.Here we shall design a CO 2 hydrogenation plant considering different hydrogen sources to optimize the methanol cost and its absolute sustainability level simultaneously.Our approach identifies the optimal operating conditions, equipment sizes, and hydrogen feed that minimize the transgression level of eight PBs and the total annualized cost.
The paper is organized as follows.First, we formally state the problem of interest, followed by the methodology and then the numerical results.Finally, we draw the main conclusions of the work and outline future research directions.

Problem statement
Process system engineering (PSE) spans multiple scales, from the molecular level to the enterprise level (Grossmann and Westerberg, 2000).Here we shall focus on process design to illustrate the capabilities of combining PBs with mathematical programming to optimize sustainable industrial systems.We are given a final product demand to meet, the final product specifications, and a set of unit operations and raw materials to cover this demand.The goal is to find the optimal design and operating conditions that minimize the cost and environmental impact while meeting technical and demand satisfaction constraints.

Mathematical formulation
The standard process design task can be formulated as a Mixed-Integer Non-Linear Programming (MINLP) problem, where continuous variables represent operating conditions, mass and energy flows and sizes of equipment units, and binary variables denote topological decisions (i.e., the existence of process units and connectivity between them).This MINLP for process design is here constructed using surrogate models trained with process simulations to predict the behavior of the chemical plant precisely.The original design problem can be mathematically formulated as follows.min x;y where w 1 x; y ð Þ; w 2 x; y ð Þ refer to the objective functions (i.e., economic and environmental), x and y denote the continuous and binary variables, hðx; yÞ and gðx; yÞ are the equality and inequality constraints, and LB; UB refer to the lower and upper bounds on the continuous variables, respectively.
As already mentioned, we here use surrogate models to build the MINLP for process design.The continuous variables are divided into independent (i.e., degrees of freedom, denoted by x I and y I ) and dependent (x O ), where the latter are computed by the surrogate (i.e., process simulation) once the independent variables have been fixed.The problem can then be reformulated as in Eq. (2).
where f 0 x I ; y I À Á represents the surrogate model, computed from the input (independent) variables, whose values are constrained within allowable lower and upper bounds (LB I and UB I , respectively).External equality and inequality constraints (i.e., h x I ; y I À Á and g x I ; y I À Á , respectively) may also be imposed on the independent variables.
We next focus on describing how the objective functions are quantified, with a strong focus on the environmental impact, which is the main contribution of this work.

Environmental performance: Planetary boundaries
The PBs represent a set of ecological limits defined on critical Earth systems that should never be exceeded to operate our planet safely (Rockström, 2009).PBs on nine Earth system processes were put forward, which define three regions for every Earth system, i.e., the safe operating space (SOS) (green zone), the uncertain zone (yellow zone), and the high-risk zone (red zone).The green zone would allow us to operate the planet safely for many years to come, while the uncertain and high-risk regions entail an increasing probability of catastrophic events that could shift the current equilibrium state of the Earth.
The PBs impose limits on a set of control variables defined for the Earth systems.In some cases, the control variables consider a natural background level, i.e., an amount of burden that is independent of anthropogenic activity.This means that although the PB would seem to give an ample margin, the anthropogenic activity must be more limited in order to not surpass the PB.The SOS then considers the difference between the value of the PB and the natural background level.For example, the PB of climate change is given by a CO 2 concentration of 350 ppm, i.e., concentrations below 350 ppm do not surpass the PB and are therefore within the SOS.However, the natural background level for this control variable is 278 ppm.Thus, the SOS regarding anthropogenic activities is obtained from this difference, i.e., 350 ppm -278 ppm, 72 ppm.Four PBs were already transgressed (i.e., climate change, biosphere integrity, biogeochemical flows, and land-system change), which reinforces the need to minimize the impact of industrial systems considering the carrying capacity of the planet.
In this work, we follow two steps to quantify the impact of a chemical process on the PBs, as described below.
Step 1: Downscaling the SOS to the process level First, we downscale the PBs to the process level.The PBs define global ecological budgets that should be shared among all economic sectors jointly.Hence, we first need to allocate part of this budget, the so-called safe operating space (SOS), to the process under study.This share of the SOS allows us to determine whether the process operates below or above its assigned ecological budget and, therefore, elucidate if it is or not sustainable.
There are multiple downscaling principles that can be applied to assign shares of the SOS to industrial systems (Ryberg et al., 2020).So far, the selection of downscaling principles lacks scientific consensus and remains controversial.In this work, we apply an egalitarian principle based on the Gross Value Added (GVA), which assigns a share of the SOS based on the GVA of the process relative to the GVA of the world, as follows: where a is the share of the SOS assigned to the system (e.g., methanol production), GVA FU refers to its GVA and GVA World is the GVA of the world.
Step 2: Calculating the impact on the control variables of the PBs The next step to quantify the PBs performance of a process is to link the life cycle feedstock requirements, emissions and waste of the process to its impact on the control variables of the Earth systems of the PBs.To this end, we first quantify the Life Cycle Inventory (LCI) of feedstock, emissions and waste following the standard LCA phase three (i.e., inventory analysis).The LCI is obtained by combining data of the foreground system (i.e., mass and energy flows associated with the process flowsheet) with data of the background system (i.e., activities supplying inputs to the main process, e.g., raw materials, electricity, steam).Background data can be retrieved from environmental databases such as Gabi (Volz et al., 1998) and Ecoinvent 3.5 (Wernet et al., 2019), as well as from technical reports and scientific literature.In contrast, data of the foreground system is provided by the process simulation model (e.g., Aspen-HYSYS, Aspen Plus, gPROMS).In essence, the LCI entries j (in the set J) can be computed from the raw materials, utilities and direct emissions of the main process as follows: where LCA RM j;r denotes the LCI entry (i.e., feedstock, emissions, or waste) per unit of mass of raw material r, LCA UT j;k denotes the LCI entry per unit of utility k, and LCA DIR j and LCA WASTE j are the natural flows directly exchanged between the process and the biosphere (i.e., direct emissions and waste).In the same equation, RM r denotes the consumption of raw material r, and UT k represents the amount of utility k consumed by the process (i.e., per unit of the chemical produced, taken as the functional unit for the LCA calculations).
We next link the LCI entries to the impact on the control variables of the PBs using the characterization factors recently developed by Ryberg (Ryberg et al., 2018), as shown in Eq. ( 5).
where J refers to the set of environmental entries (i.e., feedstock, emissions and waste) linked to each PB b, and CF b;j is the characterization factor associated with planetary boundary b and environ-mental entry j.IMP D b is the contribution to the transgression of PB b per functional unit, e.g., the yearly production of methanol.
Finally, we can calculate the transgression level attained by a design alternative in every PB as follows: where values of TL D b above one imply that the process transgresses its share of PB b, so it is unsustainable from the viewpoint of that PB.Values below one imply that the process operates within its share of the SOS, so it can be concluded that it is environmentally sustainable in that PB.
To quantify the transgression levels in the objective function, we define the transgression level excess (TLE) as follows: where in essence, we only penalize the transgression levels above one (relative to the downscaled SOS).This is because designs lying below the PB (green zone) should be deemed sustainable.We reformulate the transgression level function as follows (considering the share of the SOS previously computed following a specific downscaling principle): The environmental objective, defined as the mean transgression level (MTL), is computed as shown in Eq. ( 9), where we consider the total average transgression level across all the PBs.

MTL
In this work, we consider the following PBs: CO 2 concentration and energy imbalance, the biogeochemical flows of nitrogen and phosphorus, stratospheric ozone depletion, ocean acidification, land-system change, and freshwater use.The other PBs are omitted due either to the lack of suitable characterization factors (as in biosphere integrity), or to the lack of both appropriate thresholds and characterization factors (i.e., aerosol loading and novel entities).Note that our objective function does not rely on weighting factors because all the PBs are considered equally important.This is because transgressing any of them could compromise the resilience of the Earth.

Economic assessment
The economic performance is quantified via the total annualized cost (TAC), which is computed as follows: where ACCR stands for the annual capital charge ratio, and CAPEX; OPEX stand for capital expenditures and operational expenditures, respectively.Details on how to compute these terms can be found in standard text books (Towler and Sinnott, 2012) and are provided in the SM, in section S1.

Solution procedure using a surrogate model
We briefly describe next how the calculations are carried out.Further details are provided in the case study section.
Without loss of generality, we build the surrogate model of the original process simulation using a Bayesian Regularized Artificial Neural Network (BRANN).This choice is motivated by its robustness compared to standard artificial neural networks and its good performance in terms of overfitting (Burden and Winkler, 2009).The BRANN is built using the results of a set of simulation runs generated by varying the independent variables between their lower and upper bounds in Aspen-HYSYS.We generate a given number of data points (1000 in our case) via the Latin Hypercube Sampling (LHS).At each point, the MINLP model implemented in SYNHEAT (Yee and Grossmann, 1990) is run in order to optimize the heat exchanger network (HEN) associated with the process simulated in Aspen-HYSYS (i.e., to minimize the total cost of the HEN).As will be later discussed in more detail, the surrogate model considers eight outputs: the total cost of the inside battery limits (ISBL), the amount of chemical (methanol) produced, the requirements of heating, cooling and electricity, the direct emissions (CO 2 and CO), and the amount of waste (waste water generated).Therefore, the data to produce the ANN comprises an input matrix of c Â inputs and c Â outputs, where c is the number of points from the original sampling that converge to a solution.
The surrogate model can then be optimized using standard multi-objective optimization algorithms, such as the e-constraint method (Mavrotas and Florios, 2013).This algorithm keeps one objective as the main objective while transferring the others to auxiliary constraints that impose bounds on them.These epsilon bounds are then modified as iterations proceed, identifying in each run a different Pareto point.The computational advantage of minimizing the PBs transgression level (as opposed to minimizing several LCA metrics simultaneously) is that it leads to bi-criteria models (i.e., cost vs transgression level) that are easier to handle.

Software implementation
The overall solution approach combines a palette of software packages, i.e., Python, GAMS, Matlab, and Aspen-HYSYS, to carry out the calculations.The bulk of the calculations are performed in Matlab and Aspen-HYSYS, with GAMS solving the HEN design problem based on the MINLP developed by Yee and Grossmann (Yee and Grossmann, 1990), implemented in the SYNHEAT code.
The methodology is outlined in Fig. 1.In essence, Python sends input data to Aspen-HYSYS (i.e., values of the independent values generated via sampling within allowable limits), and reads the simulation output (i.e., values of the dependent variables).Afterward, part of this output, i.e., process streams that can be integrated and their target temperatures, are sent to GAMS where the HEN is optimized using SYNHEAT.The results of the process simulation model, together with the SYNHEAT results, are sent then to Matlab, where the neural network is built.Finally, the neural network is optimized with an interior point algorithm (Byrd et al., 1999;Byrd et al., 2000;Waltz et al., 2006) (implemented in the fmincon function of Matlab) and considering as additional constraints the bounds on the independent variables used to build the surrogate model.

Case study
Methanol is currently produced from syngas generated from natural gas, therefore relying on fossil carbon.Some alternatives are being explored using biogas (Vita, 2018), shale gas (Ehlinger et al., 2014) or captured CO 2 and hydrogen as feedstock.Shale gas is still based on fossil carbon, so it should be ultimately ruled out to close the carbon loop in the chemical industry.The CO 2 hydrogenation route, which holds good promise in the transition to renewable carbon-based chemicals, gives rise to a set of alternative pathways differing in the CO 2 and hydrogen sources.CO 2 can be captured mostly from power plants based on fossil fuels (natural gas and coal) or directly from the air, where the latter route is the only one entirely consistent with circular economy principles.
Likewise, hydrogen can be synthesized from natural gas via SMR with CCS, from biomass through gasification, and also from the water via electrolysis.Biomass gasification is constrained by its availability, strongly linked to geographical limitations and the competition for land with food and electricity generation.Biogas presents a high variability depending on the gasified feedstock.Furthermore, water electrolysis could be powered by solar, wind, nuclear, and hydropower.Nuclear is affected by social acceptability issues, while the plant's location limits the availability of hydropower.Solar and wind are both suitable renewable sources for water splitting, yet wind shows a lower carbon footprint.A more in-depth analysis of these scenarios can be found in the works by Ioannou et al. (2021) and Galán-Martín et al. (2021).
Hence, we apply our method to design a methanol plant that consumes captured CO 2 and hydrogen as feedstock, where the former is obtained from direct air capture facilities (DAC) powered with natural gas and electricity from the current mix, and the latter from either steam methane reforming (SMR) with carbon capture and storage (CCS) or water splitting powered by wind.This carbon Fig. 1.Flowchart of the proposed computational approach for sustainable process design within PBs. capture and utilization (CCU) process creates economic value from the CO 2 instead of storing it in geological sites.The scope of the process is summarized in Fig. 2.
Fig. 3 shows the process flowsheet based on the works by Van-Dal and Bouallou (2013) and González-Garay et al. (2019).The simulation was implemented in Aspen HYSYS v11, using the thermodynamic packages Peng-Robinson and NRTL -Ideal.The hydrogen and CO 2 inlet streams are compressed and then mixed before being fed to the reactor, where the reaction, highly exothermic, occurs.The reactor outlet is fed to two flashes and a distillation column that separate the methanol product from the unreacted reactants, CO, and the water byproduct.
More precisely, we consider a plant with a constant CO 2 feed of 2,000 kmol/h available at 25 °C and 1 bar.These conditions assume the existence of a DAC facility integrated with the methanol production plant.This stream is compressed to reach the reaction pressure, modeled as an independent variable to be optimized.The use of cooling mid-compression is needed due to the high temperatures reached in the compressor, while a compression train is required to maintain the compression ratio below three since the desired final pressure is in the range of 45-55 bar.
Regardless of its provenance, the hydrogen stream is available at 30 bar and also needs to be compressed to reach the reaction pressure.However, since the compression ratio lies in the range of 1.50-1.83,there is no need to implement multiple compression stages (Luyben, Dec. 2011).
After being compressed, the two gases are mixed with the recycled stream and the resultant stream is heated to reach the desired reaction temperature.Two main reactions occur in the reactor.The first one is the CO hydrogenation to produce methanol (R1), which is accompanied by the water-gas shift reaction (R2), leading to the global reaction (R3), as shown below.
These reactions take place in a plug flow reactor (PFR) with a Cu-ZnO-Al 2 O 3 catalyst .We consider the kinetic model developed by Bussche and Froment (1996).We note, however, that other alternative kinetic models could have been used instead (Huš et al., 2017;Pavlišic ˇet al., 2020;Slotboom, 2020).The reactor outlet is cooled to 35 °C and separated in a flash unit.Most of the unreacted hydrogen and CO 2 are recycled and recompressed to reach the reaction pressure.Part of the recycle stream, mainly containing CO 2 , is purged (e.g., from 0.1% to 5%).The bottom stream of the flash is depressurized to 2 bar and sent to another flash unit before entering a distillation column with a partial condenser, operating at a head pressure of 1 bar.The distillate product is methanol with a 99% purity and at 64 °C, while the overhead, containing residual traces of CO 2 and methanol, is sent to the purge.The bottom product is wastewater at approximately 109 °C, which will be treated in a wastewater plant.More details of the simulation results are shown in the supplementary material, section S3.
The HEN design is optimized based on the streams data of the Aspen-HYSYS model.The utilities considered are high-pressure steam (40 bar, 250 °C) and cooling water (25 °C to 30 °C).The minimum temperature difference is 10 °C.We assume that the compressors operate with electricity from the grid.Another option would be to use part of the heat available in the reaction outlet to generate electricity via steam turbines.Considering that the future mix will get decarbonized, it seems more sensible to use the heat available in the reactor to eliminate the need for steam in the plant, e.g., in the column's reboiler, rather than to generate electricity to power the compressors.
The purge contains direct emissions, mostly CO 2 and CO, which are appropriately accounted for in the quantification of the impact on the PBs.

Mathematical formulation
We defined six independent variables (degrees of freedom) in the optimization, i.e., volume, temperature and pressure of the reactor, molar flowrate of the inlet hydrogen, purge ratio, and reflux ratio in the distillation column.These independent variables are varied between a lower and an upper bound, as shown in Table 1.The amount of CO 2 is fixed at 2,000 kmol/h.
The provenance of hydrogen does not affect the simulation model.Therefore, the ANN is built considering six inputs, while the optimization algorithm optimizes seven variables.Hence, in the optimization step, we disaggregate the molar flow of hydrogen into two variables, i.e., SMR with CCS hydrogen and electrolytic hydrogen.This choice of the hydrogen source provides the optimization algorithm with additional flexibility to minimize the cost and the transgression level of the PBs.

Economic cost
The TAC is calculated using the design parameters of each process unit and data from Sinnott and Towler (Towler and Sinnott, 2012).
Regarding the variable costs of operation, we consider here the costs of all the mass and energy inputs, i.e., hydrogen and CO 2 feed, cooling water for the coolers and condensers, steam for the heaters and reboilers, and electricity for the compressors.We also take into account the treatment cost of the waste water.
In order to adjust the costs from any year to 2018, we consider a 1.58% inflation rate, except for the process units, where we use the CEPCI value.We apply the conversion factor from Euros to dollars, i.e., 1.1597 $/€.For the business as usual (BAU) technology, we consider a TAC equal to the market price, i.e., 410 € 2017 /t (Hank, 2018).
The operation time of the plant is 8760 h/yr.
The equations and parameters used to calculate the TAC are shown in the supplementary material, in equations Eq.S1 -S9 and in Tables S1 -S4.

Planetary boundaries
We account for the following main sources of impact: the hydrogen and CO 2 feeds, the electricity and heating requirements, and the direct emissions.The impact of the cooling water is approximated by that embodied in the electricity required to pump it through a closed refrigeration system.The LCI is quantified by combining the mass and energy flows exchanged between the plant and the technosphere and ecosphere with data from Ecoinvent (Table S5).
The share of the SOS assigned to methanol is computed based on the data in Table 2. Due to the fact that the GVA of methanol is not available as such, we assume that it is given by the world revenues linked to methanol production worldwide, which amounted to 140 million tons in 2018 (GlobalData, 2019).Our results, therefore, underestimate the transgression levels of methanol, as the GVA should (arguably) consider the net profit margin, company-dependent and unavailable in the public domain, rather than the total revenues.
The values of the SOS are taken from the works by Steffen and Ryberg (Ryberg, 2018;Steffen, et al., 2015), and are shown in Table 3.
Further details are provided in the supplementary material, in section S2 and in Tables S6 -S8.

Solution procedure
The methanol production process from hydrogen and CO 2 has been widely studied and simulated (Hank, 2018;Samimi et al., 2018;Pérez-Fortes et al., 2016;Abrol and Hilton, 2012).Here we shall optimize the methanol flowsheet by combining process simulation with optimization following a simulation-optimization approach in which an optimization algorithm interrogates the process simulation model.Optimizing process simulation models can lead to numerical problems due to convergence issues, which can be overcome to some extent by resorting to surrogate models.The use of surrogate models to replace complex process units in process design, or even whole processes, has been widely studied  in the field of chemical engineering (Gonzalez-Garay and Guillen-Gosalbez, 2018; Quirante et al., 2018;Bhosekar and Ierapetritou, 2018;McBride and Sundmacher, 2019).Potential surrogate models include artificial neural networks, kriging interpolation (Krige, 1951), and symbolic regression (Guimerà, et al., 2020).In this work, we use neural networks to build a surrogate model of an Aspen-HYSYS simulation of the methanol process, which is optimized using an interior-point algorithm.The optimal Pareto solutions are finally evaluated in the original simulation to assess their accuracy.All the calculations were implemented in an Intel Ò Core i9-9900 CPU @3.10 GHz computer.The sampling points were generated in Python 3.7.9,using NumPy 1.19.1 and pyDOE 0.3.8.Afterward, through a COM interface between Python and Aspen-HYSYS v11, these points were sent to the process simulator.The flowsheet solver in Aspen-HYSYS was stopped and re-initialized each time a variable was modified.The data for the HEN design was collected by Python and sent to GAMS 32.2.0 through a gdx file.The MINLP results were sent back to Python via another gdx file.Afterward, the data were processed in Python to determine the values of the output variables, which were used to build the neural network in Matlab R2020a.
The MINLP solved in SYNHEAT contains 222 variables and 232 constraints.It took 0.141 CPU s to solve it to a 0% relative optimality gap using the DICOPT solver on the said computer.The surrogate model, i.e., BRANN, predicts eight outputs (dependent variables) that allow quantifying the TAC and the average PBs transgression level.These are the total cost of the process units needed to estimate the inside battery limits (ISBL) expenditures, the amount of methanol produced, the heating and cooling requirements, the electricity input, the amounts of CO 2 and CO emitted, and the amount of water to be treated.
The optimization is performed using a constrained optimization with the default interior-point algorithm implemented in the fmincon function in Matlab (Byrd et al., 1999;Byrd et al., 2000;Waltz et al., 2006).The neural network does not extrapolate well; therefore, the limits imposed on the input variables during the training are also enforced as constraints in the optimization.The multi-objective optimization provides a Pareto frontier, whose points are evaluated in Aspen-HYSYS to check their feasibility and provide more accurate estimates.

Sampling and surrogate model
We generated 1000 points using the LHS sampling method inside the package pyDOE 0.3.8.From these 1000 points, 984 converged in Aspen-HYSYS.However, for six of these, the heat integration model failed to converge.
The dataset of 978 Â 6 inputs and 978 Â 8 outputs was sent to Matlab using a .matfile, where we implemented the neural network.The neural net was trained using the Bayesian regularization algorithm in the Matlab Neural Fitting app.90% of the 978 training points were used to train the neural network, while 10% were employed as testing points.The Bayesian regularization algorithm implements a specific validation method, so it does not require defining a separate validation set.The neural net implements 15 hidden neurons in one layer.
Fig. 4 summarizes the global performance of the neural net, while Fig. 5 provides further details on the quality of the fitting in each variable.Overall, R values closer to one indicate a better fit.Fig. 4 displays global performance values by expressing all the outputs using a common scale.Fig. 5 provides the R 2 together with the mean square error (MSE) and mean absolute error (MAE).
The surrogate model performs well on average, and also in almost all of the outputs (R 2 above 0.90), except for the total cost of the process units (R 2 of 0.71).This poorer approximation does not substantially affect the optimization, as the total cost is mainly dominated by the raw materials costs, not by the CAPEX term.The worse performance of the ISBL could stem from the fact that the HEN optimization problem is embedded in the sampling.From the predicted outputs, only the ISBL, heating and cooling needs are not obtained directly from the flowsheet.The heating and cooling needs, even though they are calculated from the HEN optimization, are practically equivalent to the target values (thermodynamic limits).The ISBL output combines variables from the simulation flowsheet with the SYNHEAT results, which may explain its lower accuracy.Even though the R 2 is worse than the others, it is important to notice that the MAE is around 5.47, which is deemed acceptable considering the scope of the analysis.

Results of the multi-objective optimization
Fig. 6 shows the results of the multi-objective optimization, evaluated in the surrogate model and in the Aspen-HYSYS simulation.
The gray points in the figure correspond to the neural network, while the light blue points were evaluated in the Aspen-HYSYS simulation.The dotted line links the two sets of points, i.e., surrogated -''true" value.The difference between both models lies in the range 2.80-4.56%for the economic performance and 2.81-23.72%for the environmental one.
As seen, there is a clear trade-off between economic and environmental performance.Starting from the minimum cost solution, reducing the transgression level requires increasing the cost, first marginally and then sharply after a given point.In essence, this  (Ryberg, 2018;Steffen, et al., 2015;Lenzen et al., 2013;Kanemoto et al., 2016;Kanemoto et al., 2014;Wiedmann et al., 2015;Lenzen, 2013;Oita et al., 2016) Climate change -Energy imbalance 1 Wm À2 (Ryberg, 2018;Steffen, et al., 2015;Lenzen et al., 2013;Kanemoto et al., 2016;Kanemoto et al., 2014;Wiedmann et al., 2015;Lenzen, 2013;Oita et al., 2016) Stratospheric ozone depletion 14.5 DU (Ryberg, 2018;Steffen, et al., 2015;Lenzen et al., 2013;Kanemoto et al., 2016;Kanemoto et al., 2014;Wiedmann et al., 2015;Lenzen, 2013;Oita et al., 2016) Ocean acidification 0.69 O Arag (Ryberg, 2018;Steffen, et al., 2015;Lenzen et al., 2013;Kanemoto et al., 2016;Kanemoto et al., 2014;Wiedmann et al., 2015;Lenzen, 2013;Oita et al., 2016) Biogeochemical flows -P 9.9 Tg PÁy -1 (Ryberg, 2018;Steffen, et al., 2015;Lenzen et al., 2013;Kanemoto et al., 2016;Kanemoto et al., 2014;Wiedmann et al., 2015;Lenzen, 2013;Oita et al., 2016) Biogeochemical flows -N 62 Tg N y -1 (Ryberg, 2018;Steffen, et al., 2015;Lenzen et al., 2013;Kanemoto et al., 2016;Kanemoto et al., 2014;Wiedmann et al., 2015;Lenzen, 2013;Oita et al., 2016) Land-system change -Global 25% (Ryberg, 2018;Steffen, et al., 2015;Lenzen et al., 2013;Kanemoto et al., 2016;Kanemoto et al., 2014;Wiedmann et al., 2015;Lenzen, 2013;Oita et al., 2016) Freshwater use -Global 4000 km 3 y -1 (Ryberg, 2018;Steffen, et al., 2015;Lenzen et al., 2013;Kanemoto et al., 2016;Kanemoto et al., 2014;Wiedmann et al., 2015;Lenzen, 2013;Oita et al., 2016) is accomplished by replacing SMR hydrogen with CCS with electrolytic hydrogen.The BAU solution, computed with the process ''1 kg Methanol {GLO}| market for | APOS, U" in Ecoinvent, is the cheapest alternative, yet it leads to the largest transgression level.The sharp increase in the slope of the curve found at an MTL value of 0.15 is due to the max operator used in the calculation of the environmental objective.Notably, starting from the minimum cost solution, the model reduces the transgression level by decreasing the impact on the climate change, ocean acidification and N flow PBs, as the other PBs are not transgressed.After a given point, the carbon-related PBs are all met (i.e., CO 2 concentration, energy imbalance and ocean acidification) and the only one that remains transgressed is the N flow.Hence, further reductions in the global transgression can only be attained by improving the performance in the N flow, which implies that the transgression level decreases more slowly (i.e., larger cost increases are required to attain a given change in transgression level).Fig. 7 shows the transgression level in each solution, where the green zone represents the SOS and the numbers on the axis denote the transgression level.The minimum transgression level solution substantially improves the environmental performance of the BAU (i.e., form three PBs transgressed in the BAU to only one in the minimum transgression solution).The BAU performs particularly poorly in the carbon-related PBs, i.e., CO 2 concentration, energy imbalance and ocean acidification, with transgression levels above 25 in the case of the climate change PBs and eight in the case of ocean acidification.The minimum cost solution reduces these transgression levels substantially, but not enough to operate within the SOS of these Earth systems.In contrast, the minimum transgression manages to reduce the impact below the assigned share of the SOS.Notably, the latter point attains negative transgression levels (on a cradle to gate basis) in the carbon-related PBs due to the use of CO 2 from DAC and the low carbon footprint of the electrolytic hydrogen (i.e., 2.27 kg CO 2 eq/kg H 2 , vs. 5.56 in the SMR with CCS).
No single design, including the BAU, trespasses the stratospheric ozone depletion, biogeochemical flow of phosphorus, land-system change and freshwater use PBs.The BAU does not transgress the N flow, while the other two solutions do.This PB, however, is more strongly connected to the agriculture sector due to the use of fertilizers (Sterner, 2019).Hence, this sector could offset this excess of impact in exchange for shares of the carbonrelated PBs, where the minimum transgression solution performs extremely well.
Focusing on the specific features of each design (Table 4 and further details in the supplementary material, section S4, Tables S23 -S28), we find that the extreme solutions display similar design features and operating conditions.The molar ratio of hydrogen -CO 2 is very close to the stoichiometric one (i.e., 2.80:1).The temperature in the reactor is around 215 °C, while the pressure is maintained at 49 bar.The reflux ratio and the volume of the reactor lie close to the lower bound.The optimal purge also approaches the lower bound, as otherwise the direct emissions and feed requirements would increase.The main difference lies in the source of the hydrogen feedstock, either from SMR with CCS or water electrolysis, depending on whether we seek to optimize the cost or the transgression level, respectively.
For a minimum temperature difference of 10 °C, the composite curve (Fig. 8), which is very similar in the extreme Pareto points, shows that only cooling utility (i.e., cooling water) is required.Fig. 9 displays the optimal HEN, which is the same in the extreme solutions.The heat of reaction, strongly exothermic, suffices to fulfill all the heating requirements of the system; consequently, only cooling water is required to meet the target temperatures.
Fig. 10 shows the main contributions to the environmental impacts.As seen, the hydrogen and CO 2 feeds are responsible for most of the impact in all of the cases.The contribution of CO 2 is negative in the carbon-related PBs, as this CO 2 comes from DAC and is, therefore, modeled as a negative flow from the environment.In the other PBs, however, it shows a positive contribution due to the impacts linked to the DAC facility.Electricity and direct emissions contribute mostly to the carbon-related PBs (i.e., climate change and ocean acidification).Heating has no impact since the designs do not require heat, thanks to heat integration.The impact of cooling is negligible.
In the two climate change planetary boundaries, CO 2 concentration and energy imbalance, and ocean acidification, the TL of the environmental optimum, which consumes only electrolytic hydro- gen, is negative.However, in the economic optimum, the TL is positive and this solution even transgresses the safe operating space, yet it is still lower than in the BAU process.Both extreme solutions consume almost the same amount of hydrogen (Table 4), and the differences in impacts are due to the hydrogen source.
In the stratospheric ozone depletion, biogeochemical flows of phosphorus and nitrogen, and fresh water use PBs, both optimal designs perform worse than the BAU process.In the case of the biogeochemical flow of nitrogen, the designs even operate above the share of the safe operating space.In the biogeochemical flow of phosphorus, it can be seen that the impact of the electrolytic hydrogen is much higher than that of its SMR with CCS counterpart.In the fresh water use PB, most of the impact comes from the CO 2 from DAC.
In the land-system change PB, the environmental optimum outperforms the BAU, while the economic optimum is worse than the BAU.This difference between the two optimal designs, as before, comes from the hydrogen source.
We next further highlight the advantages of using absolute sustainability metrics in process design by assessing the LCA impact of the environmental and economic optimal solutions using the ReCiPe 2016 methodology following an Egalitarian approach (results in Table 5).
As seen, the minimum impact solution performs worse in two of the endpoints, yet the PBs analysis clearly shows that it leads to the minimum transgression level.We clarify that the PBs omit human health impacts; moreover, they consider the impact on ecosystems via the biodiversity intactness index, yet this metric was omitted in our work due to methodological gaps.On the other hand, most of the impacts on the control variables of the PBs are quantified via midpoint indicators, which are summed up to yield the endpoint scores.For example, there are midpoints for freshwa-   ter use, ozone depletion, climate change, eutrophication (linked to the biochemical flows of N and P) and land use (related to land system change), and others for toxicity, particulate matter, ionizing radiation and resources, omitted in the PBs.
Hence, LCA indicators and the PBs metrics should complement each other because the standalone use of the former could lead to solutions that are unsustainable in absolute terms.Furthermore, the ReCiPe endpoints are hard to interpret in absolute terms, so at most, they could enable comparisons between technologies but fail to provide insight into their worldwide sustainability level.

Conclusions
This work introduced an approach to design chemical processes that are entirely sustainable considering the ecological capacity of the planet.Our method combines process simulation, machine learning algorithms and an absolute sustainability assessment model based on the PBs concept.The methodology was applied to the production of green methanol from hydrogen and carbon dioxide taking the BAU process as a baseline.
Green methanol from electrolytic hydrogen powered by wind cannot ensure a safe operation within all the PBs due to the transgression (according to the egalitarian sharing principle applied) of the biogeochemical flow of nitrogen.However, it can greatly reduce the global transgression level compared to the BAU and to the CO 2 hydrogenation route that uses hydrogen from SMR with CCS.Methanol from hydrogen and carbon dioxide is at least 1.62 times more costly than the BAU process (using hydrogen from SMR with CCS), and much more expensive when using electrolytic hydrogen.The optimal process conditions and sizes are quite similar in the extreme Pareto points, and only differ in the source of hydrogen, which dictates the economic and environmental performance.
Overall, our results unfold new avenues for the incorporation of absolute sustainability criteria in process design, with emphasis on the evaluation of emerging decarbonisation routes.Designing processes that operate within global ecological budgets will become more critical as we approach the biophysical limits of the Earth.In this context, our tool provides a starting point to design processes that are entirely consistent with the ecological capacity of the planet.Our method should not be seen as an alternative to   existing methods, including LCO and LCA approaches, but rather as a complement to them that can provide further insight into the sustainability level of alternative processes.Future work should focus on handling the main uncertainties affecting the calculations and enlarging the scope of the analysis to embrace other economic sectors linked to the chemical industry.

ab
Share of the total safe operating space allocated to the studied process IMP DContribution to the transgression of planetary boundary b per functional unit SOS bTotal safe operating space for planetary boundary b absolute terms, as they refer to the Earth's ecological capacity.

Fig. 2 .
Fig.2.Scope of the production system investigated in the case study.

Fig. 3 .
Fig. 3. Diagram of the flowsheet highlighting in colors the input variables in the surrogate model.

Fig. 4 .
Fig. 4. Neural network fitting performance results, as well as chosen parameters.

Fig. 5 .
Fig. 5. Comparison between the real (process simulation) and surrogate values (neural network).R 2 , MSE and MAE are shown for each output.

Fig. 6 .
Fig. 6. Results of the multi-objective optimization for ten iterations of the e-constraint method.

Fig 7 .
Fig 7.Results of the TL in the PBs using egalitarian downscaling for the BAU, the optimal environmental design and the optimal economic design.

Fig. 8 .
Fig. 8. Hot and cold composite curves for the optimal environmental solution.

Fig. 9 .
Fig. 9. SYNHEAT solution for the optimal HEN.H refers to hot streams and C to cold streams.

Fig. 10 .
Fig.10.Disaggregation of the transgression levels per source of environmental impact for the environmental and economic optimal design, as well as the total value of the transgression level of these designs and the BAU process.

Table 1
Input variables to the surrogate model and associated ranges.

Table 3
Value of the SOS for each planetary boundary.

Table 4
Results of the process simulation for the extreme Pareto points (minimum transgression and minimum cost).

Table 5
Resultsof ReCiPe 2016 Endpoint (E) for the optimal solutions obtained.