Incorporating thermoelectric power plant water use into multi-objective optimal power flow

Traditionally, power systems have been operated to minimize cost while maintaining reliability. However, extreme weather and demand events can affect traditional thermoelectric power generation operations due to their reliance on water for cooling. This paper contributes a novel multi-objective formulation of the optimal power flow (OPF) problem where cost, water withdrawal, and water consumption are minimized. Through this formulation, we assign optimization weights to water withdrawn and consumed, which can be directly incorporated into existing OPF formulations. We apply this formulation with a global mapping sensitivity analysis to a realistic case study to first demonstrate its general effectiveness under extreme climatic, hydrologic, and operational scenarios. Then, we apply a global ranking sensitivity analysis to determine the most influential generators for system performance. Through this operational scenario analysis framework, analysts can gain insights into potential system-level and component-level vulnerabilities within power systems. Such insights can be useful for informing both short-term operations as well as long-term power system planning.


Introduction
Large-scale thermoelectric power generation in the United States requires significant volumes of water to be diverted-from or impounded-in river systems. Thermoelectric water use can be categorized into water withdrawn (i.e., water directly removed from the water system and might or might not be returned) and water consumed (i.e., water that does not get returned to the water system). Currently an estimated 48 percent of fresh surface-water withdrawals can be attributed to all the uses of thermoelectric power plant electricity generation (Dieter et al 2018). Often, thermoelectric water withdrawal and consumption are quantified by a water withdrawal and consumption rate (i.e., the amount of water used to produce a given amount of electricity). These rates change over time as they are impacted by extreme climatic events such as droughts, heat waves, or upstream water users (Scanlon et al 2013a, Lydersen 2016, Meng and Sanders 2019.
Electric loads on the grid also change over time. On shorter operational timescales, electric loads vary throughout the day based on user types within an electric system (Merrick 2016, EIA 2020a. These demands also vary on longer timescales such as seasonal changes (EIA 2020a, Archibald et al 1982) as well as new loads (e.g., electric vehicle charging or new development) being added to the grid (Hatziargyriou et al 2013.
The temporal variability of electric loads, thermoelectric withdrawal rates, and thermoelectric consumption rates are particularly important during extreme natural events. These events can be broadly categorized into events that impact electric loads, events that impact thermoelectric water use rates, and events that impact both. For example, a drought event may impact the thermoelectric water use rates while a heat wave may impact the electric loading (Scanlon et al 2013a, Harto et al 2012. Such system-stressing events can also occur simultaneously, further stressing energy systems (Galbraith 2011, Freedman 2012. In these extreme operational scenarios, system operators often manually intervene in traditional operations (e.g., shutting down specific generators or granting thermal variances) (Harto et al 2012, Poumadère et al 2005, Lubega and Stillwell 2018. These extreme operational events have also prompted recent work for how to incorporate thermoelectric power plant water use into long-term power systems planning. Previous studies have proposed water-informed multi-objective formulations for decadal-scale power system expansion planning problems (Jornada and Leon 2016). Other studies have highlighted the importance of how upgrading and retrofitting generator cooling systems can impact long term system performance. On a smaller timescale, seasonal water-informed frameworks have also been proposed for granting thermal variances to generators (Lubega and Stillwell 2018) or planning for a drought event (Mu et al 2020). However, fewer studies incorporate water use into the daily and sub-daily operations of the power system. Such water-informed operations are important as they (1) promote system coordination rather than manual intervention by system operators, as well as (2) offer a rapid way for system operators to respond to current drought events, especially if long-term solutions are infeasible due to a lack of time to implement (Pacsi et al 2013).
We propose a novel operational scenario analysis framework that builds on previous works to incorporate water into grid operations. Specifically, this framework utilizes a water-informed formulation of the optimal power flow (OPF) problem. The OPF problem optimizes active power output at generators such that demands are met and the system operates within its designed limits. Typically this problem is solved by system operators on sub-hourly timesteps.
Previous work has proposed water-informed OPF formulations that minimize the sum of the monetary cost of energy production as well as the cost of water consumed for cooling needs, while constraining solutions based on the amount of water to be consumed by power plants and the heat constraints on the thermal discharge to the environment (Fooladivanda and Taylor 2015). These studies show that solutions to a waterinformed OPF, as well as overall system objective performance, is impacted by upstream thermoelectric power plant water use as well as system-wide hydrologic conditions such as streamflow and stream temperature. The novel water-informed OPF formulation extends previous work by simplifying the formulation as to keep the formulation convex. Such a convex formulation reduces the computational burden when solving the problem, which is essential for real-time operations. Our formulation achieves convexity by modeling water withdrawal and consumption in the objective function as has been done for similar power system problems (Sanders et al 2014). Another advantage of modeling water withdrawal and consumption in the objective function, as our formulation does, is that determining hard constraint values on water consumption and withdrawal are difficult to assign when looking at operations over longer time horizons where water availability and rights may be uncertain.
Solving a single case of a water-informed OPF allows for insights into a single snapshot of system performance. For example in their demonstration of their formulation, Fooladivanda and Taylor (2015) assume a single price of water consumed ($1/m 3 ), load demand case, withdrawal rates, and consumption rates. However, Lubega and Stillwell (2019) found that a thermoelectric cooling water price is a function of many dynamic economic factors. Furthermore, a $1/m 3 cost is a high price for cooling water when compared to traditional market prices Stillwell 2019, Brewer et al 2007). Similarly, loads, withdrawal rates, and consumption rates are all impacted by extreme demand and climate events as previously stated. Many studies have modeled loads as a stochastic process to attempt to capture the uncertainty of the temporal variability of loads in the grid (Yong and Lasseter 2000, Hu et al 2010. Such studies have found that by including this variability, solutions to the OPF problem are drastically different from their deterministic counterparts. Therefore, to allow insights into longer-term system performance, our operational scenario analysis framework considers the water optimization weights, generator water consumption rates, generator water withdrawal rates, and system loadings as 'exogenous parameters' that are outside of the control of the system operator, and therefore could be different every time the OPF is solved, due to system conditions. Such a novel consideration allows us to isolate the impacts of various inputs to the water-energy system. This framework also extends the previous literature by conducting several sensitivity analyses to quantify the individual significance and effects of each exogenous parameter. We show that such a novel framework allows analysts both system-wide insights into the operations of a power system as well as identifying the most critical generators that impact various aspects of system performance. This framework allows the most important generators of a power system to emerge from the underlying mathematics of short-term power systems operations, rather than assuming generator importance a priori. We first define our operational scenario analysis framework consisting of our water-informed OPF formulation, an embedding of the water-informed OPF into a larger multi-objective model, and two sensitivity analyses of the multi-objective model. We then Figure 1. Multi-objective water-informed OPF model. Exogenous parameters, those parameters out of the control of system operators, are inputs to the multi-objective model. These parameters are used in the water-informed OPF formulation to find the optimal generator output. These parameters are also used to compute objectives as indicated. Every feasible set of exogenous parameters corresponds to a set of objective values. demonstrate the insights possible through the operational scenario analysis framework on a synthetic case study located in the Midwestern United States.

Methodology
An overview of the multi-objective water-informed OPF model is presented in figure 1. First, we define all the exogenous parameter inputs to the water-informed OPF formulation. Then, the water-informed OPF finds the optimal power output at each generator. Finally, multiple objectives are a function of the optimal power outputs and the exogenous parameters.

Water-informed OPF
The water-informed OPF problem (P opf ) finds the set of generator power outputs p g that minimizes a weighted objective function while also satisfying the DC approximation of the power flow equations. The inputs of (P opf ) are exogenous parameters, while the outputs of (P opf ) are the 'state variables' of generator power outputs p g . As previously mentioned, exogenous parameters are outside of the control of the system operator, and therefore could be different every time the OPF is solved, due to system conditions. In this formulation, these exogenous parameters include the water withdrawal and consumption rates at each generator, the loads at each bus, as well as the weights of water withdrawn and consumed.
The withdrawal rate and consumption rate at each generator represent a generator's 'water efficiency', the amount of water used to produce a given amount of electricity. The withdrawal rate at generator j is denoted as β with,j . Similarly, the consumption rate at generator j is denoted as β con,j . These terms are collected into vectors β with and β con of length G, the number of generators in the system.
Value ranges of withdrawal and consumption rates can be estimated based on historic datasets of electricity production, water withdrawal, and water consumption. Additionally, general ranges of withdrawal and consumption rates have been estimated based on fuel and cooling system types (Macknick et al 2012, Rutberg et al 2011. Physically, withdrawal and consumption rates can be interpreted such that water-efficient generators would have lower rates than water inefficient generators (i.e., it takes water efficient generators less water to produce a unit amount of electricity).
Other exogenous parameters include the loads p l,i at each bus i. These terms can also be collected into a vector p l . The set of buses in the system N has the length N. Therefore, p l also has the length N (buses with no load have a zero as its corresponding entry in p l ).
The two remaining exogenous parameters in this formulation are the optimization weight associated with the water withdrawn w with and the optimization weight associated with water consumed w con . Each of these terms specifies a monetary penalty assigned to generators that withdraw or consume water. The units of these weights are a monetary cost per volume of water used. For example, a solution that assigns a low value of the withdrawal weight w with and a high value of the consumption weight w con gives preference to generators that do not consume water. Therefore, generators that consume a lot of water would have a low power output while generators that consume less water will have a higher output. Through varying values of w with and w con , we create operational scenarios based on differing levels of withdrawal and consumption preference. These preferences could be correlated with constraints in the water system. In real-world systems these weight values could be in the form of a monetary tax from a policy focused on water use, for example. A summary of all the exogenous parameters is presented in table 1. Our water informed OPF formulation (P opf ) is defined as follows: P opf (w with , w con , β with , β con , p l ) = p * g = arg min p g j∈G a j + b j p g,j + c j p 2 g,j + w with j∈G β with,j p g,j t + w con j∈G β con,j p g,j t (1) where p g is a vector of active power generation p g,j for each generator j in the set of generators G. B im is the (i, m) entry of the susceptance matrix B. In this formulation, B captures both the topology of the transmission lines in the grid as well as their line parameters, which dictate how power flows through the network. The difference in the voltage phase angles between neighboring buses i and m is written as θ im . The product of B im and θ im approximates the electricity flowing through line im. Therefore, equation (2) is an approximation of the power flow equations that relates the differences in loads and generation to the flow of electricity moving through the lines for all the buses. Equation (3) ensures that these flows are less than the line limits F max,im for all lines im in the set of lines L. The upper and lower active power capacities of each generator are collected in the vectorsp g and p g , respectively. Equation (4) limits the active power for all generators to be within these limits. For further explanation of the full power flow equations and justification for the DC approximations see (Taylor 2015, chap 3). The objective function is a weighted sum of generator costs, water withdrawal, and water consumption. The coefficients of a j , b j , and c j reflect the monetary cost as a function of active power output for each generator j. β with,j and β con,j reflect the water withdrawn and consumed as a function of active power output for each generator j. t represents the timestep conversion factor of the optimization problem, which converts the power output of generators to an energy output such that the water withdrawal and consumption rates can be applied (see supplemental information (https://stacks.iop.org/ERIS/2/015005/mmedia) for details on units). As previously discussed, w with and w con capture the monetary cost of the water withdrawn and consumed, respectively. The units of these quantities are derived in the supplemental information.

Multi-objective water-informed OPF model
We propose several objectives that each capture unique and analyst-relevant aspects of a system's performance. One such measure of the system performance is the objective function from the water-informed OPF. This objective captures the 'total effects' of the system (i.e., it captures the generator costs, costs from withdrawal, and costs from consumption). The total cost objective is defined as The total cost objective F cos is critical as it captures both the direct monetary costs of operations as well as the costs associated the generator water usage. As was visually depicted in figure 1, it is the only objective that directly uses all the exogenous parameters. Another measure of system performance is the summation of all the direct generator operational costs defined as This 'traditional' OPF objective of generator costs serves as a reference for analysts familiar with system performance under traditional OPF formulations. As is visualized in figure 1, this objective is the only objective that is solely a function of generator power outputs, such that it also serves as a proxy for measuring total system generator outputs and will not be directly impacted by varying exogenous parameters. It will only be impacted by exogenous parameter values that change the generator power output values (p g ). The final two objectives capture the total amount of water withdrawn and consumed in the system and are defined as Recall that t is a conversion factor to convert the power output to an energy output over a timestep. These objectives serve as a ground truth of the volume of water that would be used in the system over the optimization timestep (and does not consider the withdrawal weight w with or consumption weight w con ). Such information is important for operations during climatic extremes where the total amount of water usage can have detrimental ecological effects as well as operational impacts. A summary of the four objectives is also presented in table 1. The previously defined water-informed OPF formulation, as well as these four objectives, are incorporated into a multi-objective model (P moopf ) defined as P moopf w with , w con , β with , β con , p l = F cos p * g , w with , w con , F gen p * g , F with p * g , β with , F con p * g , β con (9) s.t. p g * = P opf w with , w con , β with , β con , p l .

Sensitivity analysis of uniform variability
Our operational framework focuses on two ways that exogenous parameters can vary: uniformly and nonuniformly. We define uniform variability as changes that impact exogenous parameters in all parts of the system in the same way. For example, during some extreme loading event, a system may experience increased loading across all buses. However, non-uniform variability impacts each exogenous parameter differently. For example, the water use of an upstream utility may only impact generator withdrawal rates at a single downstream plant. We introduce different coefficients on each exogenous parameter to mathematically model the different types of variability. Our operational scenario analysis framework utilizes two different types of sensitivity analyses of these introduced coefficients to gain insights into both the uniform and non-uniform variability of our multi-objective water-informed OPF model. The demands at each bus (p l ) are assumed to vary uniformly throughout system operation (i.e., all the loads increase at once). Mathematically, we present this uniform variability as a uniform loading coefficient (c load ) multiplied by p l . A uniform loading coefficient (c load ) greater than one means that the system is experiencing increased loading relative to its baseline loading case. Such uniform increases would generally be the result of higher system-wide temperatures being linked to higher system-wide loading (PJM n.d.). The generator withdrawal rates (β with ) and generator consumption rates (β con ) also vary uniformly. Similar to the loading uniform variability, the uniform water exogenous parameter variability is mathematically represented as a uniform water coefficient c water multiplied by β with and β con . A uniform water coefficient (c water ) less than one means that the generators are becoming more water efficient (using less water to produce the same unit of electricity) while a coefficient greater than one means the generators are becoming less water efficient (using more water to produce the same unit of electricity). Such variability has been categorized across cooling technology, fuel technology, space, and time (Macknick et al 2012, Peer and Sanders 2016, Huston 1975. We focus on sources of variability in time, based on assuming a system with non-changing cooling technologies and generators. Uniform variability, i.e., processes that would impact all generators in a given system, include seasonal, climatic, and hydrologic/hydraulic variability (Lee et al 2020, Clement et al 2017, Dziegielewski et al 2002. In the uniform sensitivity analysis, we test the sensitivity of the exogenous parameters w with and w con , as well as the uniform loading coefficient (c load ) and uniform water coefficient (c water ). These four parameters of our sensitivity analysis are referred to as 'input factors' to avoid confusion with the entire set of exogenous parameters (Pianosi et al 2016). The goal of this analysis is to (1) identify how the water and load uniform variability interact with one another to gain insights into what general mechanisms impact system performance, as well as (2) determine which regions of the input factor space cause extreme objective performance for each of the objectives in our multi-objective water-informed OPF model. A summary of the four input factors for the uniform sensitivity analysis is presented in table 1.
We use a tree-based global mapping sensitivity analysis to analyze the uniform variability (Pianosi et al 2016, Singh et al 2014, Harper et al 2011. First, we conduct a sampling of the input factor space as depicted in figure 2. This sampling creates a set of input factor scenarios and each scenario's corresponding four objective values. The underlying assumptions of this sampling are dependent on the analysts' access to data as well as system configurations. After the sampling, we fit a unique classification and regression tree (CART) model to predict each of the four objectives as a function of the input factors over many input factor scenarios. In section 4.1, we demonstrate how to find the bounds of these input factors, one possible way to conduct sampling of the input factor space, as well as how to interpret the resulting decision trees.

Sensitivity analysis of non-uniform variability
As previously discussed, the non-uniform variability of generator withdrawal rates (β with ) and generator consumption rates (β con ) can be attributed to plant-to-plant discharging interactions, a change in cooling system operational policies, or other water discharges in the system (Tidwell et al 2019, Miara et al 2018, Yang and Dziegielewski 2007, Fooladivanda and Taylor 2015. Physically, these processes all impact different fuel/cooling type combinations differently and pose a key operational challenge. The non-uniform sensitivity analysis allows us insights into how these non-uniform processes impact system-wide performance under several operational scenarios. These operational scenarios are defined by the input factors already examined in the uniform sensitivity analysis: the withdrawal weight w with , the consumption weight w con , and the uniform loading coefficient c load . Unlike the uniform sensitivity analysis where a single uniform water coefficient (c water ) impacted generator water use at all plants, in the non-uniform case, we introduce unique non-uniform water coefficients for each plant p denoted as c water,p . Assuming there are P unique plants in a system, there will be P non-uniform water coefficients. This formulation assumes that generators in the same plant have the same fuel and cooling types such that no processes impact single generators differently within a plant. If a single plant location consists of multiple generator/cooling system types, a unique non-uniform water coefficient is assigned to each set of unique generators/cooling system types. This condition is demonstrated in the case study in section 4.2.
The P non-uniform water coefficients constitute the input factors of the non-uniform sensitivity analysis. Just as in the uniform sensitivity analysis, the input factor space must be sampled for each operational scenario. The sampling of the non-uniform input factor space is visualized in figure 3. Starting with the outermost layer, we see the operational scenario defines the values for the withdrawal weight w with , the consumption weight w con , and the uniform loading coefficient c load . Then, the input factors, whose sensitivity is tested, are the nonuniform water coefficient of each plant in the system from p to P. Again, just as is in the uniform case, each set of input factor scenarios produces a set of objective function values. This entire sampling is repeated for each operational scenario.
A variance-based sensitivity analysis is implemented for the non-uniform variability analysis because it allows us to rank and quantify significance of input factor relationships. Variance methods assume that the input factors are stochastic variables that produce distributions of each objective value. Therefore, they assume that the variance of the output distribution is a suitable proxy for uncertainty and that the contribution of an input factor to the output variance is a measure of its uncertainty (Pianosi et al 2016). Several indices have been proposed to quantify this process, with the Sobol index being widely accepted. The first order Sobol index (S First i,j ) for objective j captures the 'main effect' of an input factor x i and is defined as where y is the objective value for objective j, V is variance, E is the expected value, and x ∼i is all the input factors except the ith value. Conceptually, S First i,j is the expected reduction in output variance that can be obtained when fixing x i (Pianosi et al 2016, Sobol' 1990. Physically, S First i,j quantifies how important input factor x i is to objective j. For example, S First c water,1 ,F cos quantifies how important the water use at plant 1 is to the total cost of the system.
In section 4.2, we demonstrate a nonparametric method for both sampling the non-uniform input factor space as well as estimating the first order Sobol index. However, many methods for estimating Sobol indices could be implemented for a non-uniform sensitivity analysis as long as all assumptions are satisfied.

Case study
The following section describes the synthetic yet realistic case study chosen to demonstrate our operational scenario analysis framework. We show how we use EIA-reported thermoelectric water use information to ensure realistic water use for our synthetic test system. This water-use-assigned test case is then used to demonstrate the insights possible from both the uniform and non-uniform sensitivity analyses. An overview of the case study is provided in figure 4.

Synthetic grid: Illinois 200-bus system
Due to national security concerns, models of the actual electricity grid contain sensitive data, so many frameworks have been developed to generate synthetic grids that retain many valuable properties of the true grid while not revealing sensitive data (Birchfield et al 2017). The framework generally consists of two phases: substation placement and transmission topology specification. The substation placement phase specifies the locations and specifications of loads and generators. This phase is based on a weighted distance clustering algorithm of loads in each postal code and known generator locations. One constraint in the algorithm is that hydroelectric, nuclear, or renewable energy resources cannot be grouped with nonrenewable energy sources. The final results of this phase are synthetic but realistic substation locations containing loads and generators. The second phase then connects these substations based on an iterative algorithm with several criteria based on line property constraints as well as line topology constraints. Finally, a variety of other devices are inserted in the network to regulate reactive power and voltage, such as shunt capacitor banks, transformer taps, and synchronous condensers.
We use the Illinois 200-bus synthetic grid currently hosted in the Electric Grid Test Case Repository (Birchfield et al 2017). Instead of simulating economic dispatch, we manually dispatch all generators in the system. However, simulation of other operational processes is discussed in section 5. Recall, that these generators do not explicitly correspond to true generators but are rather synthetic generators based on an aggregation. Fortunately, these synthetic generators do have coordinates on a coordinate reference system specified by the distance weighting algorithm and are available in the PowerWorld formats of the synthetic case. Figure SI1 shows these locations in the supplemental information. The locations of the generators reported to the United States Energy Information Agency (EIA) are also depicted. These generators are called 'EIA generators' for the remainder of this analysis.
The Illinois 200-bus synthetic grid does not initially include withdrawal rates (β with ) or consumption rates (β con ) for the synthetic power generators. Therefore, we infer realistic withdrawal rates (β with ) and consumption rates (β con ) in a two stage process. As explained below, we infer the synthetic generator cooling system type based on location and characteristics. Then, we determine withdrawal rates (β with ) and consumption rates (β con ) from regional aggregations of publicly available power systems data (EIA 2020b).

Assigning synthetic generator cooling systems
We first assigned each synthetic generator a cooling system. We manually matched the synthetic generators to EIA generators who reported data on Forms 860 and 923. Figure SI1 in the supplemental information shows both synthetic and EIA generators. From inspection of this figure, we assign cooling systems to the synthetic generators. For each synthetic generator, we assign a cooling system based on an EIA generator that is spatially close, has a matching fuel type, and has a similar capacity to the synthetic generator. This process ensures that the synthetic generators are assigned representative and realistic cooling systems.
3.3. Regional estimates of generator water use We gathered regional estimates of withdrawal rates (β with ) and consumption rates (β con ) for the fuel type/cooling system combinations present in our case study. Specifically, we use data from EIA-reported generators in the state of Illinois, as this approach ensures that water use data come from generators that have the same political regulations and climatic conditions as the synthetic generators. Because the EIA generator dataset is generated from human-reported data, this dataset is notoriously noisy and often requires preprocessing (Tidwell et al 2019, Meng and Sanders 2019, Grubert and Sanders 2018, Harris and Diehl 2017, Peer and Sanders 2016, Scanlon et al 2013b, Averyt et al 2013. We build upon the frameworks presented in these studies and implement the following filtering criteria for the regional EIA dataset: • Generators must be located in Illinois. • Only consider generator fuel type/cooling system combinations present in synthetic generators • Only keep observations that reported finite withdrawal and consumption rates (i.e., remove zeros or division by zero errors) • Only keep generators that reported water use at least half of the months between the period of observation (2014-2019) • Remove outlier withdrawal and consumption values based on the modified Z-score test (Iglewicz and Hoaglin 1993) for each generator, which was previously used to filter the EIA dataset in Peer and Sanders (2016) The distributions of the resulting regional EIA-reported estimates of withdrawal rates (β with ) and consumption rates (β con ) are presented in figure 5. This plot is consistent with the findings from Peer and Sanders (2016), where we see that once-through cooling systems have a more pronounced impact on withdrawal rate than fuel type.
With regional estimates of withdrawal rates (β with ) and consumption rates (β con ) for the generator fuel type/cooling system combinations in our synthetic grid, we find representative values of water use for each synthetic generator.

Assigning synthetic generator water use via K-means clustering
We use the regional estimates of withdrawal rates (β with ) and consumption rates (β con ) to assign the withdrawal rates and consumption rates for each synthetic generator. For the nuclear and natural gas regional generators, the median value is justified due to the symmetry of the distributions and lack of variation. However, the coal generators exhibit notable spread in withdrawal and consumption rates, as shown in figure 5. The variation of coal generator withdrawal rates (β with ) and consumption rates (β con ) can be attributed to the variation of coal generator capacities for this region. This trend emerges by plotting the withdrawal rates (β with ) and consumption rates (β con ) against the capacities for the coal plants, shown in figure SI2 in the supplemental information. From visual inspection, several clusters emerge for each cooling type based on the capacities of the generators. The median value of each corresponding cluster informs the synthetic generator water use.
This regional-estimate-with-clustering methodology is not used for wind and solar generators, which do not require cooling systems, and have zero values for their exogenous water parameters. Similarly, small capacity natural gas combustion turbines are also assumed to use negligible amounts of water, so zero values are assumed for their exogenous water parameters (Stillwell et al 2011). The model of the synthetic grid with all thermoelectric water use information is available for public download (Kravits et al 2021).

Results
The following section demonstrates our operational scenario analysis outlined in section 2 for the case study. We perform a global sensitivity analysis of the uniform and non-uniform variability of our multi-objective water-informed OPF model. This section provides an example of both sensitivity analyses based on publicly available data. However, applying these sensitivity analyses to different datasets or system configurations may require modification of the specific implementation demonstrated in this analysis.

Impacts of uniform variability on system performance
We apply the uniform sensitivity analysis outlined in section 2.3 to our case study to first demonstrate the efficacy of the water-informed OPF formulation as well as examine how system-wide sources of variability impact system performance. Table 2 outlines the bounds of the four input factors for our case study. Starting with the loading in the synthetic grid, the algorithm in Birchfield et al (2017) creates realistic load profiles based on population densities. Additionally, we assume that the default loading supplied by this test case is generally reflective of normal loading conditions (i.e., c load = 1.0). We tested feasible uniform loading coefficients c load by increasing values until the traditional DC OPF became infeasible due to thermal line limits or lack of generation capacity. This simulated increase resulted in a maximum feasible c load value of 1.6. To ensure that this bound is reasonable given actual conditions, we also examined historic loads. The synthetic grid for this case study would be located in the Midcontinent Independent System Operator (MISO). MISO only provides load data for its entire system and not for individual hubs. The distribution of historic load data as well as the uniform  figure SI3 in the supplemental information. From this figure, we see that 1.6 is just out of the historic range of load coefficients so we set the max uniform load coefficient value to be 1.5 to ensure a value that is both feasible and has occurred in the historic record.
To ensure that this value is realistic, we observe that a value of 1.5 is roughly the 99th percentile of historic loading. Studies have shown that the frequency of the historic 99th percentile loading occurring is expected to increase up to 1500 percent depending on the region and climatic predictions utilized (Auffhammer et al 2017). Therefore, we set 1.5 as the upper bound of the uniform loading coefficient as it is both a historically realistic value, feasible for our system, as well as a potential planning metric for analyzing a power system's preparedness for uncertain future conditions. A similar process was followed to assign the ranges of the uniform water coefficient c water . First, well-behaved ranges were found to be between 0.5 and 1.5 (i.e., beyond these values the water-informed OPF becomes illconditioned and fails to converge). The distribution of historical c water values (figure SI4) suggested that our variable is in the range of historical bounds. Additionally, we use empirical estimates from standard thermoelectric water use models to show that such changes can result from reasonable changes to inlet cooling temperature water as is done in the supplemental information section SI6.
The water withdrawal and consumption optimization weights (w with and w con ) ranges are selected to prevent against the previously mentioned ill-conditions for convergence.
To grid sample the input factor space, each of the four input factors is divided into an equal number of steps between its bounds as outlined in table 2. This approach produces a uniform sampling of the input factor space yielding 10 000 combinations of the input factors (w with , w con , c water , and c load ), and the multi-objective waterinformed OPF model was run for each combination producing four objective values (F cos , F gen , F with , F con ). The sampling was run on the RMACC Summit supercomputer across one node with two CPUs (Intel Xeon E5-2680 v3, 2.50 GHz, 24 cores/node) with 4.84 GB RAM per core, taking a wall time of roughly 12 min.
A subset of this sampling is depicted in figure 6. This figure shows both how our water-informed OPF reduces system-wide water withdrawals as well as how a water-informed system is able to take on additional loads with marginal increases to overall water use. This subset is chosen to demonstrate extreme values for the withdrawal and consumption coefficients. The horizontal axis of figure 6 represents increasing electricity demand. The vertical axis depicts the total water withdrawal objective F with , which was defined in equation (7). The color of the lines reflects the uniform water coefficient (c water ). Darker lines correspond to more severe system-wide water constraints (e.g., drought). Finally the marker styles reflect a subset of the minimum and maximum withdrawal weight (w with ).
Comparing system withdrawals between the same water constraint/increased electricity demand scenarios but different withdrawal weights in figure 6, we gain insights into the efficacy of the water-informed OPF. Visually, the two lightest lines (uniform water coefficient equal to 0.61) correspond to the highest water efficiency. As the system experiences increases in electricity demand (i.e., moves to the right along the horizontal axis), the withdrawal weight has an increasing effect on reducing overall system withdrawals. In fact, for loading coefficients less than 1.2, by considering a withdrawal weight, the system adapts to rely on plants that either have low or zero withdrawal weights to cover added loads. This outcome results in the system having only marginal increases to water withdrawal as the system experiences larger loads. Across varying system-wide water efficiencies (line color in figure 6), we see that the effect of the withdrawal weight increases as the system gets more water inefficient. To demonstrate the largest system-wide water inefficiencies (e.g., drought; uniform water coefficient equal to 1.5), we see that by assigning a withdrawal weight (e.g., implementing a policy focused on withdrawal), we reduce system withdrawal to a level expected at a better system-wide water efficiency (uniform water coefficient of 1.06).
Another trend from figure 6 is that linear increases in electricity demands do not cause linear increases to water use. This result might be expected given that water withdrawal is modeled as a linear coefficient on generator power output. The reasons for this nonlinearity arise from other power flow considerations in the water-informed OPF that cause individual generators to generate electricity with nonlinear increases even under linear increases to total system loadings. Figure 6. Effect of withdrawal weight on total system withdrawal. This figure shows both how our water-informed OPF reduces system-wide water withdrawals as well as how a water-informed system is able to take on additional loads with marginal increases to overall water use. This figure subsets the data such that the consumption coefficient w con = 0. Figure 7. Effect of withdrawal weight and uniform water coefficient on plant output. Samples subset such that consumption weight (w con ) equals 0.0. Subsets of the uniform water coefficient (c water ) and the withdrawal weight (w with ) depicted to right. Abbreviations: recirculating with induced draft cooling tower (RI), recirculating with cooling ponds (RC), once-through with cooling ponds (OC), no cooling system (NCS), natural gas (NG).
To gain insights into what is happening on the plant-level, we depict plant output capacity ratios in figure 7, with each plant shown as a different line color and type. Generator output capacity ratios are defined as the ratio of the current power output to its total power output capacity, shown on the vertical axis of figure 7. As in figure 6, the electricity demand (uniform loading coefficient) is plotted on the horizontal axis. Each panel moving from left to right reflects a different group of plants with the same fuel/cooling system type (with unique legends for each row on the bottom). Each row of the figure reflects a different subsetting of the data as denoted by the labels on the right. From figure 7, we see the general trend that as loading increases (moves along the horizontal axis), plants either output the same or more power (vertical axis). For example, there is no case when system loading increases but an individual generator outputs less. This behavior is expected given the nature of power systems and their governing equations. However, a more interesting insight from figure 7 is which plants output more power under linear increases to system load. Looking at the top row of figure 7 (lowest water constraint with traditional OPF), we see that the increased electricity demand is generally supplied by the coal and nuclear plants. We see initial electricity demand increases are first supplied by the nuclear generator, which quickly reaches maximum generation as uniform loading coefficients increase beyond 1.2. Continued electricity demand increases are then supplied by the coal plants with once-through and recirculating induced draft cooling systems. Within the coal/once-through plants, the cheaper plants reach maximum capacity quickly (such as plants 9 and 12), whereas the more expensive plants reach capacity much slower (plant 6). Certain plants, like 1 and 6, never reach maximum generation, even at the highest uniform loading case. This finding suggests that something else about this network, like a line constraint, is leading to infeasibility of solutions beyond uniform loading coefficients of 1.6 as there is extra generation capacity.
Comparing across rows in figure 7, we see the impact that the withdrawal weight has on individual plant power output. Recall that by putting a high weight on the withdrawn water, the system operates to minimize power output of plants with high withdrawal rates. This condition is clear when comparing the coal/oncethrough plant output. In the top row (traditional OPF), the system relies on coal/once-through plants to cover the initial increases in electricity demand. However, when these water-inefficient plants are penalized, they do not output power until the most severe loading conditions (first column, second row of figure 7). In fact, the system is able to rely on plants with lower withdrawal rates (coal/recirculating induced) or plants with zero withdrawal rates (natural gas/no cooling system) to cover the initial load increases. Beyond loading coefficients of 1.2 the nuclear plant is selected, which causes the increase in water withdrawal. Specifically, this result is why system withdrawals did not increase until a loading coefficient of 1.2 as was observed in figure 6.
We can further subset the sampling data to observe the effect that applying a withdrawal weight has on the actual flows through the system. We choose two solutions that both have a consumption weight of 0.0, a uniform loading coefficient value of 1.5, and a uniform water coefficient of 0.5. The only difference between the two solutions is that one is a traditional OPF formulation with no withdrawal weight and the other is a water OPF formulation with a withdrawal weight of 0.1. Each solution produces a value of line loading for each line that we can express as a percentage of each line's limit. We visualize the difference in line loading between the traditional and water OPF as is done in figure 8. The convention for this figure is relative to the traditional OPF, such that negative changes mean that under the water OPF flows went down in that particular line. Busses are depicted as blue circles and transformers are depicted as black circles. All line flows values are tabulated in table SI2.
From figure 8, we see that by incorporating water into the OPF formulation, flows throughout the network change. We see that there are certain sections of the area that experience an increase in line loading (depicted as dark red), while other areas experience significant decreases in loading. We discuss how such information can be useful for system planning in section 5.
In figures 6, 7, and 8 we gain insights by subsetting and plotting different parts of the sampled input factor space. However, to facilitate interpretation of the entire sampled input factor space we rely on fitting unique CARTs to predict each of the four objectives based on the set of input factors. Several hyperparameters were specified to ensure these trees were properly pruned, such as having a maximum depth of five and maximum leaf nodes of twelve. We visualize the CARTs for the four objectives, shown in figures SI5-SI8 in the supplemental information.
From figures SI5-SI8, we gain insights into the general behavior of the system across all four input factors. Generally, not all the trees look the same. For this system, this dissimilarity means that each of the four objectives from our model are contributing different information. Figures SI8 and SI7 show that parts of the input factor space with the most severe system consumption and withdrawal are the result of significant splits between both the uniform water coefficient and uniform loading coefficients. Physically, this result means that both increased electricity demands and water inefficiencies must be present for this system to use the most water. Figure SI6, which depicts the splits based on the traditional generator cost objective, shows that increases to electricity demands more than water inefficiencies generally impact traditional generator costs (i.e., the most splits occur on the uniform loading coefficient). Figure SI5 shows that the total cost objective is generally dominated by the linear water terms as is evident by the general linear behavior. Further looking at the total cost tree in figure SI5, we see that increased electricity demand, decreased water efficiency, and a high withdrawal weight lead to the highest total cost. In other words, high withdrawal weights, uniform water factors, and uniform loading coefficients are all required for the system to experience the highest total costs.
Comparing the second and last rows of figure 7 we see the impact that generator water efficiency has on a system governed by our water-informed OPF. We see that the main differences come from plant 9 and plant 16. Under more severe water inefficiencies, the coal/once-through plant 9 reaches capacity before the recirculating/nuclear plant 16. However, what happens if local processes impact the water withdrawal rates of these plants? The non-uniform sensitivity analysis allows us to quantify the sensitivities of each individual plant's water use on the various aspects of system performance.

Impacts of non-uniform variability on system performance
In the previous section, we assumed that water changed uniformly (i.e., all generators experience the same increase/decrease) in the system. In this section, we investigate how allowing individual plant water usage to vary impacts system performance under a variety of operational scenarios. We apply the non-uniform sensitivity analysis outlined in section 2.4 to several operational scenarios. The operational scenarios as well as their mathematical representations are presented in table 3. We see that these scenarios capture different loading cases and policy cases. For loading cases, there are normal loading conditions (c load = 1.0) as well as increased scenarios (c load = 1.5). For policy cases, we consider the traditional OPF (which does not incorporate water use). We also consider the decreased withdrawal case, which puts a high optimization weight on the water withdrawn and low weight on consumption (e.g., policy focus on reducing withdrawal). The converse is true of the decreased consumption case. Recall that the choice of these values were based on both feasibility and historical validity as discussed in section 4.1.
Within each operational scenario, we test the sensitivity of the water use coefficient corresponding to each plant for this system. Our synthetic grid contains multiple generators of different fuel/cooling types located in the same plant. As previously discussed in section 2.4, such groups of generators are given separate plant designations. Therefore, there are a total of 17 unique plants. However, only 10 of those plants have cooling systems. We assign non-uniform water coefficients to each of these ten plants C water,plant9 , C water,plant6 , C water,plant12 , C water,plant15 , C water,plant1 , C water,plant10 , C water,plant3 , C water,plant17 , C water,plant13 , and C water,plant16 (see table SI1 for plant information). These coefficients are the input factors of the non-uniform sensitivity analysis.
Recall from section 2.4, the non-uniform sensitivity analysis used Sobol indices to quantify input factor sensitivity. However, many of the traditional methods for estimating Sobol indices are dependent on assumptions on the input factor distributions (Cousin et al 2019, Hart andGremaud 2018). However, from the visualization of the historic distributions of the 10 input factors for the regional dataset in figure SI9, we see that the input factors do not exhibit behavior of common parametric distributions. Fortunately, previous works have  . First order Sobol indices non-uniform sensitivity analysis. Each row is a different objective. Each column is a different input factor. Each panel column represents a grouping of plants based on their fuel and cooling system type as indicated on the top. Each panel row indicates an operational scenario as indicated by the labels on the right. Under traditional OPF formulations, the generator cost and total cost objectives do not consider the input factors, and are undefined in the first order sensitivities indices; thus, they are not colored.
proposed methods to estimate Sobol indices based on samples alone (Li and Mahadevan 2016). These samplebased methods allow us to bootstrap sample the historic input factors. The number of bootstrap samples was determined to be 22 528 based on n = N * (2D + 2) with D being value of ten for the number of input factors and N being 1024 based on literature recommendations Usher 2017, Saltelli 2004). This sampling is repeated for each of the four operational scenarios. This sampling was run on the RMACC Summit supercomputer across ten nodes with two CPUs (Intel Xeon E5-2680 v3, 2.50 GHz, 24 cores/node) with 4.84 GB RAM per core which took a wall time of approximately 30 min. We compute approximate first-order sensitivity indices (see equation (11)) for each of the four objectives using the methods proposed in Li and Mahadevan (2016), for each of the ten input factors, repeated for each of the four operational scenarios. We arrange them in a heatmap as is visualized in figure 9. Each column in the figure represents the non uniform water coefficient associated with the plant labeled on the bottom. Each row indicates a unique objective labeled on the left. Each panel column represents a grouping of plants based on their fuel and cooling system type as indicated on the top. Each panel row indicates an operational scenario as indicated by the labels on the right. The color of each cell represents the first order sensitivity of the input factor on the objective.
From figure 9, we see which objectives are most sensitive to which plant's water efficiency. For example, in the panel on the second row in the first column we see that system-wide water consumption is sensitive to the water efficiency of plant 12 (depicted as bright yellow). Furthermore, we see that plants 10, 12, and 16 are almost always sensitive for all objectives and operational scenarios. This result makes sense given that these are the plants with the highest capacities (table SI1). We also see that plant 3 becomes more important during increased loading scenarios, confirming the insights about plant 3 found from the uniform sensitivity analysis. We discuss how these insights can inform operations and planning in the following section.

Discussion
Examples of extreme climate events that have impacted thermoelectric operations are common throughout the early 21st century. Perhaps the largest and most catastrophic event was the 2003 heat wave in France that forced a shutdown of 25 percent of its nuclear generators due to high cooling water inlet temperatures (Harto et al 2012, Poumadère et al 2005. Less severe events include the droughts in 2006 that required shutdown and curtailment of a number of generators in Europe (Union of Concerned Scientists 2007), as well as in Illinois and Minnesota, United States (Harto et al 2012). The following year, a unit at Browns Ferry nuclear plant was forced to shutdown and the other two units reduced capacity during a heat wave and record-setting demand (Union of Concerned Scientists 2007). Such reduced capacity operations cost the grid operators an estimated additional $1 million per day to purchase power from neighboring systems (Freedman 2012). Similarly, droughts in 2011 caused reductions in power output in Texas (Galbraith 2011, Reuters Staff 2011, with similar reductions occurring in the Northeastern and Midwestern United States the following year (Eaton 2012). Recent studies estimate generator reductions/shutdowns during the 2013-2016 droughts in India prevented 30 TWh of potential electricity generation (Luo 2017).
One commonality among all these events is that they marked instances where system operators needed to manually intervene in the traditional power system operations to accommodate external hydrologic, climatic, or demand factors. Instead of such reactive measures, our proposed water-informed OPF operational scenario framework provides an opportunity to take proactive measures. We discuss how the insights gained through the operational scenario framework are useful to system analysts and policy makers. This section demonstrates how these insights would be helpful had this operational scenario analysis been applied to an actual system.
Through the uniform sensitivity analysis, we showed that our water-informed OPF formulation can effectively mitigate impacts to water consumption and withdrawal from changes in system water inefficiencies or increased loading. These water-informed operational policies are especially useful if systems are experiencing sudden water constraints as these operational policies would be much quicker to deploy than more long-term solutions like a water pipeline or upgrades to plant components (Pacsi et al 2013), and they could avoid ecosystem damage from widespread granting of thermal variances (Lubega and Stillwell 2018, Micha 2014, Logan et al 2021. The proposed water-informed power system operations compliments the efforts to include water in long-term planning operations (Jornada and Leon 2016). These efforts are also analogous to the efforts to incorporate air pollutants into both long-term grid planning and short-term grid operations (Kumar et al 2020, Flores and Brouwer 2018, Ren et al 2010, both of which generally fall under the need to create a grid that meets environmental targets (Ipakchi and Albuyeh 2009).
We demonstrated insights on the system-wide, plant/generator, and transmission system levels in the uniform sensitivity analysis. The system-wide insights through the CART analysis demonstrated how analysts can quickly gain insights into the relative importance of water inefficiencies or increased loading on the various aspects of system performance. The CARTs divide the possible operational space into well-defined regions. Such regions could be used as signposts for extreme system events to allow analysts to quickly categorize current system operations relative to possible extreme events. Such signposts could also be used to create adaptable policy where different actions are taken as the system experiences differing levels and types of stress (Herman and Giuliani 2018).
Our framework is also helpful for determining what actions system operators or policy makers could implement by showing which individual plants/generators are impacted by various external stressors. For example, targeted maintenance or generator-level operational policies could be put in place so that critical generators are ready to operate when needed. Similarly, through analysis of line flows, analysts can identify areas of the transmission system that would be most impacted by water-informed operational policies. Because increasing line loading closer to capacity can result in increased losses, such parts of the transmission system could be preemptively upgraded. Such insights address broader issues of creating policy that integrates day-to-day component-level operations to long-term system planning and behavior (Battey et al 2014).
Non-uniform loading allows analysts to see which plants most impact various aspects of system performance. For example, our analysis showed that plant sensitivity to generator cost, water consumption, and water withdrawal were often different. Such insights allow for the creation of policy that directly targets the most important plants, or types of plants, with respect to a policy maker's values. For example, our framework could inform a policy that targets water consumption under increased loading, which might be different from a policy that targets water withdrawal under normal loading.

Conclusion
This analysis contributes and demonstrates a novel operational scenario analysis framework to analyze power system performance under loading and water stresses. This framework is built on a novel water-informed formulation of the OPF problem whose solutions are a function of generator water efficiency, system loadings, and water use penalty terms. A multi-objective model built on this formulation allows analysts insights into relevant aspects of system performance. Sensitivity analyses of this multi-objective model facilitates both system-level and component-level insights; such insights are relevant to both system analysts as well as policy makers.
This operational scenario analysis framework also has possibilities to be expanded in future work. One extension would be to more explicitly model additional power systems decisions and hydrologic processes that change over various timescales (Oikonomou et al 2021). For example, generator ramping constraints, cooling system idling, flows of water throughout the river system, and market conditions are all processes that would likely impact subsequent timesteps when solving a multi-period OPF (Tidwell et al 2019, Rosenkranz et al 2016. One possible way to extend the water-informed OPF formulation proposed in this paper to a multi-period OPF model, would be to link the exogenous parameters in each timestep via external models. For example, the exogenous generator water efficiency parameters could be linked via a hydrologic/hydraulic model and exogenous load parameters could be linked via a system loading/climate model (Miara et al 2018, Li et al 2018. Pairing our formulation with physical-based models could allow analysts to better understand the physical causes of external power system stresses. Such information could be helpful for creating warning systems for system failure. Possibilities for future implementation of a water-informed OPF is another extension of this work. As has been noted in previous analyses to incorporate water into power systems operations, there still exists the fundamental 'who pays' question in regard to how assigning water a value could be incorporated into various energy market structures (Sanders et al 2014). One possible solution would be to follow the emerging market structures for pricing carbon (Narassimhan et al 2018). Additionally, the multi-objective water-informed OPF model could be optimized to find Pareto-optimal sets of operating policies and optimization weights. Participatory frameworks could also be implemented to support decision makers in the process of selecting such operating policies (Smith et al 2017).